Cooperative Operation for Integrated Multi-Energy System Considering Transmission Losses

The increasing complexity of energy internet integrated with multiple energy systems requires a more flexible energy management system. To exploit the potential of energy complementarity, maximize the consumption of renewable energy sources (RES) and minimize the operating costs, a cooperative operation model for integrated multi-energy systems (IMES) with transmission losses is proposed in this paper. On this basis, a novel multi-energy management strategy is developed to improve the energy utilization efficiency of the systems for regional operation by integrating the alternative decomposition-based decoupling method for multi-energy flow and the dynamic wolf pack algorithm (DWPA). Furthermore, the operating costs are formulated as a multi-objective optimization problem based on the Pareto efficiency to obtain more diverse solutions. Finally, case studies considered perform sensitivity analysis and demonstrate that the proposed cooperative operation model and strategy can effectively improve the economy and the flexibility of a regional energy system.

With the growing demand of diversified energy resources and the excessive consumption of fossil fuels, the integrated energy system (IES), which is capable to coordinate diversified energy resources to generate desired energy and benefit the economy, efficiency and flexibility of the multi-energy system, has been recognized as an effective alternative to traditional energy utilization methods in the context of energy internet [1]- [4]. The complex energy couplings inside an IES are simulated by using the energy hub (EH) model [5], [6]. Compared with the individual IES, the integrated multienergy system (IMES), which couples energy networks such as utility grid (UG), gas network (GN), heat network (HN), and several neighboring IESs, can be scheduled simultaneously to supply energy demands. Thus, the IMES can further improve the reliability and flexibility of the regional energy supply.
One of the most concerning issues is to develop a dispatch model for the IES. The objective of such a dispatch model is to minimize the IES operating costs while satisfying control constraints, such as encompassing energy balance, device capacity, ramping rate limitation, etc. [7]- [9]. Advanced modeling techniques are needed to fully simulate the energy flow inside an IES [10], [11]. The two-level economic dispatch model was proposed in Ref. [12], [13]. Ref. [12] presented the decision-making conflict among different market players. The market players must find a long-term equilibrium through price signals in a multi-energy retail market. Ref. [13] proposed a security-constrained model for integrated gas and electricity networks. The lower-level model is formulated based on Karush-Kuhn-Tucher (KKT) conditions, thus the problem can be solved by using mix-integer linear programming (MILP). Ref. [14] proposed a scenario-based stochastic IES model, an uncertainty matrix was employed to tackle the uncertainties of the RES and loads. A linear coupling matrix was proposed in Ref. [16], which facilitated the computerized calculation for an EH with arbitrary configurations. Ref. [17] proposed an optimal multi-energy flow for an IES in the carbon trading market. Ref. [20] presented a model to provide multi-energy demands and to manage continuous on/off controllable equipment. To facilitate interactions between power companies, residential buildings, and industrial consumers, a system for interconnected infrastructures based distributed multi-energy resources was established [18], [21]. A model considering energy trading between the IES and energy networks, as well as the limits of battery systems, gas stations, and thermal storage systems, was proposed in Ref. [19]. Ref. [22] provided a flexibility assessment for IES operations. The flexibility indicator was assigned as the suppression of the uncertainty of the RES and the minimization of the operating costs. Ref. [23] investigated the effects of ice storage on the economic performance of the IES while satisfying demand response in an uncertain environment. Ref. [24] provided a management framework which consists of IES and microgrid. Moreover, the impact of the wind turbine (WT), photovoltaic (PV), and tidal generation on the planning of the IES was investigated. Ref. [25] presented an IES model that minimizes the operating costs subject to the power flow constraints and limits of technical indexes in electric vehicles (EVs) parking lots, and CHP constraints. A two-level problem for optimal day-ahead scheduling of a distribution system, which trades electricity with multiple IESs was presented in Ref. [38].
In addition to the research works that focus on the dispatch model, optimization algorithms for the operation of the IES have been proposed in some studies, such as the Monte Carlo method [26], the probability statistical method [27], the grasshopper optimization algorithm (GOA) [30], the evolutionary PSO algorithm [33] and the stochastic optimization method [28]. The ε-constraint method was employed to solve the IES coordinated operation in Ref. [15]. Ref. [29] proposed a MILP based on a multidimensional piecewise linear approximation method to calculate the optimal power flow for the IES. The proposed method can effectively solve a non-convex control problem for the gas networks. Ref. [31] introduced the coordinated scheduling strategy and a multi-objective optimization algorithm was used to obtain a Pareto optimal solution. Then the interval evidence reasoning was used to analyze the multi-attribute decision-making to determine an optimal operation solution of the IES. Ref. [32] used the distributed optimization method to find out a bargaining solution for the cooperative system ensuring the autonomous scheduling and information privacy. A decentralized algorithm was proposed in Ref. [34] to schedule the cogeneration for an IES by using the Benders decomposition. The calculation methods for the optimal IES power flow are discussed in [35]- [37]. Ref. [39] proposed a multi-objective function solving algorithm based on the tent map chaos optimization, which can reduce the operation costs and risks associated with an IES. The schedule plan of each energy device was optimized by the moth flame optimization algorithm in Ref. [40], which can effectively reduce operating costs, and the optimal output plan of each device can be obtained.
It is worth noting that most of the research works in the aforementioned literature focus on the optimal operation of a single small-scale IES. Moreover, most of the energy resources in the IES have limited capacity allocation and self-coordination capabilities. Therefore, it is difficult for a single small-scale IES to have better economic benefits. On the other hand, it is also very hard for large-scale energy storage in the IES to be fully utilized due to the high operation cost and geographical restrictions. Thus, the UG is often used as a compensatory source to solve the power exchange problem. However, such a solution will lead to overburden and fluctuating operation for the UG. The technologies focusing on energy complementarity in the IMES can provide an innovative way to accommodate energy storage. This local community as a local energy resource can ease the cooperative operation of IESs.
To exploit the potential of the multi-energy complementarity among various IESs and achieve higher efficiency, this paper develops a cooperative operation IMES based on the EH concept and proposes a cooperation based control strategy. The main contributions are summarized as follows: 1) Based on the linearized coupling matrix of the EH model, a cooperative operation model is presented for the IMES to combine the residential EH (R-EH), the commercial EH (C-EH), the industrial EH (I-EH), the UG, the GN and the HN. Furthermore, the transmission losses is developed, which can further study the actual supply-demand balance and avoid serious economic losses.
2) The operating costs are formulated as a multi-objective optimization problem based on the Pareto efficiency to obtain more diverse solutions.
3) A novel multi-energy management strategy that integrates the alternative decomposition-based decoupling method of multi-energy flow and the dynamic wolf pack algorithm (DWPA) has been further developed.
4) The DWPA is presented, where dynamic scout direction, Fibonacci sequence, avoiding inferiority, and Levy-flight function are introduced as an improved function. Compared with other algorithms, the DWPA has advantages in terms of operating costs and time.
5) The case studies demonstrate that, compared to the non-cooperative mode, the proposed control strategy can achieve better cost-effective performance. Whilst, the sensitivity analysis is performed to explore the impact of the confidence of the objective function and constraints on the operation of IMES.
The rest of the paper is arranged as follows. Section II presents the coupling model for EH, the energy network models, and the cooperative operation model for the IMES. The proposed multi-energy management strategy is introduced in Section III. Section IV demonstrates the case studies and the effectiveness of the proposed strategy, while Section V concludes the paper.

II. MODELING APPROACH FOR THE IMES A. LINEARIZED COUPLING MODEL OF EH
The flexibility of the EH mainly comes from the devices and the bifurcations of energy flows. Fig. 1 illustrates the topology of a typical EH, which consists of a combined heat and power (CHP), an electric boiler (EB), a gas boiler (GB), a compression chiller (CC), an absorption chiller (AC), an electrical storage (ES) and a thermal storage (TS).
The linearized coupling relationship for the EH can be expressed according to the energy efficiency matrix (1), as shown at the bottom of the next page. The energy ES , η ch TS , and η dis TS are the conversion efficiencies of energy devices, respectively. The cooperative operation model for the IMES is shown in Fig. 2, which couples five EHs, a UG, a GN, and a HN. As shown in Fig. 2, the IMES contains multiple EHs and the R-EH1 contains the PV, and the CC; the R-EH2 consists of the WT, the HP, and the EB; the R-EH3 contains the CHP, the AC, and the ES; the I-EH consists of the CHP, the EB, and the AC; and the C-EH contains the GB, the CC, and the TS. The surplus energy resources can be locally shared between the EHs to supply energy demands. This cooperative operation model can be extended to a very big network with numerous EHs.

B. UG MODEL
In the IMES, the UG can be modeled using an optimal power flow model [41]. The integration of RES and energy coupling makes the power flow more complicated.
96936 VOLUME 8, 2020 The relevant calculations can be expressed as where P n (Q n ) is the active (reactive) power injection of the node n; V n is the voltage amplitude of the node n; V m is the voltage amplitude of the node m;θ n,m is the phase difference between the node n and m; G n,m is the conductance;B n,m is the susceptance;N is the node number. The node energy balance and power flow constraints can be described as where P n,m (Q n,m ) is the active (reactive) power of the power line between the node n and m; P max n,m (Q max n,m ) is the transmission power constraints of the power line; P max n , P min n , Q max n , Q min n , V max n , V min n are the upper and lower limits of the active power, reactive power and voltage amplitude, respectively.
The voltage stability constraints can be described as where V ref is the reference voltage value; THR is the threshold, which limits the allowable voltage range

C. GN MODEL
The model for gas flow balance can be formulated as gas flow into node k from pipe l −1 gas flow out of node k from pipe l 0 no flow between node k and pipe l (8) where A G is the node-branch incidence matrix; F is the gas flow in the pipe; G is the gas injection of the node. The relationship between gas flow and node pressure can be described as [42].
where f k,ρ is the gas flow in the pipe between node k and node ρ; κ k,ρ , ϑ are the pressure drop and friction factor of the pipe; d, L are the diameter and length of the pipe; ν is the specific gravity of the gas. The node flow balance, node flow and pressure constraints can be expressed as where f ρ is the flow of the node ρ;p ρ is the pressure of the the node ρ;f max ρ , f min ρ , p max ρ , p min ρ are the upper and lower limits of the flow and pressure, respectively.

D. HN MODEL
Generally, the HN model can be divided into two parts: the hydraulic circuit and the thermal circuit. The hydraulic circuit can be modeled as [43] A H = M k (12) water flow into node k from pipe l −1 water flow out of node k from pipe l 0 no flow between node k and pipe l (14) where A H is the node-branch incidence matrix; is the water flow of the pipe; M k is the water flow of the node; i is the loop label; j is the pipe label; B is the loop incidence matrix;R is the resistance coefficient. The thermal circuit can be formulated as where Q k is the heat energy of node k; C H is the specific heat capacity; L is the length of the pipe; m k is the water flow of node k; T is the ambient temperature; T I k , T O k are the inlet and outlet temperature of the node k;T in , T out are the inlet and outlet temperature of the pipe; T I ς , m I ς are the temperature and water flow into the node from the pipe ς;T O ς , m O ς are the temperature and water flow out of the node from the pipe ς ; is the thermal conductivity.

E. TRANSMISSION LOSS
To address the energy complementarity between different EHs, the distance between EH i and j is donated as L i,j . The power transmission losses E loss and the heat transmission losses H loss of L i,j that can be expressed by [47], [48] where P i is the electricity output of the EH i; B i is the loss coefficient; ρ r is the resistivity;U 2 r is the rated voltage;S r is the cross-sectional area of the resistance; R h is the thermal resistance of the pipe per kilometer from the heat medium to the surrounding medium; T w,k is the supply-water temperature at node k; T avg,ς is the mean temperature of the medium around the pipe ς. Assuming that the characteristics of the power lines among EHs are the same, thus the energy loss is proportional to the length L i,j .

III. CO-SCHEDULING METHOD A. OBJECTIVE FUNCTION
To fulfill control demands for different energy resources simultaneously and achieve the minimum operating costs f day for the IMES, the objective function is formulated as where f 1 is the cost of the energy procurement; λ t E,net , λ G,net , λ H ,net are the prices for purchasing electricity, gas, and heat, respectively; E t net , G t net and H t net are the energy procurements from the UG, the GN, and the HN; f 2 is the environmental costs; α f 1 , β f 2 are the weight coefficients; χ G is the penalty factors.
The probabilistic constraints can be described as where ∂ is the confidence of the objective function;f day is the objective function when the confidence is ∂.

B. OPERATIONAL CONSTRAINTS
The energy balance for all resources in the IMES can be described as where E t i,j , H t i,j are the exchanged electricity and heat energy between EH i and j; (21) describes the balance of the electricity, (22) formulates the balance of heat energy, and (23) calculates the balance of cooling energy.
The probabilistic constraints can be expressed as where φ is the confidence of the constraints; ζ 1 , ζ 2 are the maximum exchanged energy The CHP needs to meet the operational constraints and the ramping rate limitation, which can be described as where Q mid CHP is the heat energy output when the power generation of the CHP is E min CHP ; E max CHP is the maximum power generation; E mid CHP is the minimum power generation of the CHP; λ m is the electro-thermal elastic coefficient; E up CHP , E down CHP (Q up CHP , Q down CHP ) are the upper and lower limits of increase and decrease of electricity (heat) per unit time; t is the ramp response time.
The energy devices must operate within its rated capacity. Their capacity limits are formulated as The energy storage devices are constrained by the state of charge (SOC), which can be obtained as The proposed IMES model requires compliance with various non-linear operational constraints. However, the increasing system complexity brings considerable computational requirements to traditional approaches, making them technically less appealing. Besides, with the increasing scale of future IMES, the traditional approaches will face a series of problems, such as lower converging speed or deviation from the global optima. On the other hand, novel intelligent meta-heuristic algorithms have been employed to solve different economic dispatch problems and achieved satisfactory results [30], [33], [44]. In this section, a power flow optimization method based on DWPA is proposed to optimize the operation and dispatch the energy resources within the IMES. The wolf pack algorithm (WPA) selects the best wolf as the lead wolf, and then, the WPA divides other wolves into the scout wolves and ferocious wolves [44]. The prey behaviors include scout, summon, siege, and renew. The lead wolf summons the wolves to capture the prey according to the prey odor smelled from the scout wolves.
We first define x 1×H as the position of any wolf and it can be used as inputs for the objective function. The position x 1×H can be formulated as where the associated costs f day calculated using (19) is defined as the fitness value for x 1×H . The scout wolves explore the environment to hunt for the prey. The x α i,h can be updated by [44] x where α = (1, 2, · · · , s) denotes the prey direction;τ a h is the scout step; s is the scout direction; x i,h represents the position of the i-th scout wolf in the h-th dimension h = (1, 2 · · · , H ), which corresponds to the electricity, gas, and heat energy procurements.
However, the basic WPA will use a fixed s that leads to insufficient diversity and randomness. In order to improve the global search efficiency, a dynamic scout direction is proposed in this paper and the new scout behavior can be updated as where I h is an identity matrix of scout direction in the h-th trace; s 2h is the dynamic scout direction. If the fitness value of the i-th scout wolf is greater than the lead wolf, the scout wolf will replace the lead wolf and summon other wolves. The position of the ferocious wolf can be updated as where w ∈ [0, 1] is the odor weight; τ b h indicates the raid step; x j+1 i,h is the position of the i-th ferocious wolf at the j + 1 iteration, g j h is the position of the lead wolf. The step size of the ferocious wolf will be slowly reduced to approach the lead wolf. The siege radius is defined as where ε is the i-th distance determinant coefficient; ψ is the distance between the ferocious wolf and the lead wolf when initiating the siege behavior; ξ max h and ξ min h are upper and lower values.
The siege behavior is formulated as [44] x j+1 where τ c h is the siege step and γ ∈ [−1, 1]. If the fitness value of the i-th ferocious wolf is better than the lead wolf, the lead wolf will be replaced. VOLUME 8, 2020 The basic WPA uses a fixed step to update the scout wolf position and if the scout step is too large, the convergence precision will be attenuated. However, if the scout step is too small, it will slow down the convergence speed. To balance the convergence precision and speed, a dynamic scout step based on the Fibonacci sequence where χ > 1 is the learning rate.
In the context of the basic WPA, only the lead wolf has affects on the siege behavior of the whole wolf pack. However, in the DWPA, ferocious wolves can interact with each other and the siege behavior of the wolf pack will be updated using interactions from all ferocious wolves. The intention for that is to quickly escape from the area where it has a lower probability for the optimal prey. Therefore, the proposed avoiding inferiority can improve the convergence speed. The pseudo-code is given as follows.
Pseudo-code for Avoiding Inferiority 1. for each ferocious wolf i 2. if the distance between the ferocious wolf i and the lead wolf is smaller than ψ 3.
Run with step size τ c h using (48).

4.
if the distance between the ferocious wolf i and another wolf with the poorer fitness value is smaller than ψ in (47).
Run with step size τ c h using (48). 8. end 9. end 10. end for To avoid premature convergence, σ wolves which have relatively lower fitness will be renewed to maintain the diversity of the wolf pack. σ is a random value between [n/2δ, n/δ], where δ is the population renewing proportional coefficient. The Levy-flight [45] is also taken to enhance the diversity, which can be attained as where β ∈ [0, 2] is a fixed number and used to adjust the stability;µ, v are normal distributions; r is a random number between [0,1]; is the Gamma function.

D. DECOUPLING METHOD
The implementation flowchart for the decoupling method for the IMES is shown in Fig. 3. The traditional energy flow calculation method, such as the Newton Raphson method [43], needs to consider different system components (e.g., the UG, the GN, the HN, and the EH) as a whole, which complicates the calculation and increases the dimension of the system equations. The proposed IES model only considers the energy efficiency matrix calculation of the EH. Therefore, an alternative decomposition-based decoupling method can be applied for the optimal interaction of EHs and energy networks [46]. The decoupling method can calculate the energy flow for the UG, the GN, and the HN, respectively. Thus, a more accurate dispatch for the multi-energy flow can be obtained.

A. CASE SETTINGS
To demonstrate the effectiveness of the day-ahead IMES model and the proposed method, a test system is developed which consists of a UG (e.g., an IEEE 33-bus system), a GN (e.g., a 9-node low-pressure gas system), the HN (e.g., a 6-node heat system) and 5 EHs, as depicted in Fig 4. The EHs cover an urban area which consists of 1800 residences, 10 commercial units, and two industrial units. The Node E5, E13, E20, E23, and E30 are the coupling points between the UG and the EHs. The Node G2, G7, and G9 from the GN are connected to the R-EH3, the I-EH, and the C-EH, respectively. Except for the Node H1 of the HN, each network node of the HN is connected to one EH. It is assumed that the maximum heat loss in water pipe is 10%, the allowable voltage range is [0.95p.u. 1.05p.u.] and the allowable pressure range is [0.9p.u. 1.1p.u.]. Fig. 5(a), (b), (c) show the data for the energy demands of 5 EHs. The hourly outputs of the RES are obtained from day-ahead forecasts, as shown in Fig. 6(a). The curves of time-of-use price for electricity, natural gas, and hot water are shown in Fig. 6(b). The parameters of energy devices are listed in Table 1-2.

B. SENSITIVITY ANALYSIS AND PARETO FRONT
In order to compare the operating results, the total operating costs of the IMES under different confidence values of the objective function and constraint are selected.
As shown in Fig. 7, with the increase in the confidence value of the objective function and constraint, the operating costs of the IMES increase. It is demonstrated that as the confidence value increases, the uncertainty input is suppressed, the reliability requirements of the objective functions and constraints of the IMES will increase, and the reliability will be compensated with a higher operating cost.
To further explain the Pareto efficiency in IMES models, the Pareto optimal sets of costs are compared with other methods. The weights α f 1 , β f 2 are set to α f 1 ∈ [0 : 0.02 : 1], and β f 2 ∈ [1 : 0.02 : 0].   As shown in Fig. 8, the Pareto optimal sets can effectively reflect the relationship between the energy purchase costs and environmental costs. The point A indicates that the energy purchase costs are dominant, that is, the weight of the energy purchase costs is higher than that of the environmental costs. The point B indicates that the environmental costs are dominant. Other points can be served as the compromise solutions for decision-makers under different environmental and economic requirements. When the energy purchase costs decrease, the environmental costs increase.

C. OPERATIONAL SCHEDULING FOR THE IMES
Simulation results are summarized in Fig. 9 to show how the energy distribution, and conversion within the IMES. The objective function confidence and constraints confidence are both set to 0.9. The weight coefficients are both set to 0.5. VOLUME 8, 2020  With the expansion of energy couplings and complementarities, it is necessary to analyze the operational scheduling of different energy devices.
The electrical loads of the R-EH1 are mainly supplied by the PV and UG. The excess PV power can send to other EHs. Because the R-EH1 doesn't own any thermal supply, it needs to obtain heat energy from other EHs, as shown in Fig. 9 (a). Since the electricity produced by the WT is sufficient in R-EH2, it will be converted to heat and cooling energy, as depicted in Fig. 9 (b). It is cost-effective to purchase electricity from UG to supply electrical loads due to the lower price during the evening. In the cooperative mode, the ES is heavily used because it is efficient to charge or discharge to accommodate the fluctuations of the RES, as illustrated in Fig. 9 (c). For the I-EH, the electrical load is mainly supplied by the CHP, as shown in Fig. 9 (d). The extra heat generated by the CHP will be sent to other EHs. The EB can be used as an auxiliary heat source. The cooperative operation improves the operating efficiency of the CHP, which leads to an increase in gas consumption. Moreover, it is more cost-effective to choose the AC over the CC for supplying cooling loads because the AC can use the extra heat to supply cooling loads. The electrical loads in the C-EH can be supplied by the surplus electricity from other EHs, which can reduce the power procurement from the UG, as depicted in Fig. 9 (e). The TS is fully utilized to store the extra heat generated by the GB, which improves the system flexibility to handle the heat supply.

D. COMPARISON WITH OTHER METHODS
The computational results are given in Fig. 10 and Table 3. The result shows that in terms of solution quality, the DWPA can provide a sub-optimal solution that is very close to the optimal value obtained by CPLEX [13]. However, the CPLEX requests a longer time to solve the economic dispatch problem. Compared with the genetic algorithm (GA) [39], the GOA [30], and the WPA [44], the proposed DWPA can reduce the operating costs by 1.96%, 0.74% and 0.55%, respectively, which illustrates that the proposed method has better global search capability and convergence performance than other methods within a shorter computational time.

E. ENERGY COMPLEMENTARITY IN THE IMES
The simulation results of the energy procurements (e.g., electricity, gas, heat) of the IMES with or without cooperation are shown in Fig. 11 and the detailed operating costs are compared in Table 4. The results show that, compared with the non-cooperation mode, the IMES with cooperation needs lower electricity and heat procurements. The gas procurement is increased due to the increased operating efficiency of the CHP. When operating under cooperation mode, the IMES has more flexibility, which means the proposed model is capable of utilizing different energy devices to flexibly deal with the diversified energy resources to fulfill desired energy demands. The internal energy exchange among different EHs can lower the electricity purchase from the UG, which reduces the burden of the UG power. It can be concluded that the EHs of the IMES can rely on the local energy provided by  neighboring EHs to serve some of their energy demands and realize a higher utilization rate of RES.
The hourly operating costs is shown in Fig. 12. The cooperative operation outperforms the non-cooperation in terms of overall operating costs, as shown in Table 5. In addition, the operating costs of the R-EH1, the R-EH2, and the C-EH are lower because of the cooperation and internal energy exchange. However, the operating costs of the R-EH3 and the I-EH are increased. That this because the CHPs of the R-EH3 and the I-EH increase their energy output to compensate for the energy shortage of other EHs. Instead of solving the scheduling problem for each downward individual EH, the total operating costs of the all involved EHs are taken into account in the upper ward IMES. The costs are reduced by 4.3% leading to a more efficient energy utilization.

F. TRANSMISSION LOSS
Assume that the distance between the EHs of the IMES is depicted in Fig. 13. The extra heat and surplus electricity can be exchanged in the IMES.
As depicted in Fig. 14, the positive values mean that electricity and heat are inputted into an EH, and vice vasa.  Taking C-EH as example, compared to the results without considering transmission loss (the dark blue and green bar in Fig. 14), the electricity input (the blue bar in Fig. 14) is slightly decreased and heat output (the yellow bar in Fig. 14) is increased due to the transmission loss (TL). The result illustrates that an EH with energy outputs needs to increase its energy supplies rather than simply meet the energy demands while the extra energy is used to compensate for the transmission loss. Because the imported energy cannot satisfy the internal energy demands due to the transmission loss, EHs who need external energy supplies will also request more internal energy supplies. The operating costs of each EH are summarized in Table 6. All EHs have higher operating costs when considering the transmission loss for electricity and heat.

V. CONCLUSION
To alleviate the pressure of energy depletion and enhance the multi-energy complementarity, this paper presented a cooperative operation model for the IMES with transmission losses. The energy distribution and conversion within the IMES were modeled using linear energy coupling matrices and energy flow equations. The operating costs are formulated as a multi-objective optimization problem based on the Pareto efficiency to obtain more diverse solutions. On this basis, a novel multi-energy management strategy that integrates the alternative decomposition-based decoupling method for multi-energy flow and the DWPA was developed to optimizing the joint operating costs for the IMES. Compared with the non-cooperative operation mode, the case study demonstrates that the energy complementarity can improve the flexibility for energy systems and accommodate more RESs to reduce the power purchased from UG. Besides, with increasing confidence, the operating costs of the IMES increase. The proposed cooperative operation model and strategy can achieve more rational energy scheduling and better cost-effective performance, which can be further used in the optimization operation of energy internet.