Operation optimization of distributed multi-energy multi-microgrid considering system robustness

Carbon emission trading is regarded as an effective way to combine energy economy with green and low-carbon, which brings new vitality to the traditional multi-micro grid day-ahead dispatch. In this paper, a robust decentralized energy management framework is proposed for monitoring a collaborative structure of gas turbines, gas boilers, ground source heat pumps, energy storage and electrolyzers, etc for microgrid in the presence of power to gas and carbon capture systems. The price sensitivity of the power market results in the fluctuations of multi-micro grid dispatching. The worst-scenario uncertainty of multi-micro grid is managed by adopting trading prices at different conservativeness levels. The kalman filter distributed algorithm based on iteration is used to decompose the dispatch problem to minimize the total daily overhead of the multi-micro grid system while protecting microgrid data privacy. Finally, the simulation results represent the effectiveness of the proposed decentralized model of trading prices to meet the demand for electricity and heat. At the same time, the kalman filter distributed algorithm is compared with the alternating direction multiplier method algorithm to ensure accuracy and speed.


Introduction
Power generation mainly relies on fossil fuels, and the burning of fuels brings global warming and environmental damage Wu et al., 2021). The replacement of renewable energy (RE) i.e., wind and photovoltaics (PV) can effectively reduce carbon emissions. Nevertheless, the Table 1 compares the proposed model with the previously reviewed literature. The paper establishes an MMG robust decentralized energy management model based on electricity and heat under the condition of electricity price uncertainty. In the process of modeling, the DR of the load is analyzed, and the technology of converting P2Gis added to the CCS, so the efficiency of the system is significantly improved. The contributions of this paper are as follows: • In order to meet the DR requirements of thermoelectric load, a new MMG structure based on electric energy and thermal energy is established. The amount of carbon dioxide released is controlled by CCS system. • The optimal dispatch problem of MMG system is analyzed in detail. Considering the price uncertainty, a worst-scenario robust decentralized optimization method is proposed and its performance is tested. • A decentralized energy management method is proposed for efficient control of thermoelectric MG. KF was applied in modeling to disperse the model.

Frame structure
In this section, the MMG system structure is introduced in detail. The MMG system comprises three MGs, and each MG is based on electrical and heat, as described in Figure 1. Load requirements of MG, i.e., electricity and heat, are met through its own generation units and the energy market. GT, GB, P2G technology, RE, ES system (electrical and thermal storage), CCS, heat, and electrical load are applied to each MG, as shown in Figure 1. In addition, the P2G technology, consisting of electrolyzer and methanator, serves as a key component to fully convert the power in each MG into gas. CCS that provide CO 2 in the proposed model is coordinated operation with P2G. Noted that, carbon capture power plan, absorption tower and stripper make up the CCS system and the P2G technology consists of an electrolyzer and a methanator. Once the MG modeling is done, an appropriate management controller is required. Each MG is separately dispatched in decentralized model, and each MG is not informed of other MGs. Figure 2 exhibits the Bi-layer distributed framework of the MMG system. As shown in the figure, a coordinator is responsible for managing the exchange of information and energy power between MG and the energy market. Each MG dispatches optimal energy power in its control center after taking into account energy prices and local loads. First, MG operators contact MSO to submit their energy power. After trading power, MSO trade with the energy market to meet the demand of MGs. The distributed optimization in Bi-layer is mainly to protect the privacy of each MG.
Bi-layer problem formulation under the deterministic model

Mg model
Each MG as a unique agent contains a controller responsible for control and management for itself, optimally. For the specified MG, consider the objective function and relevant constraints when operation in grid-connected mode Feng et al. (2022); Mansour-Saatloo et al. (2021):  • Objective function of each MG: ( p tra HL (P tra, + ,t HL + P tra,-,t HL )) (6) where p gas is the price of nature gas. PRE is the price of abandoned RE (P cur,t RE ), i.e., PV and wind. p ES is the maintain price of P dis,t ES and P cha,t ES for ES. p cut EL , p tra EL , and p tra ]HL are the price of curtailed electrical load, i.e., P cut,t EL and deferrable thermoelectric load, i.e., P tra, + ,t EL , P tra,-,t EL , H tra, + ,t

HL
, and H tra,-,t HL , respectively. V t GT , V t GB , and V t P2G are the nature gas volume for GB, GT, and P2G. p i v is the price of carbon emissions volume E v . Noted that E 0 is the allocated rated volume for carbon emission. σ v is the increment. v is interval, v =0,1,2,3,…, V. t is the time, which is used for hourly dispatch, t =1,2,3,…, 23, 24.
Seven terms constitute the objective function of the χth MG. Term one refers to the cost of nature gas from the GT, GB and P2G. Second and third terms are the cost of abandon RE and ES and fourth term is the cost of curtailed electrical load. Terms five and six indicate the cost of sheddable thermoelectric load. The last term illustrates carbon transaction cost. Carbon price is an important factor affecting the carbon trading market. When the emission reduction task is small, the tiered carbon price can restrain the power generation in high-emission areas and reduce the total cost of the system Hu et al. (2020). The constraints relating to each MG mechanism are as follows: • Constraints for GT, GB, and GSHP: where P max ,t GT , H max ,t GT , H max ,t GB , and H max ,t GSHP denote the maximum power of P t GT , H t GT , H t GB , and H t GSHP for GT, GB and GSHP. LHV is lower heat of combustion for natural gas. η P GT , η H GT , and η GB are the gas efficiencies for GT and GB. η GSHP is the energy conversion efficiency for GSHP.
Equations (9)-(11) indicate the volume of natural gas required to produce thermoelectric power for GT and GB. Energy conversion constrain is used by equation (12). Equations (13)-(16) are the limitation of power output for GT, GB, and GSHP.
• ES constraints: where S t ES shows that ES allows the value of the state of charge (SOC) for MG χ , which is required not exceed the S MAX ES of ES. P dis,t ES and P cha,t ES are the charge/discharge power of ES, bounded by P max ES . μ cha,t ES and μ dis,t ES are the boolean variables, which indicate that charging and discharging power cannot be performed for MG χ at the same time. η cha ES and η dis ES are the efficiencies for the charging and discharging power.
The constraint expressing the SoC of the ES and the permissible limit is shown in equation (17). Equations (18) and (19) represent the allowable power limits for the charging and discharging states of the ES. Equation (20) holds cannot charge and discharge at the same time of ES.
• P2G and CCS coordinated operation constraints: where P max CCS $ is the maximum value of electric power consumed(P t CCS ) for CCS. Q t CC refers to CO 2 quality required for P2G, bounded by P max CCS . η P2G , and L P2G are P2G efficiency and calorific value of natural gas for P2G, respectively. α CC , α GT , and α MGχ are CO 2 quality consumed per unit for P2G, GT, and MSO. K CC , K GT , and K χ MG are CO 2 emission intensity of electric power purchased for CCS and GT, respectively.
• DR constraints: where μ tra,+,t EL , μ tra,−,t EL , μ tra,+,t HL , and μ tra,−,t HL are binary variables indicating the up and down adjustment volume of shifted thermoelectric load for EL and HL. P tra, max EL and P tra, max HL are the maximum value for thermoelectric shifted load. T cut, min EL and T cut, max EL are the maximum and minimum time of curtailed electrcial load. μ cut,t EL is binary variable indicating electrical cutting load.
• RE and transaction power with MSO constraints: where P P,t RE indicates the predicted RE output of P t RE . P max MGχ is the limitation of P t MGχ .
• Thermoelectric power balance constraints: Equation (45) is an electrical power balance constraint. Similarly, the heat power balance constraint is referred to in equation(46).

MSO model
As the medium between the power market and the MG, MSO is mainly responsible for the coordination in between them (Armendariz et al., 2017;Bendato et al., 2017). Therefore, in the case of energy trading, MSO, as the leader of MG, seems to be able to obtain more benefits by giving the transaction volume of MG participating in the market (Karimi et al., 2021). However, any time energy trading is uneconomic, it needs to be done in an optimal way from an economic and technical point of view. Therefore, the operation cost should be determined according to the objective function of the MSO.
• Objective function of MSO: where p t buy and p t sell show the price of buying and selling electrical power between MSO and electricity market, i.e., P t buy and P t sell .
The objective function of MSO includes the term referred to the cost of buying power from the electricity market. In addition, when the MG decides to sell power to an MSO, the cost of purchasing power from the MSO needs to be subtracted from the benefit in certain time periods, and it is clear that electric power cannot be bought and sold between the MG and the MSO at the same time.
The six-part formulation is similar to the MG χ constraint mentioned in the previous section.
• Transaction power with electricity market constraints: where P MSOχ indicates the electric power between MG χ and MSO. Equation (48) is the consistency constraint of transaction power for MG χ and MSO. Equation (49) shows the transaction power participation in the electricity market. The limits for P t sell and P t buy are refered to in Equation (50).

Robust optimization
In this paper, the electricity price is considered to be fluctuating (Lu et al., 2020;Mansour-Saatloo et al., 2021;Nikmehr, 2020). The electricity price is mainly effectively mitigated through RO, which is mainly related to its impact on decision-making, but leads to operating costs. Among these models, RO method is mainly used to analyze the fluctuation of electricity price. The RO method considers the case of low transaction price to obtain the optimal solution of this kind of problem.
In addition, the robustness of the dispatch problem is controlled by an integer parameter Γ, called the indeterminate budget, which determines the amount of risk the coordinator wants to take. In this procedure, the RO method is applied to a specific amount of Γ, and then the number of uncertain hours in the worst-scenario is determined by optimizing Mansour-Saatloo et al. (2021). Taking ϵas the deviation of the trading price, the objective function can be restated as: In the min-max model of the objective function (51), the inner item can be defined as a penalty term considering the power price is provided in worst-scenario and operation cost of outer item is minimized. The sum for index of t is limited to the Γ. The complex min-max problem is transformed into a linear min problem by using duality theory Shahidehpour (2002). The reformulate (51) is translated into a mixed-integer linear problem (52)-(56).
where α and β t are dual variables in equation (51) that makes the problem non-linear, y t is added as an auxiliary variable in MSO to avoid nonlinearity. Γ is uncertainty budget.

KF distributed algorithm
Algorithm description. The Kalman theory needs to synthesize the system state at the next moment according to the observations of two independent linear/non-linear systems. Now, The MMG containing MSO system is divided into two subsystems, i.e., the MSO and the observation system (MG). The coupling variable between the two is the energy transaction Power (57). Although MG itself cannot be linearized, its decision-making behavior and optimization process can be described by a linear system based on the KKT condition. Figure 3 shows the MMG architecture of the KF distributed algorithm. It is emphasized that the proposed KF distributed algorithm in this paper will use the KF theory in the "selection domain", i.e., predict the k + 1th value according to the value of the kth generation: (57) MMG based on MSO for KF distributed algorithm. In the iterative process, the initial solution (P t,0 MSO , P t,0 MG ) of the coupled variables may not satisfy the newly added constraints, which the models (1) and (52) are processed by Lagrangian relaxation.
where u t MSO and u t MGχ are the dual variables for the lagrange functions (58) and (59). θ t is an intermediate variable for distributed optimization.
The KF distributed algorithm continuously adjusts the Kalman gain K t to achieve convergence in the iterative process. Assuming that the white Gaussian noise in the linear system is ignored, the coupled variable is described as: Energy Internet operation framework Where μ t,k and ν t,k are the input signals of the linear system MSO and MG in Figure 3, respectively, where the initial value is 0. F t,k and H t,k are state transition matrices. Λ (•) denotes a diagonal matrix with "•" as elements. When the optimal solutions of models (58) and (59) are obtained depends on the KKT conditions, the penalty terms in the objective function are 0. Models(58) and (59) are derived: Compared with equations(60)-(61), the optimal solution is obtained with a state transition matrix: F t,k = H t,k = I. Meanwhile, θ t,k χ needs to be updated according to the KF theory, i.e., where k|k − 1 is the previous step (k-1)th. The Kalman gain K t,k can be updated by equations (66)- (70): where Z t,k refers to the covariance matrix of the updated state variableK t,k = Z t,k|k−1 (H t,k ) T (H t,k Z t,k|k−1 (H t,k ) T + R) −1 s. Q and R are the systematic error for the linear systems MSO and MG. KF distributed algorithm architecture is solved by distribution in Table 2. Algorithm: KF distributed algorithm Input: θ t,0 , Z t,0 , P t,0 MGχ , P t,0 MSO , Q and R 1: k ← 1 2: for t = 1 to 24 do 3: Solving Lagrangian functions (58) and (59) in parallel, to obtain P t,0 MGχ and P t,0 MSO 4: Update θ t,k based on P t,0 MGχ , P t,0 MSO and equation (60); 5: While θ t,k χ − θ t,k−1 χ ≥ ς do 6: k ← k + 1 7: Update matrix F t,k , H t,k , Z t,k and K t,k based on equations(66)-(70); 8: End while Output: θ t,k , P t,k MSO and P t,k

Case study
The decentralized robust day-ahead dispatch optimization was described in multi-energy system composed of three MGs and an MSO, as shown in Figure 1. Each MG is equipped with three units for GT, GB, and HB; P2G technology; ES as energy storage system; P2G-CCS operation system; wind and PV as RE sources. Figure 4 shows the electric trading prices for the next day. The parameters of generation units and ES are referred to the literatures (Gao et al., 2017;Rueda et al., 2013). The model in this paper is tested on a PC with a main frequency of 2.8GHz and a memory of 32GB using the Matlab R2020a simulation platform, and the YALMIP optimization toolbox is embedded to call CPLEX for optimization.
The day-ahead dispatch deterministic result of MMG Figure 5 indicates the optimal dispatch of thermoelectric power for each MGs separately. Figure 5 use a deterministic method to solve optimization problem, ignoring the uncertainty of the energy market. Figure 5 (a1)-(a3) are the electric power, (b1)-(b3) are the heat power, and 1-3 are the labels of the MG. MGs 1 and 3 are equipped with wind and MG 2 is equipped with PV.
Due to rising electrical load, GT of MG1 is called after hour t = 6, and operators prefer to use domestic technology to supply the loads. At the same time, the buying of electrical power is also rised to the operator. The electric load of MG2 is at its peak during t = 17-21, and the operator purchases electric power from the power market. This is due to the operation cost of GT is higher than power market prices. In MG3, the output of wind in a 24-h period is relatively large, and the electricity purchased from the electricity market is 2-3 times smaller than that of MG1. Every MG operator benefits from ES because of its ability to save electrical power during low-price periods and its ability to deliver power in peak hours. According to Figure 5(a1), ES charges during t = 2-4, 20, 23-24, Where the trading price of electricity is low and discharge when the price rises, i.e., t = 6-7, 11, 13. Due to the low power purchase price, CCS mainly obtains CH4 from the electricity market and GT, and P2G  Figure 6 shows the result of DR implementation for thermoelectric loads, respectively. (a1)-(a3) are the electric load, (b1)-(b3) are the heat load, and 1-3 are the labels of the MG. The operators encourage demand-side departments to play a role in DR dispatch, reducing operational costs and promoting RE utilization through peak cutting and valley filling. According to Figures 6(a1) to (a3), electricity DR to successful supply, performance during this period is positive. Take MG1 as an example, the MG shift the loads from on-peak power market hours, i.e., t = 6-19, to off-peak hours, i.e., t = 1-5, 20-24. As a result, MG can deliver a load of 5.6%at a lower cost. In addition, according to the heat DR results in Figure 6(b1), for MGl during hours t = 9-18, heat demand was shifted up, which is need to shift up the heat load at the moment.

Tiered carbon price analysis
This part compares the operating costs of carbon Price using two cases: Case1: The tiered carbon price method mentioned considering P2G-CCS system in this article. Case2: The uniform carbon price adopted in each MG considering P2G-CCS systerm. Case3: Regardless of P2G-CCS technology in the tiered carbon price method i.e., keep CCS and P2G separately. Figure 7 shows the carbon balance for each MG in Case 1. The CO 2 required by each MGs comes from the GT, which is supplied to P2G and CCS to convert it into natural gas. Take MG as an example, the carbon required for CCS is at its peak during t = 1-7,18-24. Just in time, CO 2 produced by GT is also in a peak state. At t = 10-11, CCS requires less carbon.
As can seen from Table 3 that the CO 2 discharged into the air in case3 is (17914, 20948, 6828)(m 3 ) less than (68, 81, 93.1)(m 3 ) in case1(17846, 20867, 6734.9)(m 3 ).However, the proportions of CO2 emissions into the air are 50.15%, 54.83%, and 46.0% in Case 1, and the proportions are 52.0%, 56.68%, and 50.34% in Case 3. It can be seen that adding CCS to the system can effectively reduce carbon emissions. Figure 8 gives two cases about carbon price. Compared with Case1, the carbon emission of Case 2 is slightly higher by (0.37%, 0.51%, and 0.49%) for each MGs, which utilize the total carbon in 24-h. When the emission reduction task is small, compared with the uniform carbon price, the tiered carbon price can restrain the GT to reduce the total cost and carbon emissions.

Robustness analysis of trading price
The uncertainty of energy market affects MMG trading power, as shown in Figure 9, with different deviations ɛ. During the worst-case scenarios, the larger ɛ, the larger the transaction volume of MG, i.e., t = 18-22 in MG2, t = 19 in MG3.
The impact of wholesale power market fluctuation on the operating cost of MGs has been analyzed, as shown in Table 4. Since operators are pessimistic about the predicted price one day into the future, handling uncertainty in a conservative manner increases operation costs in the case of RO. Consider worst-case operation costs for 24 h more than 8 h (1.5%, 1.8%, 0.6%, and 6.8%). As can be seen from the table, the more the uncertainty budget is increased, the more worst-case hours to consider the dispatch, the more operating costs are obtained.

Convergence analysis of KF distributed algorithm
In the last part of the results, the performance of the KFDA technique is analyzed. Figure 10 presents the total running cost of each iteration of MMG at Γ = 24, while comparing with the ADMM algorithm. As described below, KF distributed technology converges after 6 iterations, while ADMM requires 7, demonstrating its ability in terms of computational speed and scalability. In addition, using the centralized problem for the MMG model, its total operation cost is 88943$, which is 2.5% more than the KF distributed algorithm, which confirms the accuracy of the KF distributed algorithm.

Conclusions
This paper established a decentralized energy management platform based on robustness, and makes a comparative study of its performance. At the same time, it also studies the optimization of thermoelectricity-based MG dispatching when the P2G-CCS system exists. The operator of MGs can effectively control the dispatching result of the next day by using the RO method to deal with the price fluctuation problem. In addition, MG operators can reasonably distribute carbon emissions, which plays an important role in promoting community economic development. The comparative analysis results show that the MG saves 10.12%in cost. At the same time, P2G-CCS technology can also promote the MG, and the reduction of carbon emissions can reach 0.39%. Based on the RO method, the risk avoidance of the MG is carried out. The simulation results show that the KF distributed algorithm converges in a short time after it is introduced. It can be judged that it is of great significance to improve the accuracy of the model, and it is also an important object for the future development of the distributed system.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Changzhou Science and Technology Project: Study of the control mechanism and application of constant radiation temperature prediction model under radiation cooling, the National Natural Science Foundation of China (NSFC) (grant number CJ20220178, 61873336).