Corresponding to the cloud computing center in the IES, the ECU proposed in this article is the local energy data management unit, which can reduce cloud computing burden by finishing a part of the computing tasks at the local (Weisong et al., 2017). The structure of the ECU model presented in this article is shown in Figure 1.

EC, HP, CHP, and GB, as the coupling equipment convert power and natural gas into cold, heating, and power. Cooling and power are directly supplied to the local load, while the heating is transmitted to the HS or the DHN to jointly support the heating demand of the local load.

The modeling method of the upper radiant heat network is similar to that of the power grid and natural gas network. Therefore, the ring heating network is only established at the lower layer in this study for heat circulation. The ring network has high security and reliability (Wang et al., 2016), and the computing method for the ring network also meets the requirements of the radial heating network.

The ECU supports the multizone expansion, where the units are connected to each other through a power grid, natural gas network, and DHN, becoming the main body of the MRIES.

The ECU obtains the data of the devices and nodes through the local awareness layer and submits the boundary information to the cloud after internal optimization. Then, the cloud feeds back the consensus information obtained from the centralized computing to each ECU. After that, the ECU adjusts the internal operation according to the consensus information. The MRIES achieves global optimization through several iterations.

The centralized optimization needs to submit all the device data to the cloud, which results in vast costs of communication facility construction, data transmission, storage, and processing. It deviates from the energy development idea of energy conservation and emission reduction. Meanwhile, the information security and privacy cannot be guaranteed, as shown in Figure 2A.

Figure 2B shows the distributed optimization based on the energy system decoupled. Each computing unit needs to process the data of the whole system, which is the mainstream distributed computing method at present. Although the amount of cloud workload is reduced by two to three times, it is still affected by the wide geographical range of the MRIES. Meanwhile, the cost of the communication facility construction is still relatively high. In addition, this method deconstructs the EH model, resulting in the separate operation of each coupling device. As the number and types of the coupling devices increase, the computational complexity of this method will increase and the cost of data transmission and storage will also increase significantly.

The decoupled method based on the ECU is shown in Figure 2A, where the data are transmitted within the RIES and only submits the boundary information at the coupling point between the ECUs to the cloud. This decoupled method dramatically reduces the cost of communication facility construction, data transmission, storage, and processing, while ensuring the security and privacy of the data. Furthermore, the ECU collects the data of three types of energy and processes them in a centralized manner, which can better play a synergistic role in optimization. In addition, this method lays a foundation for future research work of the unified energy management.

As previously mentioned, EH is the main body and the energy management unit of the RIES. Therefore, the EH serves as the main body of the ECU to decouple the energy system. The partition principle should be based on the whole system’s lowest communication cost. The communication cost is proportional to the physical distance between the nodes. Therefore, we decouple the energy networks separately according to the adaptive method of P-median problem as shown in Eq. 1. min ∑ r ∈ R ∑ f ∈ F l r f y r f  s.t.  ∑ r ∈ R y r f = 1 , ∀ f ∈ F y r f ∈ 0,1 , ∀ r ∈ R , ∀ f ∈ F . (1)

The objective function minimizes the information transmission distance within the network partitions, while the constraints ensure that the nodes are completely partitioned without repetition and omission. The IES model contains an electrical network, natural gas network, and DHN in the study. Therefore, the objective function should be applied in the three networks. After that, we can get the indices of the nodes in each unit of every network, but we need to further determine the boundaries between the units to provide the boundary data for cloud computing. As a result, section 3.2 provides a decoupled method to distinguish the boundaries between the units under different EH connection conditions.

As previously mentioned, the units contain different nodes but are connected by an energy transmission line. The transmission line also falls within the scope of unit internal optimization, so the virtual node needs to be inserted into the transmission line to establish clear boundaries between the units and serve as a data collection point of the connecting transmission line.

This study proposed a decoupled method, where the EH serves as the core in each ECU. Based on the connected position of the EH to the energy network, the whole system is decoupled by the partition method as given below. In Figure 3, the source represents the source of the energy network, which is the generator in the electrical network and the compressor in the natural gas network. EH represents a real hub. The virtual node is located in the middle of the transmission line, and the node data are the energy data in the middle of the transmission line. The virtual transmission line is lossless, aiming to reduce the coupling degree between the nodes.

In this study, the Distflow model is used to simulate the power flow of the distribution network in the ECU (Baran and Wu 1989), and the SOCP is used to deal with the nonlinear model. The power balance of the node described in Eqs 2–4 is the voltage loss equation, and Eq. 5 represents the relationship between the branch current, node voltage, and power. P j , t in = − ∑ g ∈ Π j P g , t G − ∑ p v ∈ Π j P p v , t PV − ∑ w t ∈ Π j P w t , t WT − P i j , t − I ̂ i j , t R i j + ∑ o ∈ Ω j P j o , t + ∑ d ∈ Π j P d , t LD + P j , t EH . (2) Q j , t in = − Q i j , t − I ̂ i j , t X i j + ∑ o ∈ Ω j Q j o , t + ∑ d ∈ Π j Q d , t LD . (3) V ̂ j , t = V ̂ i , t − 2 P i j , t R i j + Q i j , t X i j + I ̂ i j , t R i j 2 + X i j 2 . (4) 2 P i j , t 2 + 2 Q i j , t 2 + I ̂ i j , t − V ̂ i , t 2 ≤ I ̂ i j , t + V ̂ i , t 2 . (5)

The capacity limitation of the transmission line is shown as Eq. 6. I ̂ i j , t min ≤ I ̂ i j , t ≤ I ̂ i j max V ̂ i min ≤ V ̂ i , t ≤ V ̂ i max . (6)

Power consumption of the ECU is shown as Eq. 7. P n , t ECU = ∑ i ∈ I n upper P i , t in − ∑ j ∈ J n under P j , t in . (7)

The natural gas network model in the ECU is similar to the electrical network, including the node energy balance constraint and pipeline transport constraint. Eq. 8 represents the balance between the inflow and outflow of node b. G b , t in = − G a b , t + ∑ ω ∈ Λ b G b ω , t + ∑ d ∈ Υ b G d , t LD + G b , t EH . (8)

The natural gas network in this article is a medium–low pressure network without considering the compressor model (Hu et al., 2020). The second-order cone programming is used to relax the Weymouth equation of the natural gas network, as shown in Eq. 9 (Liu et al., 2020). Eqs 10,11 are the upper and lower limits of the gas transmission volume and natural gas pressure of the pipeline, respectively. G a b , t 2 + K a b π b , t 2 ≤ K a b π a , t 2 . (9) 0 ≤ G a b , t ≤ G a b max . (10) π a min ≤ π a , t ≤ π a max . (11)

The natural gas consumption of the ECU is shown as Eq. 12. G n , t ECU = ∑ a ∈ A n upper G a , t in − ∑ b ∈ B n under G b , t in . (12)

This study uses the available heating H ^av to describe the thermal flow in a pipe. Two strong-coupled variables, temperature T and mass flow q in the DHN model are decoupled to form the thermal flow network model and the basic flow-temperature equation.

The linearized heat loss equation shows as Eq. 14, and Wei et al., (2017) prove that it has good accuracy when the initial temperature is 88 ∼92°C. Section 4.2 will introduce the improved distributed algorithm based on linear heat loss. H av = μ q T − T rw (13) Δ H u v av = 2 p i T sw − T e ∑ R l u v . (14)

The setup of the reference direction of heating medium flow in the DHN is shown in Figure 4. The EH _n represented the nth EH. The thermal flow network model is shown as Eq. 15. H n , u , t EH + ∑ m ∈ Ψ u H u m , t av = 0 H u m , t av = − H m u , t av − Δ H m u av ,  if  H m u , t av > 0 H u m av , min ≤ H u m , t av ≤ H u m av , max ,  if  H u m , t av > 0 . (15)

A mixed-integer model of the thermal flow network is established based on the abovementioned equation to obtain the available heating distribution. Combined with the basic flow-temperature equation below, we can calculate the transmission temperature and thermal flow in the pipe. H u , t EH = k q u , t EH T u , t EH − T rw H u v , t av = k q u v , t T u v , t − T rw T u , t EH = T u , t ,  if  H u , t EH > 0 T u , t EH = T sw ,  if  H u , t EH < 0 T u v , t = T u , t ,  if  H u v , t av > 0 q u , t EH + ∑ m ∈ Ψ u q u m , t = 0 q u v , t + q v u , t = 0 . (16)

The DHN operation cost is caused by the operation cost of circulating water pumps in the network. Therefore, the EHR is introduced to describe the cost of the circulating water pumps. Eq. 17 represents the heat network operation cost of the nth ECU. Each pipe is configured with a circulating water pump in this study. C n heat = ∑ t ∑ p = 1 W EHR p c e , t H p , t . (17)

The EH is the energy conversion main body of the ECU. In this article, the standard matrix model is adopted for the EH model (Wang et al., 2017), and the piecewise linearization method is used to improve the accuracy of the EH model (Huang et al., 2019).

The energy flow direction of the coupling equipment is shown in Figure 5. Eqs 18–20 describe the standard model of the EH. The heating output from the EH and the heating obtained from the DHN jointly support the heat load in this region. V = P in  Grid  P in  EC P in  HP P in  CHP P in  GB P out  EC P out  HP P out, e  CHP P out, h  CHP P out  GB P in  HS P out  HS T . (18) V 1 = P EH G EH L cool L elec L heat − H EH 0 0 T (19) Z × V = V 1 . (20)

To simplify the nonlinearity of the coupling equipment efficiency, this study uses the piecewise linearization method as shown in Eqs 21–23 dealing with the nonlinear efficiency function of the coupling devices instead of the constant efficiency. X = X 0 + ∑ γ ∈ Γ σ γ . (21) F L X = F X 0 + ∑ γ ∈ Γ η γ σ γ . (22) I γ + 1 X γ ̄ − X γ ̲ ≤ σ γ ≤ I γ X γ ̄ − X γ ̲ . (23)

The model of the ECU has been completed. The internal power grid and the natural gas network build up the SOCP, while the DHN and EH build up the MILP problem. This study takes the minimum economic cost as the optimization objective, and the objective function and the constraints of each ECU are shown in Eq 24. min ∑ t c e , t P n , t ECU + c g , t G n , t ECU + ∑ p = 1 W EHR p c e , t H p , t s . t . 2 − 16 , 18 − 23 . (24)

According to the partition method, the insertion of virtual nodes into the energy networks and partitioning them is based on the principle stated in section 3.1. The ECUs can completely contain all the nodes without repeat, as shown in Figure 8A. The partition details are shown in Figure 8B.

This article uses the improved C-ADMM algorithm to complete the distributed optimization of the MRIES in two scenarios.

Scenario A: Optimizing the MRIES containing electrical, gas, and heating networks based on the ECUs.

Scenario B: Optimizing the MRIES containing the electrical and natural gas networks based on the ECUs.

In both the scenarios, the MRIES have the same loads. Nevertheless, in Scenario B, the coupling equipment produces the heating to meet the heating demands.

Conversion efficiency functions of the coupling devices are shown in Table 1. The types of coupling devices in the MRIES are shown in Table 2. The energy prices and loads of each EH are shown in Figure 9.

The codes were written in MATLAB 9.4.0.813654 (R2018a), and all the experiments were conducted on a desktop with 3.00 GHz Intel Core i5 and 16.0 GB RAM.

The improved algorithm proposed is used to optimize this MRIES in this article. The energy inputs to the equipment of the four ECUs in the different scenarios are shown in Figure 10.

As shown in Figure 10A, to support the heating load concentrated in midday, the HP is constantly working during low loading time of energy consumption when the power is cheaper. The heating is stored in the HS at night, which releases when the electricity price is high. In addition to this, the GB supports the remaining heating demand. Since there is no GB in EH3, it can only generate heat through the CHP. Furthermore, EH3 has the largest electric load at midday. As a result, the CHP runs at full capacity at all the time. In addition, the CHPs in EH2 and EH4 also run at full capacity at midday when the electric and heating load is large. On the other hand, the electric and heating load is relatively small in EH1. Therefore, the CHP is mainly used for power and heating supply in the peak load period, and the HP will start up for heat supplement during the flat power demand period.

As shown in Figure 10B, EH4 uses HP and GB with a high cost to supply heat due to the lack of DHN for heat exchange in scenario B. EH3 no longer needs to provide the heating to the DHN, which produced significant heating in scenario A. As a result, the utilization rate of HP and HS becomes lower. The overall energy cost response of the MRIES weakened in scenario B.

The heating generated by the EH, which is beyond the local heating load, is transmitted to the DHN. As shown in Figure 11, the energy represented by the blue block is equal to the energy of other color blocks higher than the heat load line (except the purple block that indicates heating storage). It should be noted that for the intuitive description, the load curves in Figure 11 are superimposed of electricity, gas, and heating load curves.

Because the coupling of the ECU in the electric network and natural gas network is only reflected in energy transmission, the energy interaction of each ECU in the DHN is mainly analyzed. The energy outputs of each EH device and the energy interaction with the DHN are integrated into Figures 11A–D. Evidently, the DHN realizes the heating interaction between the ECUs and dramatically increases the flexibility of energy utilization of the system. In scenario A, the EH1 internal equipment composition is relatively simple, which mainly absorbs the heat from the DHN. Although EH4 has many types of internal heat generation equipment, it still mainly absorbs the heat from the DHN due to its large heating load.

From Figure 11E–H, we can see that the heating support of EH3 is significantly reduced in scenario B. EH1, EH2, and EH4 have to convert energy through the coupling equipment in their units to meet the local loads, confirming the previous analysis conclusion. The total energy cost of scenario A is 12,281.24 dollars while that of scenario B is 15,234.56 dollars. For now, we can get the importance of the multi-energy cooperative scheduling.

The accuracy and convergence performance of the algorithm proposed are improved in this article. As analyzed previously, the improved algorithm can realize the MRIES global optimization through the unit internal optimization and information interaction between the units. The distributed optimization is not based on the whole system data; it can only approach the result of the centralized algorithm infinitely but cannot be better than that, and the insufficient iteration of the MISOCP will aggravate this error. The centralized algorithm which optimizes based on the whole system data can provide a precision reference for the improved distributed algorithm. The improved C-ADMM algorithm in this article adds the nested subalgorithm, which can not only overcome the nonconvergence problem caused by the MISOCP but can also approach the feasible region more closely than the C-ADMM algorithm, making the results more accurate. The results of the centralized algorithm, C-ADMM algorithm, and improved C-ADMM algorithm are shown as Table 3. The total energy cost error of the improved C-ADMM algorithm is about 0.1%. In contrast, the error of the C-ADMM algorithm is about 1.3%. In addition, the electricity and gas cost error are also relatively less. This shows that the improved algorithm proposed in this article achieves better optimization results than the original C-ADMM.

By comparing the convergence performance of the improved algorithm proposed in this article with that of the C-ADMM, the improved algorithm meets the accuracy requirement after 40 iterations. In comparison, the original algorithm needs 90 iterations under the same parameters as shown in Figure 12. In this article, the unit internal optimization and the information interaction between the units with the cloud are integrated and simulated on a computer. According to the simulation statistics, the C-ADMM algorithm takes 79.71 s, while the improved C-ADMM algorithm takes 40.76 s. Therefore, the convergence performance of the improved algorithm is better than that of the original algorithm.

The rationality of the distributed optimization with the ECU as a subsystem was verified by the improved algorithm. At the same time, the improved algorithm is proven to be superior to the original algorithm. The ECU provides a standard unit form for the improved consensus-ADMM algorithm, and their strengths are closely combined to achieve good optimization results.