Calculation method of carbon flow distribution in load-intensive regional energy centers

Abstract

With the development of the economy, people’s demand for green energy has increased significantly. However, the traditional single fossil energy supply system cannot meet the needs of low-carbon. Therefore, this study employs energy hub to establish a multi-energy flow network that enables the integration of carbon flow within the network. Additionally, by utilizing the multi-energy flow trend, a carbon flow tracking method is adopted to achieve real-time carbon flow calculation. Results show that this network calculates the electricity cost of 20043 yuan, gas cost of 67253 yuan, and carbon emission cost of 3152 yuan. Compared with the traditional energy flow system, gas cost is reduced by 4.3% and 1.7%, electricity cost by 21.3% and 15.0%, and carbon emission cost by 8.7% and 6.6%. The two-way sharing carbon flow calculation model calculates that the user side and power supply side of the node each bear half of the network loss, proving two-way sharing effectiveness. Test results on IEEE5 machine 14-node system show that the calculation method can accurately find high-emission and low-emission areas, making the carbon emission allocation between power generation and user more fair and reasonable. This research can effectively reduce emissions cost, accurately calculate emissions flow in real time, and facilitate reasonable emission reduction planning.

1. Introduction

The multi-energy flow network is the various energy sources interconnection network in energy center. In this network, different types of energy sources are supplied, transmitted, and utilized through various energy devices [1]. In the calculation of carbon flow (CF), it is necessary to evaluate and calculate the carbon emission (CE) of energy in this network. For different energy types, its carbon emission factor needs to be considered. By conducting comprehensive calculations of energy consumption and carbon emission factors, it is possible to determine the carbon emissions associated with different energy sources. By comprehensively considering the carbon emissions and the carbon flow between energy sources in the multi-energy flow power flow network, carbon emission reduction potential of energy centers can be evaluated and calculated more comprehensively. By analyzing the path and characteristics of carbon flow, we can equationte corresponding carbon emission reduction strategies, optimize energy flow and carbon emissions, and achieve low-carbon development and sustainable energy utilization [2]. Therefore, based on the energy center multi-energy flow network (MEFN) and carbon flow mapping relationship model, this study constructs a calculation method that can be used for real-time calculation of carbon flow distribution and couples with the above mapping model to realize real-time calculation and analysis of energy center carbon flow distribution. In addition, the current combination of multi energy flow network and carbon emission calculation is relatively rare, and it focuses on modeling and optimizing the operation of multi energy flow systems, as well as overall analysis of the multi energy flow system. There is a lack of overall research on the operating characteristics of other aspects of the system, such as renewable energy and residual energy circulation. Therefore, the carbon and sulfur distribution calculation method proposed in the study is innovative. It innovatively expands the research on the relationship between multi energy flow networks and carbon flow mapping, and provides theoretical support for the calculation method of carbon flow distribution in load intensive regional energy centers. This study is divided into four parts. First is the research findings regarding multi-energy flow networks and the calculation of carbon flow. Second is the construction of multi-energy and power flow network coupling model and carbon flow calculation model. The third part tests and analyzes the model, and the fourth part summarizes the article and points out the shortcomings.

2. Related works

The multi-energy flow network can realize the coupling between different energies to meet different energy demands. Some experts and scholars have related research based on this. Xie et al. addressed the related issues of energy coupling in multi energy flow power flow networks in power transportation, and constructed a corresponding coupling network based on considering the node electricity price of road charging stations, thereby effectively improving the efficiency of energy transfer [3]. Ding et al. comprehensively discussed the coupling method of energy hubs based on the consideration of multi energy flow networks to address the issue of conversion efficiency of different renewable energy sources, thereby reducing carbon emissions while strengthening energy control [4]. Yang et al. addressed the relevant issues of promoting energy marketization in multi energy coupling and proposed an indirect multi energy trading method based on the multi energy flow network, effectively improving energy utilization efficiency [5]. Wang et al. proposed an architecture for energy planning in the power system by considering the high penetration rate of renewable energy and considering the coupling of multiple energy sources, which effectively improves the reliability and safety of the energy system [6].

Jiang et al. addressed the relevant issues in power control in the Pengdian network and comprehensively discussed energy conversion and industrial time in distribution networks based on multi energy flow networks, providing data support for improving the efficiency and flexibility of energy utilization in distribution networks [7]. Wang et al. proposed a multitask architecture model based on the analysis of energy coupling relationships in integrated energy systems to address related issues in multi energy load forecasting, thereby improving the accuracy of multi energy load forecasting [8]. Yang et al. found that existing carbon allocation methods are mainly based on the principle of maximum efficiency or the principle of constant efficiency, however, these two principles are extreme cases and are particularly difficult to achieve in reality. Based on this, an alternative approach based on the principles of step-by-step efficiency planning and emission reduction planning is proposed. Applying the proposed method to the allocation of carbon emission quotas in various provinces in my country, some conclusions with practical guiding significance have been drawn [9].

To sum up, reducing carbon emissions has become an important research object in the field of environmental protection, and MEFN can realize coordination and dispatching of various energy sources. However, existing methods seldom combine the multi-energy flow network with the calculation of CE flow. Therefore, the CE calculation method using multi-energy flow constructed in this study has important implications in the fields of energy and environmental protection. D.

3. Carbon flow distribution model construction in load-intensive areas

To calculate CE in load-intensive areas and realize the requirements of emission reduction, this chapter is divided into two parts for model construction. The first part is optimized using MEFN to realize the coupling between carbon flow and MEFN. The second part realizes the real-time calculation of carbon flow using the carbon flow tracking method of multi-energy flow.

3.1 MEFN optimization model construction

MEFN is the interconnection of various energy sources in the energy center. In MEFN, different types of energy are supplied, transmitted and utilized through various energy devices [10]. Its system operation mainly depends on the coupling between various energy sources, which is usually described by an energy hub. Its schematic system diagram is shown Fig. 1.

Figure 1.

MEFN system schematic diagram.

In the energy hub, conversion, transmission, and storage among various energy forms mainly rely on the coupling matrix to express [11], as shown in Eq. (1).

$\{\begin{array}{c}L=CP\hfill \\ \left[\begin{array}{cccccccccccccccccccc}\hfill L_1\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill L_2\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill L_n\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]=\left[\begin{array}{cccccccccccccccccccc}\hfill c_{11}\hfill & \hfill c_{12}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill c_{1n}\hfill & & & & & & & & & & & & & & & & \\ \hfill c_{21}\hfill & \hfill c_{22}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill c_{2n}\hfill & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & & \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & \\ \hfill c_{n1}\hfill & \hfill c_{n1}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill c_{nn}\hfill & & & & & & & & & & & & & & & & \end{array}\right]\left[\begin{array}{cccccccccccccccccccc}\hfill P_1\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill P_2\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill P_n\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]\hfill \end{array}$(1)

In Eq. (1),

$\left[\begin{array}{cccccccccccccccccccc}\hfill L_1\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill L_2\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill L_n\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]$

represents the output variables of various energy components, and

$\left[\begin{array}{cccccccccccccccccccc}\hfill c_{11}\hfill & \hfill c_{12}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill c_{1n}\hfill & & & & & & & & & & & & & & & & \\ \hfill c_{21}\hfill & \hfill c_{22}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill c_{2n}\hfill & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & & \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & \\ \hfill c_{n1}\hfill & \hfill c_{n1}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill c_{nn}\hfill & & & & & & & & & & & & & & & & \end{array}\right]$

represents the coupling matrix of energy, and the elements in it $c_{ij}$ represent coupling conversion coefficients, which are expressed as output energy ratio of first $i$ to first $j$ . The coupling matrix describes the energy distribution and conversion relationship of various energy sources during inputting various energy sources to the load output port [12]. All kinds of energy can be converted into energy according to a certain ratio through energy distribution, and energy conversion can convert a certain type of energy into another type of energy through certain equipment [13]. Figure 2 shows the energy hub.

Figure 2.

Energy hub schematic diagram.

The calculation equation of energy conversion efficiency is shown in Eq. (2)

$C=\eta N$(2)

In Eq. (2), $N$ represents the energy distribution matrix $N=(v_{ij})$ , $v_{ij}$ represents the energy distribution coefficient, $\eta$ represents the conversion efficiency matrix $\eta =({\eta }_{ij})$ , ${\eta }_{ij}$ represents the energy conversion coefficient. In practical applications, there are various nonlinear elements in the energy matrix, so corresponding factors can be added according to practical applications [14]. In MEFN, power, natural gas, and thermal systems are the three core systems. The node calculation method constructs the power system flow equation [15], as shown in Eq. (3).

$\{\begin{array}{c}{P}_{ij}=(g_{si}+g_{ij})V_i^2-g_{ij}{V}_i{V}_j\cos {\theta }_{ij}-b_{ij}{V}_i{V}_j\sin {\theta }_{ij}\hfill \\ Q_{ij}=-(b_{si}+b_{ij})V_i^2+b_{ii}{V}_i{W}_j\cos {\theta }_{ij}-g_{ij}{V}_i{V}_j\sin {\theta }_{ij}\hfill \end{array}$(3)

In Eq. (3), $P_{ij}$ represents the active power of the transmission branch $ij$ , and $Q_{ij}$ represents the reactive power of the transmission branch. Therefore, the active power and reactive power of the total branch of the circuit can be obtained, as shown in Eq. (4).

$\{\begin{array}{c}{P}_i={\displaystyle \sum _{j\in i}}{P}_{ij}\hfill \\ Q_i={\displaystyle \sum _{j\in i}}{Q}_{ij}\hfill \end{array}$(4)

In Eq. (4), $P_i$ represents the active power injected by each node, and $Q_i$ represents the reactive power injected by each node. The natural gas network is controlled by fluid mechanics, and the air pressure of the nodes is the main variable. Therefore, natural gas system power flow equation is shown in Eq. (5).

$Q_i-{\displaystyle \sum _{m\in i}}{f}_{im}+{\displaystyle \sum _{n\in i}}{f}_{ni}-{\displaystyle \sum _{n\in i}}{G}_{ij}{F}_j=0$(5)

In Eq. (5), $Q_i$ is the net injection amount of natural gas at the node $i$ , $f_{im}$ and $f_{ni}$ represent the flow injected from the upstream node and output to the downstream node, respectively, $F_j$ represents compressor gas consumption, and $G_{ij}$ is the correlation coefficient when gas is taken at the node $i$ . The heating network includes hydraulic networks, whose calculation equation is shown in Eq. (6).

$\{\begin{array}{c}{A}_{s}m=m_q\hfill \\ B_h{h}_f=0\hfill \\ h_f=Km\left|m\right|\hfill \end{array}$(6)

In Eq. (6), $A_s$ it represents the node-branch correlation matrix, $m$ represents pipeline flow, $m_q$ represents the outflow flow, $B_h$ represents the loop-branch correlation matrix, and $h_f$ represents the pressure head loss vector, $K$ represents the pipeline resistance coefficient matrix. The thermal network satisfies the law of thermodynamics, and the calculation equation is shown in Eq. (7).

$\{\begin{array}{c}\mathrm{\Phi }=C_p{m}_q(T_s-T_o)\hfill \\ T_{\text{𝑒𝑛𝑑}}=(T_{\text{𝑠𝑡𝑎𝑟𝑡}}-T_a)e^{-\frac{\lambda L}{C_{p}m}}+T_a\hfill \\ ({\displaystyle \sum m_{\text{𝑜𝑢𝑡}}})T_{\text{𝑜𝑢𝑡}}={\displaystyle \sum m_{in}{T}_{in}}\hfill \end{array}$(7)

In Eq. (7), $\mathrm{\Phi }$ it represents the nodes thermal power, $C_p$ represents the specific water heat capacity, $T_{\text{𝑠𝑡𝑎𝑟𝑡}}$ and $T_{\text{𝑒𝑛𝑑}}$ represent the temperature at the beginning and end of the pipe, represents $T_a$ environment temperature, $\lambda$ represents the pipe heat transfer coefficient, and $L$ represents pipe length, $m_{\text{𝑜𝑢𝑡}}$ and $T_{\text{𝑜𝑢𝑡}}$ , $m_{in}$ , and $T_{in}$ represent the outflow and inflow The pipe water flow rate and temperature. To apply the coupled multi-energy system in practice, the above power flow is expressed by algebraic equation, as shown in Eq. (8).

$\{\begin{array}{c}0=F(x_e,x_g,x_h,x_{eh})\hfill \\ 0=G(x_e,x_g,x_h,x_{eh})\hfill \\ 0=H(x_e,x_g,x_h,x_{ch})\hfill \\ 0=C_{EH}(x_e,x_g,x_h,x_{ch})\hfill \end{array}$(8)

In Eq. (8), $F$ , $G$ , $H$ represent the energy balance equations of three systems, respectively, and $C_{EH}$ represent the input-output coupling equations of the energy hub, $x_e$ and, $x_g$ , $x_h$ , $x_{eh}$ represent the relevant variables of the electricity, gas, heat energy systems, and the energy hub, respectively. Considering that there is energy loss in the energy system, all losses are converted into energy resistance coefficients, and a multi-energy flow network optimization model is constructed, as shown in Eq. (9).

$\{\begin{array}{c}f=\min C_{\text{𝑡𝑜𝑡𝑎𝑙}}\hfill \\ C_{in}N{P}_{in}=C_{\text{𝑜𝑢𝑡}}{P}_{\text{𝑜𝑢𝑡}}\hfill \\ 0⩽P⩽P_{\max }\hfill \end{array}$(9)

In Eq. (9), $P_{in}$ represents the input pipeline matrix, $P_{\text{𝑜𝑢𝑡}}$ represents the output pipeline matrix, $C_{in}$ represents the correlation matrix, $N$ represents the energy resistance coefficient matrix, and $P_{\max }$ is the maximum capacity of the energy storage system equipment in the system. Since each node in the system must satisfy the input and output phase balance [16]. In addition, the energy storage nodes included in the system must also contain constraints, as shown in Eq. (10).

$\{\begin{array}{c}0⩽P_{e,es}⩽u_c{\gamma }_{c,\max }{P}_{e,es,\max }\hfill \\ E_{k,0}=E_{k,T}\hfill \\ E_{s,t+1}=E_{s,t}+\left[{\eta }_c{P}_{e,i}(t)-{\displaystyle \frac{P_{d,i}(t)}{{\eta }_d}}\right]\mathrm{\Delta }t\hfill \\ 0⩽P_{d,es}⩽u_d{\gamma }_{d,\max }{P}_{d,es,\max }\hfill \\ u_c+u_d=1\hfill \end{array}$(10)

In Eq. (10), $E_{s,t}$ it is the energy stored at a certain time, ${\eta }_c$ represents the charging device efficiency, $P_{e,i}(t)$ is the charging device power at a certain time, ${\eta }_d$ is the energy device discharge efficiency, and $P_{d,i}(t)$ is the energy storage device’s Discharging power $E_{s,t+1}$ represents the remaining device energy running to the next time step, ${\gamma }_{c,\max }$ represents ${\gamma }_{d,\max }$ the charging and discharging efficiency of the energy storage device, $P_{e,es,\max }$ and, $P_{d,es,\max }$ represents the maximum charging and discharging device power. In order to minimize the system operating cost, an optimization model is constructed considering the cooling, heating, electrical loads and related equipment constraints, as shown in Eq. (11) [17].

$\{\begin{array}{c}f=\min C_{\text{𝑡𝑜𝑡𝑎𝑙}}=C_e+C_g+C_c\hfill \\ C_g={\displaystyle \sum _{i=1}^T}{p}_g(P_{mt,t}+H_{gb,t})\cdot \mathrm{\Delta }t\hfill \\ C_e={\displaystyle \sum _{t=1}^T}{p}_{\text{𝑔𝑟𝑖𝑏}}{p}_{\text{𝑔𝑟𝑖𝑏},t}\mathrm{\Delta }t\hfill \\ C_c=\alpha {\displaystyle \sum _{t=1}^T}\left[{\beta }_{\text{𝑔𝑟𝑖𝑏}}{p}_{\text{𝑔𝑟𝑖𝑏},t}+{\beta }_{mt}{P}_{mt,t}+{\beta }_{gb}{H}_{gb,t}\right]\mathrm{\Delta }t\hfill \end{array}$(11)

In Eq. (11), $\alpha$ it is the carbon emission cost factor of 0.038 yuan/kg, ${\beta }_{\text{𝑔𝑟𝑖𝑏}}$ represents that of the power grid, which is 0.889 kg/(kW h), and ${\beta }_{mt}$ represents that of the gas turbine of 0.184 kg/(kW h), ${\beta }_{gb}$ is that of the gas boiler, which is 0.226 kg/(kW $\cdot$ h). Through the construction of the above optimization model, energy saving and emission reduction can be achieved. Figure 3 shows MEFN system coupling model flow chart.

Figure 3.

Flow chart of the coupling model for a MEFN system.

3.2 Construction of carbon flow calculation model based on MEFN coupling

Intensive regional energy centers refer to centralized energy supply and multi-energy complementary energy systems built in urban centers or dense urban areas. The purpose is to utilize energy efficiently, reduce its consumption and carbon Emissions [18]. Figure 4 shows the energy supply system frame diagram in the load-intensive area.

Figure 4.

Framework diagram of energy supply system in densely loaded areas.

Figure 5.

Multi energy complementary mode diagram.

In order to further realize the statistics of CE and reduce it, CF analysis method is used. This method utilizes the CF tracking method of the multi-energy flow trend, by tracking CE generated during the flow and conversion of energy in different links, as well as the carbon flow transfer between different energy sources, to comprehensively understand and evaluate energy. The carbon emissions of the system [19]. The multi-energy complementary mode diagram is shown in Fig. 5.

The carbon flow tracking method is divided into forward tracking and countercurrent tracking. In the forward tracking, the active power flow is set as the positive direction of the branch. Therefore, the power output of the node, branch inflow power, and the load contribution are as follows Eq. (12) shown.

$\{\begin{array}{c}{[A_d]}_{ik}=\{\begin{array}{c}1,i=j\hfill \\ -S_{ij}/S_j,j\in D_i\hfill \\ 0,\text{𝑜𝑡ℎ𝑒𝑟}\hfill \end{array}\hfill \\ S_{Gi}={\displaystyle \frac{S_{Gi}}{S_i}}{\displaystyle \sum _{k=1}^n}{[A_d^{-1}]}_{ik}{S}_{Lk}\hfill \\ S_{fi}={\displaystyle \frac{S_{fi}}{S_i}}{\displaystyle \sum _{k=1}^n}{[A_d^{-1}]}_{ik}{S}_{Lk}\hfill \end{array}$(12)

In Eq. (12), $S_{Gi}$ is $i$ ’s power output, $S_{fi}$ represents the input power of the branch, $S_{Lk}$ represents the node load, $S_i$ represents the node flowing power, $A_d$ represents $(n+m)\times (n+m)$ the order downstream distribution matrix, and represents $D_i$ the set of downstream nodes of the node $i$ . Therefore, the negative value of the branch loss is regarded as the virtual power supply output, and network loss borne by the load is shown in Eq. (13).

$\mathrm{\Delta }{S}_{Lk}=-{\displaystyle \sum _{i=n+1}^{n+m}}{\displaystyle \frac{S_{Gi}}{S_i}}{[A_d^{-1}]}_{ik}{S}_{Lk}$(13)

In Eq. (13), $\mathrm{\Delta }{S}_{Lk}$ it represents the load network loss borne. In the upstream power flow tracing, the calculation of each parameter is shown in Eq. (14).

$\{\begin{array}{c}{S}_{Li}={\displaystyle \frac{S_{Li}}{S_i}}{\displaystyle \sum _{k=1}^n}{[A_U^{-1}]}_{ik}{S}_{Gk}\hfill \\ S_{it}={\displaystyle \frac{S_{it}}{S_i}}{\displaystyle \sum _{k=1}^n}{[A_u^{-1}]}_{ik}{S}_{Gk}\hfill \\ {[A_u]}_{ij}=\{\begin{array}{c}1,i=j\hfill \\ -S_{ij}/S_j,j\in U_i\hfill \\ 0,\text{𝑜𝑡ℎ𝑒𝑟}\hfill \end{array}\hfill \end{array}$(14)

In Eq. (14), $A_u$ represents $(n+m)\times (n+m)$ the order upstream distribution matrix and $U_i$ represents the set of upstream nodes of the node. Branch loss is regarded as a virtual charge, and network loss borne by power supply is shown in Eq. (15).

$\mathrm{\Delta }{S}_{Gk}={\displaystyle \sum _{i=n+1}^{n+m}}{\displaystyle \frac{S_{Li}}{S_i}}{[A_u^{-1}]}_{ik}{S}_{Gk}$(15)

In Eq. (15), $\mathrm{\Delta }{S}_{Gk}$ it represents network loss borne by power supply. Since only single-factor network loss is considered in both downstream and upstream power flow tracking, in order to realize network loss sharing, the two-way network loss method is used to calculate the network loss of load and power supply, as shown in Eq. (16).

$\{\begin{array}{c}\mathrm{\Delta }{S}_{Lk}=-{\displaystyle \sum _{i=n+1}^{n+m}}{\displaystyle \frac{\beta S_{Gi}}{S_i}}{[A_d^{-1}]}_{ik}{S}_{Lk}\hfill \\ \mathrm{\Delta }{S}_{Gk}={\displaystyle \sum _{i=n+1}^{n+m}}{\displaystyle \frac{(1-\beta )S_{Li}}{S_i}}{[A_u^{-1}]}_{ik}{S}_{Gk}\hfill \end{array}$(16)

In Eq. (16), $\beta$ represents an adjustable parameter, and its value is [0, 1]. After allocating network losses to loads and power sources, the original multi-energy flow network is transformed into a virtual lossless network. In this network, the branch loss is 0, so the active power distribution of each branch will be redistributed, as shown in Eq. (17).

$P^{\prime }=B^{\prime }\theta$(17)

In Eq. (17), $P^{\prime }$ it represents the vector formed by node intersection voltage in virtual network, $\theta$ represents the vector formed by node input power, and $B^{\prime }$ represents the node admittance matrix formed by branch reactance. According to Eq. (17), the parameter values in the upstream power flow tracking can be recalculated, as shown in Eq. (18).

$\{\begin{array}{c}{P}_{Li}^{\prime }={\displaystyle \frac{P_{Li}^{\prime }}{P_i^{\prime }}}{\displaystyle \sum _{k=1}^n}{[A_u^{{}^{\prime }-1}]}_{ik}{P}_{Gk}^{\prime }\hfill \\ P_{it}^{\prime }={\displaystyle \frac{P_{it}^{\prime }}{P_i^{\prime }}}{\displaystyle \sum _{k=1}^n}{[A_u^{{}^{\prime }-1}]}_{ik}{P}_{Gk}^{\prime }\hfill \end{array}$(18)

In Eq. (18), $P_{Li}^{\prime }$ represents the node load, $P_{it}^{\prime }$ represents the branch power, $P_{Gk}^{\prime }$ represents the node net active power output, and $A_u^{{}^{\prime }-1}$ represents the updated upstream distribution matrix. Carbon emissions mainly come from the gas emitted by thermal power generating units, which are affected by factors such as fuel quality and unit performance. Therefore, the coal consumption per degree of electricity of thermal power units is shown in Eq. (19).

$w_i=(a_i{P}_{Gi}^2+b_i{P}_{Gi}+c_i)\zeta /P_{Gi}\times {10}^3$(19)

In Eq. (19), $a_i$ , $b_i$ , $c_i$ represent the characteristic parameters of the coal consumption of the unit under normal operating conditions, respectively. $\zeta$ The correction parameter representing the coal consumption per degree of electricity is 1.0 under normal conditions, 0 under shutdown conditions, and 1.01 under load conditions. Unit CE intensity calculation is shown in Eq. (20).

$E_{Gi}={\eta }_i\times {\xi }_i\times (1-{\mu }_i)\times {\displaystyle \frac{M_{C{O}_2}}{M_C}}\times {\displaystyle \frac{w_i}{{10}^3}}$(20)

In Eq. (20), ${\eta }_i$ , ${\xi }_i$ , ${\mu }_i$ respectively represent the unit coal content, carbon oxidation rate and carbon capture rate, $M_C$ and $M_{C{O}_2}$ represent the molar mass of carbon and carbon dioxide, respectively, and the values are 12 g/mol and 44 g/mol. Combining Eqs (18)–(20), the day-ahead basic node load, branch power and unit carbon flow rate (CFR) are obtained, as shown in Eq. (21).

$\{\begin{array}{c}{R}_{Li}={\displaystyle \frac{P_{Li}^{\prime }}{P_i^{\prime }}}{\displaystyle \sum _{k=1}^n}{[A_u^{{}^{\prime }-1}]}_{ik}{P}_{Gk}^{\prime }{E}_{Gk}\hfill \\ R_{it}={\displaystyle \frac{P_{it}^{\prime }}{P_i^{\prime }}}{\displaystyle \sum _{k=1}^n}{[A_u^{{}^{\prime }-1}]}_{ik}{P}_{Gk}^{\prime }{E}_{Gk}\hfill \\ R_{Gk}=Re[\mathrm{\Delta }{S}_{Gk}]E_{Gk}\hfill \end{array}$(21)

In Eq. (21), $R_{Li}$ it represents the load CFR. $R_{it}$ Indicates the branch CFR, and $R_{Gk}$ indicates the CFR of the unit network loss borne, that is, $\mathrm{\Delta }{S}_{Gk}$ the real part after taking the complex number. Since the bidirectional carbon flow calculation method lacks consideration of fuel consumption in the preparation stage of thermal power units, the power deviation is considered in the countercurrent power flow tracking method, and the deviation equation is shown in Eq. (22).

$\{\begin{array}{c}{\delta }_{Gk}^{\prime }={\delta }_{Gk}\hfill \\ {\delta }_{Lk}^{\prime }={\delta }_{Lk}+\mathrm{\Delta }{\delta }_{Lk}\hfill \end{array}$(22)

In Eq. (22), ${\delta }_{Gk}^{\prime }$ is the unit actual power deviation, and ${\delta }_{Lk}^{\prime }$ represents the actual power deviation of the load. Therefore, CE intensity of mechanism power deviation is shown in Eq. (23).

$\mathrm{\Delta }{E}_{Gi}={\displaystyle \frac{{\stackrel{\sim }{E}}_{Gi}{\stackrel{\sim }{P}}_{Gi}-E_{Gi}{P}_{Gi}}{{\stackrel{\sim }{P}}_{Gi}-P_{Gi}}}$(23)

In Eq. (23), $P_{Gi}$ represents the planned output of the thermal power unit, ${\stackrel{\sim }{P}}_{Gi}$ represents the actual output, $E_{Gi}$ represents CE intensity corresponding to the planned output, and ${\stackrel{\sim }{E}}_{Gi}$ represents that corresponding to the actual output. Therefore, the modified carbon flow rate correction amount is shown in Eq. (24).

$\{\begin{array}{c}\mathrm{\Delta }{R}_{Li}={\displaystyle \frac{{\dot{P}}_{Li}^{\prime }}{{\dot{P}}_i^{\prime }}}{\displaystyle \sum _{k=1}^n}{[{\dot{A}}_u^{{}^{\prime }-1}]}_{ik}{\dot{P}}_{Gk}^{\prime }\mathrm{\Delta }{E}_{Gk}\hfill \\ \mathrm{\Delta }{R}_{it}={\displaystyle \frac{{\dot{P}}_{it}^{\prime }}{{\dot{P}}_i^{\prime }}}{\displaystyle \sum _{k=1}^n}{[{\dot{A}}_u^{{}^{\prime }-1}]}_{ik}{\dot{P}}_{Gk}^{\prime }\mathrm{\Delta }{E}_{Gk}\hfill \\ \mathrm{\Delta }{R}_{Gk}=0\hfill \end{array}$(24)

In Eq. (24), $\mathrm{\Delta }{R}_{Li}$ represents the load CFR correction amount, $\mathrm{\Delta }{R}_{it}$ represents the branch CFR correction amount, and $\mathrm{\Delta }{R}_{Gk}$ represents the CFR correction amount of unit, and a value of 0 means that the unit does not bear the CFR change caused by deviation. Add the correction amount and the CFR before the day to get the actual CFR as shown in Eq. (25).

$\{\begin{array}{c}{\stackrel{\sim }{R}}_{Li}=R_{Li}+\mathrm{\Delta }{R}_{Li}\hfill \\ {\stackrel{\sim }{R}}_{it}=R_{it}+\mathrm{\Delta }{R}_{it}\hfill \\ {\stackrel{\sim }{R}}_{Gk}=R_{Gk}+\mathrm{\Delta }{R}_{Gk}\hfill \end{array}$(25)

In Eq. (25), ${\stackrel{\sim }{R}}_{Li}$ , ${\stackrel{\sim }{R}}_{it}$ , ${\stackrel{\sim }{R}}_{Gk}$ respectively represent the load actual CFR, branch and unit that include the embodied carbon emissions. Through the calculation of the above equation, CFR in the unit can be calculated in real time. The flowchart is shown in Fig. 6.

Figure 6.

CF calculation flowchart.

4. Case analysis of MEFN optimization model and carbon flow calculation model

To test the emission reduction performance of MEFN optimization and the correctness of the carbon flow calculation model, this chapter is divided into two parts for testing. The first part tests and analyzes the multi-energy flow optimization model, and the second part tests and analyzes the carbon flow calculation model.

4.1 Practical case analysis of MEFN power flow optimization model

The industrial park of a certain city is selected as the research object, which includes a residential area with an area of 1.5 km², an industrial area with an area of 1.8 km² and a commercial area with an area of 0.8 km². To test the generality of MEFN optimization model, three different cases were selected for analysis. The first type adopts the traditional cogeneration system, and the second type introduces renewable energy systems, including photovoltaic and wind power generations, on basis of the first type. The third type introduces an energy storage system on basis of the second type. Assuming that the system is connected to the grid, the supply of internal power generation and renewable energy is insufficient to purchase electricity from the external grid. Meanwhile, there is no natural gas source in the park, which is supplied by the urban natural gas pipeline network. Related equipment parameters are shown below.

Table 1

Related equipment parameter settings table

Device	Rated power (MW/h)	Conversion efficiency/charging and discharging efficiency/refrigeration system
Photovoltaic power generation board	0.40	0.97
Wind turbine	0.50	0.96
Gas turbine	Electricity: 4.00/Heat: 6.00	Gas-Electricity: 0.30/Gas-heat: 0.45
Waste-heat recovery unit	4.00	0.80
Gas Boiler	3.00	0.88
Storage battery	1.80	Charge: 0.95/Discharge: 0.96
Heat storage tank	3.60	Thermal storage: 0.96/Exothermic: 0.97
Electric refrigerator	2.00	3.60
Absorption chiller	4.00	1.40

Assuming that the day-ahead power consumption of electric heating and cooling loads before optimization and the day-ahead prediction of photovoltaic, wind power, and electric heating and cooling loads are accurate, the day-ahead forecast power and day-ahead output power are demonstrated in Fig. 7.

Figure 7.

Daily predicted power and daily output power.

Figure 7a shows the predicted power of cooling, heating, and electric loads. It can be seen that the cooling load reaches a peak value of 4600 kW at 11:00, the heat load reaches a peak value of 4000 kW at 12:00, and the electric load reaches a peak value of 4000 kW at 17:00. When it reaches a peak value of 5100 kW. Figure 7b shows the daily output power of photovoltaic and wind power generation. Their peak values are 250 kW at 12:00 and 350 kW at 21:00, relatively. The multi-energy flow optimization model is divided into energy sources on each energy node, and the daily average value of a standard month is used as data. During peak hours, electricity price is 0.65 yuan/(kW h), and natural gas price is 0.36 yuan/(kW h). During off-peak hours, the price of electricity is 0.26 yuan/(kW h), and natural gas price is 0.36 yuan/(kW h). During normal hours, electricity price is 0.45 yuan/(kW h), and natural gas price is 0.36 yuan/(kW h). The results of electric power optimization in three different cases are shown in Fig. 8.

Figure 8.

Electric power scheduling results of multi energy flow systems in three scenarios.

Figure 8a shows electric power dispatching results of the traditional cogeneration system, which mainly uses the functions of the external power grid and natural gas pipeline network. It is calculated that the electricity purchase cost is 25,482 yuan, the gas purchase cost is 70,259 yuan, and the carbon emission cost is 3,455 yuan. Figure 8b introduces renewable energy, and the optimized power part decreases and the peak value is also flat. It is calculated that the electricity purchase cost is 23,592 yuan, the gas purchase cost is 68,429 yuan, and the carbon emission cost is 3,376 yuan. Figure 8c adopts a multi-energy flow optimization model, which reserves and transmits the idle power consumption to the peak power consumption stage, which increases the utilization efficiency of renewable energy and reduces carbon emissions. It is calculated that the electricity purchase cost is 20,043 yuan, the gas purchase cost is 67,253 yuan, and the carbon emission cost is 3,152 yuan. Compared with the first and second cases, the gas cost is reduced by 4.3% and 1.7%, electricity cost by 21.3% and 15.0%, and carbon emission cost by 8.7% and 6.6%.

4.2 Verification of carbon flow calculation model

A node system is selected to verify the carbon flow calculation model. A total of 4 nodes are set in the node system, and the units connected to each node are different. Node 1 is a balance node, the coal consumption characteristic parameters are set to 0.0014, 0.12 and 0.25, the voltage amplitude is 1.050pu, the correction parameter is 1.01, the carbon content rate is 0.80, and the carbon oxidation rate is 0.98 $\zeta$ . The carbon capture rate is 0.80. Node 2 is connected to a wind farm with a voltage amplitude of 1.045pu and a planned output of 114MW. Both other nodes are load nodes, and day-ahead load prediction values are 200 $+$ j80MVA and 300 $+$ j120MVA. The multi-energy flow calculation method can calculate system’s power flow, and the upstream and downstream tracking algorithms are used to allocate network losses to loads and units. The results of different methods are shown in Fig. 9.

Figure 9.

Unit, load, and branch allocation.

From Fig. 9, units-shared total load is 515MW. When the network loss is allocated to user, the load allocation increases to 536MW, and G1 and G2 remain unchanged. When the network loss is allocated to the power generation, the units G1 and G2 bear the active power loss of 18MW and 3MW. When two-way sharing is used, the active power losses borne by loads L1 and L2 are 3.2MW and 8.1MW respectively, and the active power losses shared by units G1 and G2 are 9.4MW and 1.9MW respectively. The user and power sides each bear half of the network loss, which proves the effectiveness of two-way sharing. Set the predicted load carbon flow rate L1 to 34 t/h, L2 to 66 t/h, unit carbon flow rate G1 injection to 82 t/h, G1 network loss to 1.5 t/h. The current actual carbon flow rate is shown in Fig. 10.

Figure 10.

Actual CFR before the day.

Figure 10, static analysis method 1 represents network loss distribution to the load, and the static analysis method 2 represents the distribution of network loss to the power supply. Static analysis method 1 is compared with static analysis method 2, the G1 injection of static analysis method 1 is 76 t/h, and the network loss is 0. The G1 injection of static analysis method 2 is 73 t/h, the network loss is 3.6 t/h, and CFR has increased. The incremental sum of load actual CFR and predicted CFR shows that the change in system carbon emissions is shared by unit and load. To further prove method’s effectiveness, a test is carried out on an IEEE5 machine 14-node system, and the results are shown in Fig. 11.

Figure 11.

Load CFR and unit injection CFR.

Figure 12.

CFR and active power flow of each branch.

Figure 11 is the CFR of each node calculated from the power drawn by the load node. CE generated by the network loss is attributed to the power supply side, and the load side does not bear this part of CE. The total load CFR of grid loss is 138.5 t/h, and the total CFR of unit injection is 139.8 t/h. Therefore, CFR generated at the grid loss site is 1.4 t/h. CFR and active power flow of each branch are shown in Fig. 12.

From Fig. 12 that generator set’s CE intensity at node 1 is relatively large, and the generated electric energy is also more, which belongs to the area with heavy carbon emission. Therefore, the 16 branches and 12 branches connected to node 1 are the power to withstand is also greater, and CFR is also greater. The branch 23 bears CE caused by the transmission power of the generator sets 1 and 2, so branch CFR is also relatively large. Since CE from the node 5 generator is 0, CFR is also 0. This calculation method can accurately find high-emission and low-emission areas, and provide a basis for subsequent emission reduction planning.

5. Conclusion

To realize the real-time calculation of carbon flow in the load energy center, this research is based on the MEFN optimization to realize the coupling between carbon flow and MEFN, and achieve the purpose of real-time calculation of carbon flow. The results show that the traditional cogeneration system calculates the electricity purchase cost of 25,482 yuan, the gas purchase cost of 70,259 yuan, and the carbon emission cost of 3,455 yuan. After the introduction of renewable energy, the optimized power part has declined, and the peak value is relatively flat. The calculated electricity purchase cost is 23,592 yuan, the gas purchase cost is 68,429 yuan, and the carbon emission cost is 3,376 yuan. The multi-energy flow optimization model calculates that the electricity purchase cost is 20,043 yuan, the gas purchase cost is 67,253 yuan, and the carbon emission cost is 3,152 yuan. Compared with the former two, gas cost has been reduced by 4.3% and 1.7%, electricity cost has been reduced by 21.3% and 15.0%, and CE cost was reduced by 8.7% and 6.6%. Reserve and transmit the idle power consumption to the peak power consumption stage, which increases the utilization efficiency of renewable energy and reduces CE. The test results on the IEEE5 machine 14-node system show that the calculation method can accurately find high-emission and low-emission areas, and find the key points in the low-carbon link, making CE sharing between the power generation side and the user side fairer and more reasonable. There are deficiencies in this study. In MEFN system construction, the complex disturbance at the source was not considered. In the next study, the relationship between the system load and photovoltaic and wind power generation should be considered.