Two-stage robust optimization of regional integrated energy systems considering uncertainties of distributed energy stations

For a regional integrated energy system (RIES) composed of an energy supply network and distributed energy station, the uncertainty of distributed photovoltaic (PV) output and the fluctuation of various loads pose significant challenges to the stability of system operation and the accuracy of optimal scheduling. In order to enhance the operational reliability of regional integrated energy systems and reduce the impact of photovoltaic and load uncertainties on distributed energy stations, this study proposes robust optimization method of regional integrated energy systems that takes into account the uncertainty of the distributed energy station. First, the regional integrated energy system is divided into an upper electric-gas energy supply network and a lower distributed energy station. The upper model mainly realizes energy transmission, while the lower model is a two-stage robust optimization model of distributed energy stations in the form of min–max–min, which mainly realizes flexible energy supply of different types of energy. Then, the lower two-stage robust optimization model is simplified and solved using a column and constraint generation (CCG) algorithm. After that, an alternating direction method of multipliers (ADMM) is used to solve the upper and lower models of the regional integrated energy system, and the solution scale is reduced while ensuring the correlation between the energy transmission network and the distributed energy stations. Finally, a test example is provided to illustrate the effectiveness and usefulness of the proposed method. It follows from simulation results that the robust optimization method can effectively reduce the instability of the system operation caused by uncertainty factors and improve the system’s anti-interference ability, and in addition, systems with high penetration levels of photovoltaic output will benefit more from robust optimization.


Introduction
In order to comply with the world's green and low-carbon development trend, it has become particularly significant to promote the integration of sources, networks, loads, and storage and to build a multi-energy-coupled regional integrated energy system (RIES). In addition to the exploration and research of new energy sources, the improvement of the energy supply form is also of equal significance (Zhao et al., 2022a). The role of the regional integrated energy system, as a promising solution of bearing energy for human society in the OPEN ACCESS EDITED BY future, is rapidly developing in the field of efficient energy utilization, energy conservation, renewable energy consumption, and emission reduction (Balsamo et al., 2020). Followed by the foundation of the RIES concept, countries worldwide have responded positively and vigorously by carrying out extensive research on the RIES (Zhang et al., 2022). Some key issues of the RIES have already been studied from various perspectives in China Zhou and Zhang, 2020;Zhao et al., 2022b;Rezaei et al., 2022;Yang et al., 2022), and many theoretical and practical results have been obtained.
A RIES can fully exploit the complementary characteristics of multiple energy sources to achieve multi-layer energy efficiency, which could also reduce operating costs and carbon emissions effectively, and further consume renewable energy that cannot be consumed by traditional single energy systems through the conversion between different energy sources (Sioshansi et al., 2021). However, a RIES is much more complex than the traditional energy system and has various coupling forms, which need higher requirements on optimization and planning techniques . In the study by Li et al. (2021a), a second-order cone method for a RIES with electricity-gas interconnection was proposed, by which the non-convex electricity-gas optimal power flow model was transformed into a convex optimization model, and the scaling degree was limited by setting a penalty term, which increases the model convergence. In the study by Li et al. (2022a), an energy network was modeled in matrix form, where energy conversion and coupling could be described by the network graph theory, and the complexity of the optimization operation could be reduced. In the study by Li et al. (2021b), a method for the optimal design of RIES partitions was proposed by combining the K-means and genetic algorithm, by which the local resource data were clustered to evaluate the availability of all resources. In the study by Li et al. (2022b), an improved third-order quantized state system (QSS) method was proposed for the simulation of a district heating system (DHS) by combining a quantized state system and timediscrete integration. In the study by Zhu et al. (2022), a dynamic programming algorithm was introduced to replace a multi-stage decision problem with a series of single-stage subproblems to obtain the global optimal scheduling solution. It is pointed out that the aforementioned results investigate the optimization of the RIES from several aspects. Nevertheless, few results study the operating characteristics and solution methods after the distributed energy station (DES) is connected, where the influences of its uncertainty factors on the system are not sufficiently considered.
On the other hand, the uncertainties of various loads and renewable energy sources seriously affect the safe and reliable operation of the RIES (Wang et al., 2020), which is an important problem that cannot be ignored in the optimal operation process of the RIES. In the study by Xu et al. (2022), an efficient wind speed prediction model based on phase space reconstruction and BLS was proposed, which could evaluate the regularity of wind speed effectively without the overfitting issue. In the study by Xiang et al. (2021), an uncertainty model was proposed for the integrated demand response (IDR) of energy prices based on fuzzy and probabilistic variables. In the study by Feng et al. (2018), based on the imprecise Dirichlet model (IDM), a non-parametric fuzzy set of wind power output distribution probabilities was obtained from historical data and used in the uncertainty optimization of the RIES. In the study by Wang et al. (2017), chance-constrained objective planning was structured to model wind power uncertainty, and the developed results could guarantee the reliability of the system under a certain probability of satisfaction. Furthermore, based on the aforementioned probabilistic models (Wang et al., 2017;Feng et al., 2018), the probabilities were dealt with using the point estimation method, and the objective function was then set as the expected value in Li et al.'s (2021a) study. In the study by Li et al. (2022c), an optimal scheduling strategy of an IDR-enabled RIES in uncertain environments was proposed to improve the flexibility of system operation and the comprehensive satisfaction of users, while battery degradation, electric vehicles' interaction, and load uncertainty were ignored. In the study by Fu et al. (2021), a two-stage stochastic programming model was proposed to deal with the uncertainty of PV generation and energy demand in the RIES. However, the aforementioned models deal with uncertainties that are all based on probabilities and rely on the distribution models with uncertain parameters. However, in some cases, the actual operation results have a certain probability of constraint violation; therefore, the robust optimization methods are expected to be used to solve these issues effectively.
Based on the aforementioned discussions, for the RIES containing a DES, it is necessary and important to take the uncertainties of the DES into consideration and thoroughly investigate its optimal and analyzing methods. Hence, in this paper, the robust optimization of the RIES considering the uncertainties of the DES will be discussed in detail. First, the overall model is divided into two layers according to their specific functions: the energy supply network (ESN) is set as the upperlayer model, and the DES is set as the lower-layer model. Second, the lower-layer DES model is formulated as a two-stage robust optimization model in the form of "min-max-min," and it is split into two parts: a main problem and a sub-problem. Then, the column and constraint generation (CCG) algorithm (Zeng and Zhao, 2013) is applied to iteratively solve these two problems in the lower-layer DES model. Finally, considering the correlation between the upper-layer ESN model and the lower-layer DES model, the alternating direction method of multipliers (ADMM)  is used to solve the overall model. The main innovative contributions of this paper are summarized as follows: (1) A two-layer model of the RIES is constructed innovatively according to the functions of each section of the model, where the electricity-gas transmission network model is taken as the upper-layer model and the multi-energy coupled distributed energy station is taken as the lower-layer model. (2) The two-stage robust optimization method is applied to matrix the lower-layer model, which can dramatically enhance the stability of the RIES to resist the uncertainty of the photovoltaic and load. More importantly, it can be easily solved using the CCG algorithm. (3) The ADMM algorithm is used to solve the RIES with distributed energy stations, which can largely reduce the solving difficulty in an iterative way, where a connection between the upper-and lower-layer models is established.
The remaining sections of this paper are arranged as follows: Sections 2 and Sections 3 present the mathematical models of the upper-layer ESN and lower-layer DES which include a two-stage robust optimization equivalent model with its constraints and objective functions, respectively. Section 4 describes the solving methods used to solve the proposed model, whereas the simulation example and obtained results are demonstrated and discussed in Section 5. Finally, the conclusions are drawn in Section 6.

Upper-layer energy supply network model
The ESN is the upper layer of the RIES, which plays the role of connecting various energy-using nodes and supplying energy to the lower-layer DES. Heat is transmitted in the form of water flow, and the energy loss is serious when the transmission distance is long. Moreover, the time lag phenomenon of heat is more notable than that of electricity and gas. Since the heat transmission network has evident disadvantages compared with electricity and gas transmission networks, the RIES studied in this paper only considers electricity and gas transmission networks, while the heat transmission network is not considered.

Steady-state model of the power system
The power flow equation of radial distribution networks is formulated as follows (Farivar and Low, 2013a): where P t j and Q t j denote the active and reactive powers, respectively, injected into node j at time t; P t jl and Q t jl are the active and reactive powers, respectively, flowing from node j to node l at time t; R ij and X ij are the resistance and reactance, respectively, of line ij; I t ij is the current flowing through line ij at time t; n(j) is the set of sending nodes with node j as the end node; and m(j) is the set of the end nodes with node j as the sending node.
In addition, the relationship between the current I t ij flowing through line ij and the node voltage is given as follows: where U t i is the voltage of node i at time t. In the process of transmitting power in the distribution network, reactive power compensation can effectively improve the power factor and reduce line losses. Common reactive power compensation devices include continuous-type Static Var Generators (SVG) and discrete-type capacitors (CP), where CP satisfies where Q t CP is the compensated reactive power of the CP at time t; Q CP,U is the reactive power output corresponding to each step of the CP; θ t CP,m , θ t CP add , and θ t CP del denote the 0-1 variables at time t; γ is the maximum number of steps that can be adjusted; and θ max is the maximum number of adjustments per day.
It should be pointed out that in Eq. 4, a 0-1 variable is set for each grade, and the sum of each grade is the current grade. Equation 5 limits the gears to incremental steps, and it is not possible to have a high gear of 1 and a low gear of 0. Equation 6 ensures that the increase and decrease of gears cannot occur simultaneously, and Eq. 7 shows that the number of gear adjustments per day is limited. Equations 8, 9 together restrict the upward or downward adjustment of gears to the maximum adjustment range of gears.
For the SVG, the model is simpler than the CP and satisfies where Q SVG min and Q SVG max are the lower and upper limits of the regulation range, respectively, and Q t SVG is the reactive power of the SVG at time t.
The input and output energy of a power system should satisfy the law of energy conservation, so the following relationship should be ensured: where Q t L,j and P t L,j are the reactive and active loads of node j at time t, respectively. P t G,j and Q t G,j are the active and reactive powers of the generator of node j at time t, respectively.

Steady-state model of the gas system
The gas system is mainly composed of a gas source, gas transmission pipelines, a compressor, and loads, and the relationship between its flow and pressure is formulated as follows: where f t q is the flow of the pipe q at time t, p t m is the pressure of node m at time t, ς represents the efficiency factor of the pipe q, T 0 is the Frontiers in Energy Research frontiersin.org standard temperature, D q denotes the inner diameter of the pipe q, p 0 is the standard pressure, G is the gas relative density, L q is the length of the pipe q, T qa is the average temperature of the pipe q, and Z qa is the average compression factor of the pipe q.
There is a pressure loss during gas transmission, which needs to be replenished by a compressor. The gas flow to be consumed is formulated as where σ t u is the gas flow consumed by the compressor u at time t; α u , β u , and γ u are the energy consumption coefficients of the compressor u; P t com denotes the gas energy consumed by the compressor u at time t; G u is the coefficient related to temperature and efficiency; f com,u represents the flow through the compressor u at time t; and Z u is the other compression factor.
Similar to the law of conservation of energy in electric power systems, the input and output of gas at each node in a gas system should be equal, which satisfies where f denotes the branch flow vector, ϖ is the node gas injection flow vector, τ is the compressor consumption flow vector, A is the branch node association matrix, U represents the unit node association matrix, and T is the compressor node association matrix. The power flow model and natural gas model is the non-convex and non-linear. In order to facilitate the subsequent solution, second-order cone scaling (Farivar and Low, 2013a;Farivar and Low, 2013b) is used to convert the distribution network and natural gas system models into the convex model.

The objective function of the upperlayer ESN
The operation of the upper-layer ESN is optimized by minimizing net costs as follows: where a, b, and c denote the generators' cost coefficients; h is the cost coefficient of gas supply; G t G is the gas supply quantity of the gas source at time t; price t e and price t g represent the purchase prices of electricity and gas of the lower-layer DES at time t, respectively; P t grid and P t gas are the purchased power and gas power of the lower-layer DES from the upper-layer ESN at time t, respectively; and T is the dispatch period. P t G,k and G t G,k are the k-th generator and air source supply, respectively. N and M are the total number of generators and gas sources, respectively.

Other constraints
In the upper-layer ESN, upper-and lower-limit constraints should be added to the individual equipment output, voltage value, gas flow value, node pressure, and gas supply from the gas source, which are expressed uniformly by where Θ min and Θ max are the lower and upper limits of the variable Θ t , respectively. Θ t represents the variables such as equipment output, voltage value, gas flow value, node pressure, and gas supply from the gas source.
3 Lower-layer distributed energy station model The DES can meet multiple energy demands of consumers and interconnect different types of energy sources through energy conversion equipment. In this study, the considered DES model includes two types of energy inputs, electricity and gas, and three types of energy loads, electricity, heat, and cooling. Moreover, energy conversion equipment includes gas turbine, gas boiler, absorption chiller, and electric chiller, while energy storage equipment includes electricity and heat storage, mainly considering the uncertainty of photovoltaic inside the DES and the uncertainties of the three loads.

Gas boiler
The gas boiler can be modeled according to the following equation: where H t GB and G t GB are the heat production and gas consumption of the gas boiler at time t, respectively, and η GB is the heat production efficiency.

Gas turbine
The gas turbine is coupling equipment that contains electric, gas, and thermal devices with the following input-output relationship: where P t GT denotes the power production of the gas turbine at time t, G t GT is the gas consumption of the gas turbine at time t, H t GT represents the heat production of the gas turbine at time t, η GT is the power production efficiency, and λ is the electricity-to-heat ratio.

Electric chiller
The electric chiller is mathematically expressed as Frontiers in Energy Research frontiersin.org where C t EC and P t EC are the cooling production and power consumption of the electric chiller at time t, respectively. COP EC is the performance coefficient of the electric chiller.

Absorption chiller
The absorption chiller model is presented as follows: where C t AC and H t AC are the cooling production and heat consumption of the absorption chiller at time t, respectively; COP AC is the performance coefficient of the absorption chiller; COP 0 AC is the rated conversion efficiency of the absorption chiller; and a AC , b AC , and c AC denote the conversion factors of the absorption chiller. β AC is the workload rate of the absorption chiller.

Energy storage equipment model
Energy storage equipment is modeled as where E t i is the reserves of energy storage at time t; i is the type of energy storage device, which includes electricity and heat storage; δ i is the self-loss rate of the energy storage device; C t i and D t i are the charging and discharging powers of the energy storage device, respectively; η i,char and η i,dis represent the charging and discharging efficiencies of the energy storage device, respectively; and υ i is a 0-1 variable to represent the storage mode: 0 if discharging, and 1 if charging.
There are also various constraints on the energy storage device, for example, charging and discharging cannot be done simultaneously and the total energy storage remains unchanged in a scheduling cycle, in addition to upper and lower limits on variable parameters, which are formulated as

The model of uncertain parameters
PV and various load data are mainly derived by forecasting; however, there are always various random factors that affect the accuracy of forecasting in reality. Random factors cannot be avoided, and if they are not considered when developing the operation scheme, they may lead to inaccurate operation schemes and even cause operational accidents. To equip the system with a certain antiinterference ability and to ensure the system's safe and reliable operation in different environments, it is necessary to deal with those uncertain parameters. The uncertainty U model of PV output and load is given by where u t PV and u t i,L are the uncertainty variables of PV output and load power at time t after introducing uncertainty, respectively.û t PV andû t i,L are the expected values of PV output and load power at time t, respectively. Δu max ,t PV and Δu max ,t i,L denote the maximum allowable deviations of PV output and load power at time t, respectively; i represents the energy types including electricity, heat, and cooling; and B t PV , D t PV , B t i,L , and D t i,L are 0-1 variables, while Γ PV and Γ i,L are the uncertainty conservativeness adjustment parameters.

Energy balance constraints
The lower-layer DES also needs to satisfy the energy conservation law, and the balance between the input and output quantities of each type of energy should be maintained as where D t e , C t e , D t h , and C t h are the variables of charging and discharging powers corresponding to the electric and thermal energies, respectively, while u t e,L , u t h,L , and u t c,L are the variables of load uncertainty corresponding to the electricity, heat, and cooling, respectively.

Upper-and lower-bound constraints
Like the upper layer, there are also some limit constraints on the output of each device in the lower layer, which are omitted here due to the length of the paper.

The objective function of the lowerlayer DES
After presenting the mathematical models of both layers, the following objective function can be formulated: To facilitate subsequent simplification and solution, the lowerlayer DES model is written in a matrix compact form: The objective function can be split into two stages: the first stage is to optimize the charging and discharging behavior of the power and heat storage equipment, called the main problem (MP), and the second stage is to search for the optimal equipment output and energy purchased under the worst uncertain parameters, called the subproblem (SP). Uncertainty is expressed by the box uncertainty set. The sub-problem is to solve the worst scenario in the interval according to the conservatism parameter and complete the two-stage robust optimization together with the main problem. The main problem and sub-problem of the objective function are described as follows: The sub-problem is a max-min type problem, which is not easy to solve. By the strong dual principle, the inner min problem can be converted into a max problem and merged with the outer one to obtain the following form of the sub-problem (Ding et al., 2016): where ω 1 , ω 2 , and ω 3 are dual variables. In Eq. 32, u T ω 3 is a bilinear variable, which is difficult to solve. Therefore, it is linearized by using the big M method: ], B′ is the auxiliary variable, and M BIG is a large positive real number, which is taken as 10 9 in this study.

The solving method 4.1 The flow of the CCG algorithm
In this study, the RIES is divided into two layers. The upper layer is a traditional ESN model, while the lower-layer DES is a complex two-stage robust optimization model. The CCG algorithm is used to obtain the optimal solution to the lower layer. The process is described as follows: 1) The initial uncertainty parameter u 1 * , the iteration convergence threshold ε, the upper and lower bounds for the iterative solution as LB −∞, UB +∞, and the number of iterations as k 1 are set.

FIGURE 1
Overall flowchart of the proposed solving method.

Frontiers in Energy Research
frontiersin.org 2) The MP is solved as shown in Eq. 30 according to the uncertain parameters u * k , then the optimal solution (x * k , σ * k , y 1 * , y 2 * , /, y * k ) is obtained, and the lower-bound LB σ * k is updated.
3) The storage charging and discharging energy-related variables x * k obtained from the solution of the MP (Step 2) are substituted into the SP as shown in Eq. 33; then, the objective function Q(x * k ) and the corresponding uncertain parameters u k+1 * for the worst-case scenario are obtained, and the upper-bound UB UB, Q(x * k ) is updated. 4) According to the set threshold, it is determined whether the iteration converges; if it satisfies UB − LB ≤ ε, the iteration is stopped and the results are output; otherwise, the variable y k is added with the following constraints: σ ≥ C T y k+1 Dy k+1 ≥ d Ly k+1 + Tx 0 I u y k+1 u k+1 5) We assume k k + 1 and return to step 2.

FIGURE 2
Model of the system used for simulation.

FIGURE 3
Predicted values of various loads and PV.

FIGURE 4
Time-of-use electricity price.
Frontiers in Energy Research frontiersin.org

The flow of the ADMM algorithm
It is worth pointing out that the sum of the objective functions for the upper and lower layers is the complete objective function of RIES optimization, as the two layers are not separate. However, if the upper and lower layers are structured as one model, its complexity will greatly increase, and it will be more difficult to solve it. For this reason, the ADMM algorithm is used to maintain the correlation between the upper and lower layers while solving the models of the upper and lower layers independently. The flow of the ADMM algorithm is described as follows: Step 1: The augmented Lagrange function is established for the upper and lower layers, where it is formulated for the upper layer as For the lower layer, it is where L up and L dn are the augmented Lagrange function of the upper and lower layers, respectively; P t grid and G t gas are the global variables; P t grid,up , G t gas,up , P t grid,dn , and G t gas,dn represent the local variables of purchased electricity and gas in the upper and lower layers; η is the coefficient of penalty term; and λ t 1,up , λ t 2,up , λ t 1,dn , and λ t 2,dn are Lagrange multipliers.

FIGURE 5
Convergence curves of the ADMM and CCG algorithms.

FIGURE 6
Voltage distribution diagram.

FIGURE 7
Gas pressure distribution diagram.

Frontiers in Energy Research frontiersin.org
The lower layer is a two-stage robust optimization model: MP and SP. The MP can solve the electricity and gas purchase, so only the objective function of the MP is changed to the augmented Lagrange function: Step 2 Step 3: Let P t grid P t grid,up,r and G t gas G t gas,up,r are substituted into the lower-layer problem to find P t grid,dn,r+1 and G t gas,dn,r+1 .
Step 4: Let P t grid P t grid,dn,r+1 and G t gas G t gas,dn,r+1 are substituted into the upper-layer problem to find P t grid,up,r+1 and G t gas,up,r+1 .
Step 5: The Lagrange multipliers are updated according to the following equations: Step 6: The maximum deviation value is calculated. When the maximum deviation is less than the convergence threshold, the iteration is terminated and the results are output; otherwise, we return to step 3. The maximum deviation is calculated by φ max P t grid,dn − P t grid,up , G t gas,dn − G t gas,up , where φ denotes the maximum deviation.
It should be remarked that as the iteration proceeds, the Lagrange multipliers gradually become larger. Consequently, the proportion of the difference between the local and global variables in the objective function increases. So, the minimum solution will be made as small as possible until it is close to 0. Because the global variable is constantly taken from the local variable in the upper and lower layers, the difference between the corresponding local variable in the upper and lower layers will become smaller as the iteration proceeds, until it is close to 0, which is the principle of the ADMM algorithm. Because the expected value of the predicted value is not much different from the final solution value, the expectation of the predicted value is used as the initial value of the parameter of the CCG algorithm. The flowchart of the proposed solving method is given in Figure 1 for the overall optimization process.

Example parameters
The arithmetic model selected as a test example is shown in Figure 2. The upper-layer ESN includes a modified IEEE 33-node power system and a Belgian 20-node natural gas system, and the lower-layer DES model is a typical campus energy supply system. It should be pointed out that the method proposed in this study can be applied to different test examples, and only the relevant parameters need to be changed when applying to other test examples.
The forecast curves for each type of load in the lower-layer DES and PV are shown in Figure 3. Without loss of generality, the fluctuation range of the load is 0.1 of the predicted value, and the fluctuation range of PV is 0.15 of the predicted value. Figure 4 shows the time-of-use electricity price, and the gas price is set as 0.35 CNY/ kW. The gas turbine's maximum capacity is 1,400 kW, the maximum capacity of the gas boiler is 800 kW, the capacity of the electric chiller is 500 kW, the absorption chiller's maximum capacity is 300 kW, the maximum storage capacity of electricity storage equipment and heat storage equipment is 1,500 kW h, and the minimum storage capacity is 200 kW h. The maximum value of Heat-power balance.

Frontiers in Energy Research
frontiersin.org charging and discharging power is 600 kW. Other parameters are selected according to references (Ding et al., 2017;Miao et al., 2020).

Convergence analysis
The convergence of ADMM and CCG algorithms is shown in Figure 5.
It can be seen from Figure 5 that both CCG and ADMM algorithms have fast convergence. The CCG algorithm almost converges to 0 after two iterations, while the ADMM algorithm needs nine iterations to obtain the optimal solution. This is because the initial value of the CCG algorithm is the predicted value of uncertain variables, and those predicted values are not very different from the final solution's value. However, the initial value of the ADMM algorithm is set as 0, and the Lagrange coefficient increases with every iteration. At the same time, the difference between the upper and lower layers will become increasingly smaller to obtain the minimum (optimal) solution, which is the reason why the ADMM algorithm is slower than the CCG algorithm. It can be seen from the convergence diagram of the two algorithms that the ADMM algorithm is difficult to converge, which needs to be iterated nine times, while the CCG algorithm has already converged after two iterations. Therefore, the proposed method has high computational efficiency.

Voltage and air pressure analysis
The voltage and air pressure distribution diagrams of the upperlayer ESN are shown in Figures 6, 7, respectively.
As shown in Figure 6, the voltage distribution of the power system of the RIES is between 0.94 and 1.06, which is within the permissible voltage limits and reflects the system's high stability. The time and voltage magnitude relationship at node 1, for example, satisfies the load variations, with lower voltage magnitude during low-load layers and higher voltage magnitude during high-load layers. The voltage distribution corresponds to the generator access at the first hour. For example, the node closer to the generator has a higher voltage value, and the node with a higher distance from the generator has a lower voltage value. Figure 7 shows that the air pressure distribution of the gas system of the RIES is between 8.5 and 9.5 kPa, which is in line with the engineering air pressure requirements. Therefore, the gas system has good stability. Similar to the voltage node, the barometric pressure value also corresponds to the gas source access node, i.e., the closer to the gas source node, the higher will be the barometric pressure value.

Equipment output analysis
The power balance diagrams of electricity, heat, and cooling of the lower-layer DES are shown in Figures 8, 9, 10, respectively.

FIGURE 11
Uncertain parameters in the worst case for Γ 24.

Frontiers in Energy Research
frontiersin.org Figure 8 presents the electrical load of the lower-layer DES which is mainly supplied by the purchased electricity. The gas turbine works at full load during the high-load period to reduce the pressure on the grid. Storage equipment charges during the lowload period and discharges during the high-load period to achieve peak shaving and valley filling of the electrical load, while discharging in time during the period when the electrical load just reaches its peak to avoid self-loss. Figure 9 shows that the heat load mainly relies on the gas turbine, and the gas boiler is relied on only during the time when the gas turbine

FIGURE 12
Uncertain parameters in the worst case for Γ 12. is out of service, which is because the gas turbine produces electricity and heat at the same time. Furthermore, charging and releasing the waste heat from electricity production can effectively reduce costs. Heat storage equipment is not operated because the heat load is small compared to other loads, and the peak-to-valley difference is not so evident. So, there is no high demand for heat storage, and at the same time, the process of charging, releasing, and storing heat will produce losses and increase the costs. Figure 10 shows that the cooling load is mostly dependent on electric chillers to achieve the supply. This is due to the higher cooling factor of electric chillers, as electric chillers can directly convert electricity purchased from the grid to cooling energy, while heat is converted from gas and then converted again to cooling energy which will lead to higher losses. Figures 11, 12 show the worst-case uncertainty parameters for the conservative adjustment parameters Γ 24 and Γ 12, respectively.

Robustness analysis
The sub-problem of the lower-layer DES is used to search for the worst-case scenario with uncertain parameters. Figures 11, 12 show that the worst-case scenarios are all those where the PV output is as small as possible and the load is as large as possible. For Γ 24, the most conservative time occurs, where the worst-case scenario is that the PV output equals to the lower limit of the uncertainty interval, and the load equals the upper limit of the uncertainty interval in 24 h. For Γ 12, the worst case for PV is reached only during 12 h because PV output is only in the daytime and the total output time is less than 12 h, so PV output is the minimum value in the uncertainty interval. For the load, the more output in the peak-energy consumption period, the greater the cost, so the worst case of load when its value reaches the maximum limit in the uncertainty interval appears only during 12 h, while the predicted value is maintained during the remaining hours.
In order to explore the influence of PV's permeability on the results of this study, the following three scenarios are set. Scenario 1: all parameters are the same as the aforementioned example; scenario 2: the penetration level of PV is increased by 30%, and other parameters remain unchanged; and Scenario 3: the penetration level of PV is reduced by 30%, and the other parameters remain unchanged. Table 1 shows the simulation results of the scenarios obtained at Γ 0, Γ 12, and Γ 24 when the worst-case scenario is applied. It can be seen that as the conservativeness parameter increases, the energy purchase cost gradually increases. Meanwhile, the supply deficit gradually decreases, which is due to the fact that the sub-problem of the lower-layer DES can search for the worst scenario and protect it in advance so that the supply deficit gradually decreases. However, this protection also comes at a cost, i.e., the energy purchase cost increases as well. The change in the conservativeness parameter from Γ 0 to Γ 12 increases the cost of energy purchase more and reduces the undersupply less than the change from Γ 12 to Γ 24. This is because the system prioritizes the protection of the worst scenarios (uncertain parameters at maximum load and maximum PV), followed by the protection of the less severe ones. Hence, the increase in the cost of energy purchase slows down with the increase in conservativeness, and the reduction in the undersupply slows down with the increase of conservativeness. It can be seen from Tables 1-3 that the penetration level of PV does not affect the aforementioned results. However, the greater the penetration level of PV, the lower the energy purchase cost of the system. Therefore, the greater the penetration level of PV in the system, the more robustly the system should be optimized.

Conclusion
For reducing the negative impact of uncertainties brought by the DES on the system's safe and reliable operation, a robust optimization model of the RIES considering uncertainties of the DES was developed in this study. First, the RIES with the DES was divided into an upper-layer ESN model and a lower-layer DES model according to their functions. The two layers were iteratively solved using the ADMM algorithm. Second, for the lower-layer DES, a two-stage robust optimization model was presented with the objective function of minimizing the energy purchase cost. In the first stage, the charging and discharging behaviors of storage equipment were explored, and in the second state, the optimal equipment output under the worst case could be developed. Then, for the sub-problem in the second stage, the bilinear terms were reduced to linear constraints using the big M method and solved using the CCG algorithm. Finally, the effectiveness of the proposed solving method was demonstrated by several simulation cases, which showed that the developed results could effectively suppress the uncertainty effects.
In future work, more focus will be given to the contradiction between the conservativeness of the robust optimization model and the optimization results. The robust optimization model based on the probability distribution function of uncertain variables can combine the possible probability distribution function of uncertain parameters COP 0 AC Rated conversion efficiency of the absorption chiller a AC , b AC , c AC Conversion factors of the absorption chiller β AC Workload rate of the absorption chiller δ i Self-loss rate of the energy storage device η i,char , η i,dis Charging and discharging efficiencies of the energy storage device Γ, Γ pv , Γ i,L Uncertainty conservativeness adjustment parameters M BIG Large positive real number, which is taken as 10 9 in this study ε Iteration convergence threshold of the CCG algorithm Variables P t jl , Q t jl Active/reactive powers flowing from node j to node l at time t I t ij Current flowing through line ij at time t U t i Voltage of node i at time t Q t CP Compensated reactive power of the CP at time t Q CP,U Reactive power output corresponding to each step of the CP θ t CP,m , θ t CP add , θ t CP del A 0-1 variable Q t SVG Reactive power of the SVG at time t P t L,j , Q t L,j Active/reactive loads of node j at time t P t G,j , Q t G,j Active/reactive powers of the generator of node j at time t f t q Flow of the pipe q at time t p t m Pressure of node m at time t σ t u Gas flow consumed by the compressor u at time t α u , β u , γ u Energy consumption coefficients of the compressor u T Compressor node association matrix G t G Gas supply quantity of the gas source at time t price t e , price t g Purchase prices of electricity and gas of the lowerlayer DES at time t P t grid , P t gas Purchased power and gas power of the lower-layer DES from the upper-layer ESN at time t P t G,k , G t G,k The k-th generator and air source supply Θ t Other variables such as equipment output, voltage value, gas flow value, and node pressure H t GB , G t GB Heat production and gas consumption of the gas boiler at time t P t GT Power production of the gas turbine at time t G t GT Gas consumption of the gas turbine at time t H t GT Heat production of the gas turbine at time t C t EC , P t EC Cooling production and power consumption of the electric chiller at time t   Maximum allowable deviations of PV output and load power at time t D t e , C t e , D t h , C t h Variables of charging and discharging powers corresponding to electric or thermal energies u t e,L , u t h,L , u t c,L Variables of load uncertainty corresponding to electricity, heat, and cooling ω 1 , ω 2 , ω 3 Dual variables in the lower-layer model B9 Auxiliary variable in the lower-layer model L up , L dn Augmented Lagrange function of the upper and lower layers P t grid , G t gas Global variables of purchased electricity and gas P t grid,up , G t gas,up , P t grid,dn , G t gas,dn Local variables of purchased electricity and gas in the upper and lower layers η Coefficient of the penalty term λ t 1,up , λ t 2,up , λ t 1,dn , λ t 2,dn Lagrange multipliers φ Maximum deviation of the ADMM algorithm Frontiers in Energy Research frontiersin.org