Study on bi‐level multi‐objective collaborative optimization model for integrated energy system considering source‐load uncertainty

The uncertainty of both the generation side and the user side brings certain impacts on the source‐load coordinated operation of integrated energy system (IES). At first, the source‐load characteristics of the IES with large‐scale distributed energy is analyzed. Secondly, considering lower operating cost and higher renewable energy consumption rate, a bilevel optimization model of the source‐load coordination based on the IES operation structure is established, which considers the uncertainty of both the generation side, such as wind and photovoltaic (PV) power output, and the load side. The robust optimization theory and stochastic chance‐constrained programming theory are introduced to establish the upper layer and the lower layer models of the IES. Besides, the firefly algorithm is used to design the solution process of the bilevel optimization model. The results show that the uncertainty of the wind power, PV power, and the load are positively correlated with the comprehensive operation cost. Also, based on the full consideration of the bilevel uncertainties, the system operation cost is significantly reduced. Meanwhile, the improvement of the confidence level of the objective function and constraint conditions will enhance the optimistic value of the comprehensive operation cost of the system, which means gaining higher reliability by sacrificing partial system benefits. Overall, the bilevel optimization model considering the source‐load uncertainty can both improve the renewable energy consumption rate and reduce the comprehensive operation cost of the IES system.

is always affected by the randomness of the generation side, mainly from the renewable energy output, and the load fluctuation, limiting the operation safety and economics of the power grid. 4 Therefore, to optimize the power structure and the system coordination, reducing the uncertainty impacts becomes the one of the most concerns of the IES.
To deal with the uncertainty of the IES, scholars have conducted in-depth research on power generation uncertainty. Literature 5 proposed that the reduction in the uncertainty of renewable energy output in IES can improve the system efficiency. In literature, 6 the uncertainty of the generation side is effectively reduced through the double-layer model. To give full consideration of uncertainty, literature 7 proposed a two-stage robust dispatching model considering the electric boiler, heating storage, energy storage, and the combined heat and power (CHP) units. In literature, 8 based on probability distribution, the stochastic programming and chance constraint rules are introduced to reduce the uncertainty of renewable energy generating. However, too many uncertain variables are involved to accurately determine the probability distribution function, so the load uncertainty is involved in the complementary optimization. In literature, 9 on considering the uncertainty of the system, the influence of demand response, renewable energy output, and multi-load type on the energy system was analyzed. Literature 10 established a comprehensive daily scheduling model considering the bilevel uncertainty of the wind power output and the load.
To reduce uncertainty, scholars have adopted robust optimization, interval number method, and random rule method. Literature 11,12 proposed a distributed robust scheduling model for IES considering the renewable energy output uncertainty. In literature, 13 the stochastic rule is used to optimize scheduling problem considering the uncertainty of renewable energy output. Literature 14 established the optimal scheduling model for the electricity-gas coordinated system considering the randomness of energy output. In literature, 15 a distributed robust optimization method based on fuzzy forecast error set was proposed to solve the electricity-gas coupling collaborative optimal scheduling model. In literature, 16 the distributed robust optimization method was used to deal with the wind power output uncertainty. In literature, 17 the robust method was introduced to optimize the load uncertainty of the IES with electricity-gas collaborative optimization. In literature, 18 the distributed robust optimization method is used to solve the coordinated optimal economic dispatch problem of the IES considering the uncertainty of wind power output. In literature, 19 the distributional data-driven robust optimization (DDRO) model is established to solve the power system scheduling problem with wind power uncertainty. In literature, 20 the standardized modeling method for the interconnected energy hub model is designed. In literature, 21 to reduce the uncertainty from both generation and user side, a stochastic optimal scheduling model was established using robust optimization theory. In literature, 22,23 the dual uncertainty in source-load problem of the power-gas IES is modeled using robust optimization theory, and uncertain set of moment wind power output and load and the stochastic programming model were established. However, by taking the accurate probability distribution information of uncertain factors as premises, stochastic rule method requires complete distribution characteristics of random variables and large calculation scale. By constructing the uncertain parameters sets, robust optimization method does not need to obtain its accurate probability distribution, which is relatively conservative, and the probability and statistics information is hard to obtain.
Based on this, considering the uncertainty of distributed energy output and the load, a bilevel optimization model for the IES source-load coordination is proposed. The specific of this paper is as follows.
1. The multi-objective functions of the upper layer model are to minimize the operation cost of the power subsystem and to maximize the renewable energy consumption rate. The uncertainty of the power subsystem is expressed by intervals. By introducing game theory between the decision makers of both the system and the generation side, the optimal collaborative scheme of the upper layer model is solved through Nash equilibrium. 2. By introducing the stochastic chance-constrained programming theory, the lower layer model describes the objective function and constraint conditions in probability form considering the energy price uncertainty. Then, the firefly algorithm is used to design the solution flow of the upper and the lower layer models. 3. A case study is carried out to verify the security and robustness of the proposed model, which will provide reference for decision makers in benefit-risk weighing.
The rest of this paper is organized as follows. In Section 2, the coupling relations of the subsystems are analyzed. In Section 3, the MTL and LSSVM methods are introduced. In Section 4, a load forecasting model based on the MTL-LSSVM for the electricity, heat, cooling, and gas of the IES is proposed. In section 5, an industrial park in Suzhou is chosen for case study. Finally, Section 6 gives a brief summary and some conclusions of the research paper.

| Analysis of source-load uncertainty
In the IES generation and load collaborative optimization, both the wind power, PV in the generation side and various kinds of load in the user side contain uncertainty, which brings some risks to the source-load cooperative operation system.
(1) Uncertain characteristics of wind power output. In IES, the wind power output is related to the wind speed, whose randomness and volatility performance has uncertain characteristics. 24 Therefore, Weibull distribution is applied to describe the probability density of wind speed, as follows.
Wherein, the wind speed is described as v. The state shape parameter and scale parameter of f WT (v) are described as c and k. According to the generation characteristics of wind power, the relationship between the actual power output and the real-time wind speed is described as follows.
Wherein, P WT (t ) is the available output of the wind turbine at time t. P r,W is the rated power of the /wind turbine. The inlet wind speed, rated wind speed, and cut-out wind speed are described as v in , v r , and v out . The real-time wind speed at time t is described as v.
(2) Uncertain characteristics of PV power output. In IES, the PV generator output is related to the light radiation intensity that contains uncertainty. Therefore, Beta distribution is introduced to fit the probability density of radiation intensity, 25,26 as follows.
Wherein, r is the radiation intensity at time t. r max represents the maximum light irradiance at time t. and are the shape parameters of the Beta distribution, which can be calculated by the mathematical expectation ( ) and the light irradiance variance ( ) at time t.
Thus, the actual PV output power function at time t is.
wherein, P PV (t ) is the maximum PV output at time t. s PV is the light radiation efficiency is PV . the illumination irradiance radiation area is. ( t) is the radiation intensity at time t.
In the source-load cooperative IES, the supply situation will adjust in real time according to demand of electric power, heat, cold, and gas loads of users, which can coordinate and optimize the generation side and the user side. In general, the process of terminal load dynamic change includes random uncertainty and regular change, as follows.
wherein, the power load, heat load, cold load, and gas loads in the system at time t are P electricity (t), P thermal (t ), P cool (t) and P gas (t ). The regular variation functions of each kind of load and time are P electricity base (t), P thermal base (t), P cool base (t), and P gas base (t). ΔP electricity (t), ΔP thermal (t), ΔP cool (t), and ΔP gas (t) represent the random uncertainty of each kind of load.

| Structure of source-load cooperative operation system
According to the uncertainty characteristics of the generator side output and user side load of the IES, this paper considers three modules of the electrical, thermal, and gas energy subsystem. Among them, the power subsystem includes wind power plant (WPP), PV generators, and the external power grid. The demand side management resources are considered in power grid operation in the load side. The heat subsystem includes the external heat network, the heat storage device and waste boiler, and the USES double-effect absorption unit to meet the cooling demand of the user side. The gas subsystem includes the external gas network and the microgas turbine for gas supply. The operation structure is shown in Figure 1.
As shown in the Figure 1, the power subsystem is the core of the whole system to realize multi-energy complementarity and coordination, which can be used as the upper model of collaborative optimization. At the same time, through demand side management, the system benefit can be further improved by realizing cooperative operation with generator side. Considering the uncertainty of distributed renewable energy, the gaming theory is introduced to build the robust optimization model. In addition, based on the optimized operation plan of the power subsystem, the thermal and gas subsystems are optimized as the lower layer collaborative model.

| Source-load cooperative interaction
The bilevel optimization theory has been widely used in the field of decision optimization and is introduced in this paper. In a typical bilevel optimization model, all controllable variables are in different positions. As the leader of the whole model, the decision makers of the upper layer model affect those of the lower layer model and thus determine the overall optimization of the model. 27 The controllable variable formulation of the decision makers in the lower layer model is related to that of the upper layer model. Although partial objective functions of the upper layer model are influenced or even determined by the lower layer model, the decision makers of the upper layer can indirectly affect the objective functions by affecting the strategies of the lower model. The specific mathematical expression of bilayer optimization is as follows.
wherein, the controllable variable of the decision makers in the upper layer model is x. The number of the controllable variables of the lower layer model is N, which are expressed as y 1 , y 2 , ⋯, y N , separately. The objective function of the upper layer model is f( x, y 1 , y 2 , ⋯, y N ), and that of the decision maker j is f j (x, y 1 , y 2 , ⋯, y N ). G(x, y * 1 , y * 2 , ⋯, y * N ) ≤ 0 and H(x, y * 1 , y * 2 , ⋯, y * N ) = 0 are the inequality constraints and equality constraints for the decision makers of the upper layer model. f j (x, y 1 , y 2 , ⋯, y N ) } ≤ 0 and g j (x, y 1 , y 2 , ⋯, y N ) = 0 are the inequality constraints and equality constraints for the decision makers of the lower layer model. According to the above model, x is determined by the decision makers of the upper layer, which can indirectly influence the decisions of the lower layer through x, and thus makes the whole model presents a hierarchical decision-making scenario. Therefore, according to the proposed bilevel collaborative optimization model, the source-load cooperative optimization operation structure can be described as shown in Figure 2.
As shown in Figure 2, the upper layer model and the lower layer model constitute the source-load collaborative optimization operation scheme. Among the operation structure, the power subsystem, as the core of the whole system network, bears the pivotal role of energy conversion and transmission, taking up the key position in the IES source-load coordinated operation optimization. So the power subsystem optimization is taken as the upper layer model, which can directly determine the source-load operation plan. The operation of the other subsystems and their interaction with the power subsystem are determined by the strategy of the decision makers of the upper layer, which will be the lower layer model of the bilevel collaborative optimization model.

| Upper layer collaborative optimization model
To improve the economy and reliability of the IES, the optimized operation plan is made based on the analysis of the bilevel interaction relationship of the source-load coordination optimization and considering different forms of energy, energy supply equipment, and demand side management resources.
(1) Objective function. The upper layer model is the power-charge collaborative optimization module of the power subsystem, and the controllable variables of the decision makers are the distributed generation and load collaborative plan, including the WPP and PV power generation plan, and the demand side management plan. 27 In general, the upper model contains multiple objectives. The subobjective 1 represents the maximum consumption rate of wind power and PV output in the power subsystem, and the subobjective 2 represents the lowest operation cost which includes the operation coast of the source side and the management cost of the load side, as follows.
wherein, f UP,1 is the subobjective 1 of the upper layer model. r is the comprehensive consumption rate of the renewable energy during the operating cycle. T is the number of scheduling periods in the operating cycle. The operation cycle is one day and the scheduling period is 1h. Within the operating cycle, P WT (t) and P PV (t) are the consumption rates of the actual wind power output and PV output at time t, respectively. P WT (t) and P PV (t) represent the forecast ultra-short term wind power output and PV power output at time t, respectively.
wherein, min f UP,2 is the subobjective 2 of the upper layer model. C S is the operation cost of the generation side in the power subsystem, and C L is the management costs of the load side. C IES is the integrated operation cost of the thermal and gas subsystems, which is determined by the indirect impacts of the decision makers of the upper layer model on the lower layer model.

Gas energy subsystem model of source charge interaction in lower layer
Operation scheme of electric energy subsystem Operating cost of gas energy subsystem wherein, P GRID (t) is the purchased power of the power subsystem with the external network at time t and q GRID (t) is the offgrid TOU price at time t. k SB is the operation cost coefficient of the energy storage. P SB (t) is the charge-discharge power of the ESS at time t. When P SB (t) is greater than 0, the ESS is at discharging state, and when P SB (t) is less than 0, the ESS is at charging state. k WT and k PV represent the operation and maintenance cost coefficient of WPP and PV. f MT ( ) is the fuel cost function of microturbine and f FC ( ) represents the fuel cell. P MT ( t) represents the microturbine output at time t and P FC (t ) represents the fuel cell output.
wherein, c 0 delay represents the fixed fee of the translatable load that signed between the IES operator and the translatable load user. P 0 delay ( t) represents the translatable load at time t. c delay is the variable cost coefficient of the translatable load generated through real operating. P delay (t) is the load translation capacity at time t. delay is the conversion coefficient of the user's power consumption loss caused by load translation. c 0 cut is the/ fixed cost of the interruptible load signed between IES operators and users. P 0 cut (t) is the interruptible load capacity at time t. c cut is the variable cost coefficient caused by load translation and load interruption. P cut ( t ) is the interruption power implemented at time t. cut is the conversion coefficient of the user's power loss caused by the load interruption. Meanwhile, considering the demand side management cost and the influence on the user side, the externalities of DSM are contained in the objective function of the upper layer model.
To sum up, the multi-objective problem can be transformed into the comprehensive satisfaction index F through weighted coefficient method. Then, the comprehensive satisfaction index is taken as the objective function of the upper layer model on considering different attributes of the two subobjectives, as shown in Equation (16).
wherein, 1 and 2 are the weighted coefficient of the subobjective 1 and subobjective 2 of the upper layer model. f UP,1,max represents the maximum value that can be achieved when the upper layer model only applies subobjective 1, and f UP,2,min represents the minimum value that can be achieved when the upper layer model only applies subobjective 2.
(2) Constraint conditions. a) Power balance constraints of the power subsystem.
To optimize the power subsystem's source-load collaborative operation, the power balance constraint must be met. Otherwise, the system will suffer power shortage or surplus, and thus the operation reliability of the system reduces. In addition, the power balance constraint of the power subsystem is valid for any period of the operation cycle of the IES, as shown as follows.
And P L (t) can be expressed as follows: Wherein, P L (t) is the load of the power subsystem at time t after DSM. P loss (t ) is the grid power loss at time t. P E,H (t) and P E,G (t ) represent the power converted from power subsystem to thermal subsystem or to gas subsystem, respectively. d( j, t) is the elements of the load shift matrix, which represents the proportion of the translational load, P delay (j), shifted from period j to period t. b) Capacity constraints of renewable power generation. The actual consumption power must not exceed the ultrashort-term predicted output when the wind and PV power consumption plan is made in the upper layer model, which can be expressed as follows.

c) Output constraints of controllable microgenerators
In the upper layer model, the controllable microgenerators' output is supposed to meet the constraints of the output limit and the climbing rate as follows.
wherein, R MT is the rated climbing rate of the microgas turbine. R FC is the rated climbing rate of the fuel cell. Δt is the unit duration. P min MT and P max MT are the minimum and maximum allowable output of the microgas turbine. This constraint formula means that the output of the microturbine must be 0 or meet the range of [P min MT , P max MT ]. Similarly, P min FC and P max FC are the minimum and maximum permissible output of the fuel cell, which means that the output of the fuel cell must be 0 or meet the range of [P min FC , P max FC ]. In addition, as a connection point of the power exchange with the external network, the common coupling point of the power subsystem also needs to meet the power limit constraint, as follows.
wherein, P max GRID is the maximum trading power with the external grid. d) Energy conversion constraints. The sources-load coordination center of the power subsystem is the upper layer of the bilevel optimization model, and the electric-heat and electric-gas conversion formulated by the center directly affect the operation plan of the thermal subsystem and the gas subsystem, which can be described as follows.
wherein, P E,H,max is the maximum conversion power from power to heat and P E,G,max is the maximum conversion power from power to gas. e) Energy storage operation constraints.
Energy storage system (ESS) is a commonly used controllable micropower source. The charging and discharging power and the remaining power in each period of the ESS are the key information when solving the optimal operation scheme, and it can be described as follows.
wherein, S SB ( t) and S SB (t + 1) are the dump power of SB at time t and time t + 1. dis and ch are the discharging and charging efficiency. S min SB and S max SB are the minimum and the maximum value of remained power storage. P min SB is the minimum charge-discharge power of the ESS and P max SB is the maximum. D SB is the self-discharge coefficient of the ESS. Q SB is the ESS capacity. f) Constraints of demand side management.
In the source-load collaborative optimization, both load translation constraints and interruptible load constraints are considered. The DSM constraints are also considered as follows.
wherein, P max E is the maximum power supply capacity of the power subsystem. Equation (27) indicates that the load level of the power subsystem after demand side management cannot exceed the power supply capacity of the system at time t. Equation (28) indicates that the actual translational load cannot exceed the translational load power signed between the IES operators and users at any time. Equation (29) indicates that the actual load interruption power at any period cannot exceed the interruptible load power signed between the IES operators and users.

| Lower layer collaborative optimization model
In the lower layer model, the operation scheme of the thermal and gas subsystems are optimized, and the optimization operation center of the two subsystems are two decision makers whose decisions directly affected by the upper level model. Of course, the value of the objective function will affect the objective function of the upper layer model.
(1) Objective function. The decision strategy of the thermal subsystem optimization operation center is to minimize the operating cost of the thermal subsystem as follows.
wherein, P H (t ) represents the purchased power of the thermal subsystem with the external network at time t and p H (t) represents the heating price. P X (t ) is the charging and discharging power of the heat storage system at time t. When P X (t) is greater than 0, the heat storage system is at the heat-storing state, and otherwise, it is at the heat releasing state. k X is the operation and maintenance cost of the storage device per unit power.
The decision strategy of the gas subsystem optimization operation center is to minimize the operating cost of the gas subsystem as follows.
wherein, Q G ( t) represents the purchased gas from the external gas network at time t. In order to calculate the operating cost of the IES, the natural gas price is converted into the power price. C NG ( t) is the unit natural gas price in spot market at time t and has been converted to the unit power price.
In the thermal subsystem operation optimization, some constraints are considered, including the supply demand balance, the energy storage system operation, the relationship between the heat charging-releasing and the residual heat system, and the heat purchasing with the external heat network as follows.
wherein, C he is the heating coefficient of a double-effect absorption unit, which is related to the output of the microturbine to supply heat. P L,H ( t ) is the heat load of the thermal subsystem at time t, and P H (t ) is the purchased heat power from the external heat network. X(t) and X(t − 1) are the residual heat at time t and time t − 1. x is the self-loss coefficient of the heat storage system. Q min and Q max are the minimum and maximum charging/discharging power of the heat storage system. X min and X max are the minimum and maximum residual heat of the heat storage system. P min H and P max H are the minimum and maximum purchased heat while trading with the external heat network.
While optimizing the gas subsystem operation, some constraints are considered, including the supply demand balance of the gas subsystem, the operation of the gas tank, and gas purchasing with the external network. In addition, to facilitate the cost calculation, the gas amount of the gas subsystem is calculated as power as follows.
wherein, G s ( t) is the releasing power of the gas tank at time t. P E,G ( t ) is the gas generation power of the methane reactor at time t. P L,G (t ) is the natural gas load power at time t. Q s (t) and Q s (t − 1) represent the residual gas power of the gas tank at time t and time t − 1. Q min s and Q max s represent the minimum and maximum residual gas power of the gas storage tank, respectively. G min s and G max s are the minimum and maximum gas releasing power. Q min in and Q max in are the minimum and maximum purchased gas power with the external network.

| Comprehensive collaborative optimization model
According to the upper layer model and the lower layer model established in Sections 3.2 and 3.3, the source-load bilevel collaborative optimization model of the IES is obtained as follows.
wherein, P WT is the real wind power consumption vector and P PV is the real PC consumption vector. P GRID is the power purchase vector from the external grid. P SB is the charge-discharge vector of the ESS. P MT and P FC are the output vector of the microgas turbine and the fuel cell, respectively. P delay is the load translation power vector and P cut is the interruptible load power vector. P E,H is the power transforming vector of electric power to heat, and P E,G is the transforming vector of electric power to gas. P H is the heat purchase vector of the thermal subsystem with the external thermal network. P X is the heat retransmitting vector of the heat storage device. Q G is the purchased gas vector from the external gas network. G s is the charging and discharging power vector of the gas storage tank. Meanwhile, g UP ( ) and h UP ( ) are the equality and inequality constraints of the upper layer model.
among them, P H , P X , Q G , G s are determined by the following solution of the planning respectively g DOWN,H ( ), and h DOWN,H ( ) are the equality and inequality constraints for the decision makers of the lower layer model.

| The problem description
Based on the bilevel comprehensive source-load optimization model, the uncertainty of the wind power output, PV output, and energy price are considered, and the robust optimization theory and the stochastic opportunity-constrained programming are introduced in the double-layer model, as the uncertainty modeling described below.
(1) Robust optimization in upper layer model. The upper model objective functions of the power subsystem cooperative operation optimization are greatly affected by the distributed wind power and PV. In the upper layer model, the research results may be as expected under extremely adverse circumstances while the control variables are formulated without considering the distributed energy output uncertainty. In the robust optimization, the anthropomorphic contribution of the distributed energy is regarded as another decision maker whose control variable is selected with the aim of obtaining the worst result. Meanwhile, the decision makers of the upper layer model is regarded as the system decision maker, which is opposite to the decision makers of the distributed energy generation, so that the objective function can get reach the optimal. Moreover, the two control variables are coupled and influenced mutually, forming a dynamic game process. The system decision maker wants to optimize the objective function as much as possible based on the renewable energy output strategy, which is a typical robust optimization problem.
In the upper layer model, the robust optimization is obtained considering the introduction of decision makers of distributed energy, as shown in Equation (35). The formula is a constrained min-max optimization problem, in which the uncertainty of the renewable energy output is considered and an optimal strategy of minimizing the system cost is constructed.
wherein, u 1 is the control variables for system decision makers and u 2 is the control variables for distributed energy decision makers. F( u 1 , u 2 ) represents the system decision maker and is the objective function of the decision makers in the upper layer model to minimize the integrated operation cost, using u 1 as the only control variable. Similarly, F( u 1 , u 2 ) is the objective function of renewable energy output decision makers, which is to maximize the integrated operation cost, using u 2 as the only control variable. g UP (u 1 , u 2 ) = 0 is the equality constraint of robust optimization model, and h UP (u 1 , u 2 ) ≤ 0 is the inequality constraint. U 1 represents the feasible strategy set for system decision makers, and U 2 represents the feasible strategy set for the decision maker of distributed renewable output.
(2) Stochastic chance-constrained programming in lower layer model.
In the lower layer model, the energy price curve is uncertain because of the spot price mechanism of the external thermal and natural gas networks. In view of the uncertainty caused by the above energy prices, the stochastic chanceconstrained programming theory is introduced in this section to describe the uncertainty of the lower layer model, as shown in Equation (36). The objective function of decision makers in the lower layer model is to minimize the operating cost.
wherein, u is the control variables of stochastic chanceconstrained programming model. Pr{ } is the corresponding probability of an event. f( u, ) is the objective function, and is stochastic variable. is the confidence level of the objective function when the objective function is less than f . j is the confidence coefficient of the constraint j. p is the number of constraints with uncertain form. f is the optimistic values when the confidence coefficient of f( u, ) is . q is the number of constraints with certain form.
Since two decision makers are in the lower layer model, some differences in the variables of the stochastic opportunityconstrained programming models are produced. In the thermal subsystem optimization, u = [ P H , P X ] and = p H , in which p H is the energy price vector. In the gas subsystem optimization, u = [ Q G , G s ] and = C NG , in which C NG is the energy price vector.

| Objective functions
(1) Objective function of the upper layer model.
In the upper layer model, different modeling forms are applied based on the deterministic forms. Both system decision makers and distributed energy output decision makers are considered. Therefore, based on the robust optimization upper layer model, the objective function of the system decision maker is obtained.
The objective function of distributed energy decision makers is as follows.
Among them, the specific variables of u 1 and u 2 are shown in Equations (39) and (40).
(2) Objective function of lower layer model. Based on the stochastic chance-constrained programming, the objective function of the lower layer model is expressed through confidence coefficient. The objective function of the decision makers of the thermal subsystem is as follows.
wherein, H represents the confidence level of the thermal subsystem objective function, which is related to the risk preference of the decision makers. C H represents the optimistic value of the objective function under the confidence level H .
Similarly, the objective function of the decision maker in the gas subsystem can be expressed as follows.
wherein, G is the confidence level of the objective function of the gas subsystem and C G is the optimistic values of the objective function under the confidence level G . (1) Constraints of upper layer model. a) Constraints for distributed energy decision makers. Considering the gaming behaviors between the distributed energy decision makers and the system decision makers, the strategy of distributed energy decision makers needs to meet the strategy constraint set, that is, the interval of distributed energy output. Among them, the randomness of distributed PV output is described as intervals in the robust optimization as follows.

| Constraint conditions
wherein, P PV (t) is the actual distributed PV output and P PV (t) is the forecasted PV output of the IES at time t. P PV (t) is the deviation between the actual and the forecasted PV output at time t. P PV,min (t) and P PV,max (t) represent the minimum and maximum deviation, respectively.
Similarly, the randomness of the distributed wind power output is described in the form of intervals as follows.
wherein, P WT (t) is the actual distributed wind power output and P WT (t) is the forecast wind power output at time t. P WT (t) is the deviation between the actual and forecast wind power output at time t. P WT,min (t) and P WT,max (t) represent the minimum and maximum deviation, respectively.
Similarly, due to load uncertainty, the interval form can be described as follows.
wherein, P L (t) is the actual load and P L (t) the forecasted load of the IES at time t. P L (t) is the deviation between the actual and the forecast load at time t. P L,min (t) and P L,max (t) represent the minimum and maximum deviation, respectively. b) System decision maker constraints. Based on the IES upper layer model, the constraints are made considering the strategies of system decision makers.
(37) min (40) (2) Constraints of lower layer model. By introducing the random chance constraint programming modeling in the lower layer model, the power balance constraint is hard to achieve. The power balance constraint in the model should be satisfied when the unbalanced power deviation satisfies within a certain range and the probability satisfies the expected confidence.
The power balance constraint of the thermal subsystem in probabilistic form can be expressed as follows.
wherein, is the allowable deviation. H is the confidence coefficient when the thermal subsystem decision maker constraints are met.
The power balance constraint of the gas subsystem in probabilistic form can be expressed as follows.
wherein, G is the confidence coefficient when the gas subsystem decision maker constraints are met.

| Solving method
According to the robust optimization of the upper layer model, the system decision makers and distributed energy decision makers finally reach the Nash equilibrium through gaming. Assuming that the Nash equilibrium of the upper robust model is (u 1 , u 2 ), neither of the two kinds of decision makers can optimize their objective function further by changing strategies. Then, the formula is satisfied by ∀ (u 1 , u 2 ).
It can be seen from the above two equations that, based on the Nash equilibrium, the system decision makers cannot reduce the comprehensive satisfaction further by changing strategies, which also indicates that the distributed energy decision makers cannot improve their comprehensive satisfaction either. Therefore, supposing when the strategy of the system decision makers is u 1 , the response of the distributed energy decision makers to u 1 is u 2 (u 1 ), so it is F( u 1 , u 2 ( u 1 ) ) that needs to optimize, which can be translate into a single decision maker model. However, if the response of the system decision maker to strategy u 2 is u 1 (u 2 ), there is a fixed point in the gaming process, u 1 = u 1 (u 2 (u 1 ) ), which is the Nash equilibrium point.
In the lower layer stochastic opportunity optimization model, the monte Carlo stochastic simulation technique is introduced to solve the objective function and the constraint conditions through the probability form, which is expressed as follows.
wherein, is the random variable of the probability space (Ω, A, Pr). The sample s is chosen randomly from, Among Ω, s = 1, 2, ⋯, N s , the realization value f( x, y 1 , y 2 , ( s ) ) of f( x, y 1 , y 2 , ) at the amount of N S is obtained. The calculating steps of U 1 with the stochastic simulation are as follows.
1. Setting the counter at n sum = 0 and s = 1; (52) P GRID (t) + P MT (t) + P FC (t) + P WT (t) + P PV (t) + P SB (t) = P L (t) + P loss (t) + P E,H (t) + P E,G (t) 0 ≤ P WT (t) ≤ P WT (t);0 ≤ P PV (t) ≤ P PV (t) P MT (t)(P MT (t) − P min MT ) ≥ 0;P MT (t) ≤ P max MT P FC (t)(P FC (t) − P min FC ) ≥ 0;P FC (t) ≤ P max 2. s is produced stochastically from Ω based on the probability measure Pr In this step, several sets of data s is simulated in random distribution included in U 1 . 3. For each set of basic data, the function value of the constraint condition f(x, y 1 , y 2 , ( s ) ) is calculated based on the lower layer optimization model in definite form. If the value of the function is less than 0, the constraint holds for this set of data. Then, set the counter at n sum = n sum + 1 and proceed to the next step. Otherwise, the constraint condition is not satisfied, and go to the next step with no counting. 4. Judge the authenticity of s < N s . If it is established, count s = s + 1 and go back to step 2. Otherwise, proceed to the next step. 5. n sum ∕N s is the confidence coefficient satisfying the constraint conditions. If the confidence coefficient reaches the setting value, the constraint condition should be reasonably established.
The calculating steps of U 2 with the stochastic simulation are as follows.
1. s is produced stochastically from Ω as the basic data of the simulation. s = 1, 2, ⋯, N s . 2. Calculating the value of function h s = f(x, y 1 , y 2 , ( s ) ), s = 1, 2, ⋯, N s . The integrated operation cost data {h 1 , h 2 , ⋯, h N s } of N s can be obtained. 3. Set N ′ as the integer part of N s and represents the confidence level at which the objective function follows. 4. Output the minimum value at N ′ of {h 1 , h 2 , ⋯, h N s }, which is the desired value and the optimistic value of the lower layer model at the confidence level .
To sum up, the firefly algorithm is introduced to solve the bilevel uncertainty optimization model. The algorithm was first designed by Professor Xin-she Yang of the Cambridge University. In the firefly algorithm, the optimal object is the firefly population. Each firefly needs to meet some basic characteristics, namely, one is that any two fireflies will attract each other within the range of mutual perception. The other is that the fireflies with higher fluorescein value can attract more other fireflies. The attraction degree is related to distance. Farther the distance is, less the attractive there will be. [28][29][30] In the optimal firefly population, the operation of the individual firefly needs to meet some basic characteristics. Within the mutual perception range, any two fireflies will attract each other, which enhances the moving tendency between them. Then, fireflies with high luciferin value can attract more other fireflies to move toward themselves, but the attraction degree is related to the distance. The farther the distance is, the less the attraction becomes. Therefore, the solution process of the specific bilevel model is shown in Figure 3, and the pseudo-code of firefly algorithm is shown in Table 1.

| CASE STUDY
To verify the feasibility and the accuracy of the proposed model, a typical day of an IES is selected to analysis. The operation structure of the system follows Figure 1. Figure 4A shows how the distributed PV and wind power output of an IES changes in a typical day, and Figure 4B shows the loads conditions of the power load, the heat load and the gas load. The time-of-use (TOU) tariff mechanism is applied in the energy transaction between the power subsystem and the external power grid, as shown in Table 2.

| Basic data
C NG , the gas price, equals 2.06 yuan/m 3 . The load interruption compensation coefficient ( ) is 3. The energy storage discharging efficiency ( dis ) and charging efficiency (

| Scenario settings
Considering the renewable energy output uncertainty and the load distribution, four scenarios are set to uncertainty have an important impact on system operation cost and renewable energy consumption. Therefore, four different scenarios are set up to analyze the source-load collaborative optimization operation plan, as shown in Table 3. As shown in Table 3, in scenario 1, the source-load uncertainty is not considered. In scenario 2, the uncertainty of the source side is considered. In scenario 3, the load uncertainty of various users is considered. In scenario 4, the dual uncertainty in both source and load side is considered.

| Collaborative optimization results of scenario 1
According to the bilevel optimization model, the collaborative optimization results of scenario 1 are obtained, as shown in Figure 5.
As shown in Figure 5, the coordinated optimization operation scheme of the IES in scenario 1 can better meet the energy supply demand of the system. From period 1 to period 7, due to the lower load level of the power subsystem and the lower TOU price of the external grid, the system purchases power from the external grid to store by the ESS for discharge during the peak load periods. During this time, the load is mainly met by the external grid with lower power price, and the microgas turbine and fuel cell remain stopped. The heat energy subsystem from the first period to the seventh period is the peak load period. The heat energy subsystem mainly meets the heating demand by purchasing from the external thermal network and heat storage devices. Period 8 to 15 comes the first load peak. Since the renewable energy output reached the peak, the power supply mainly consists of distributed wind power and PV. Due to the higher TOU power price of the external grid, the power subsystem could gain profits through the surplus power and the fuel cell with relatively lower generation cost operates more. During period 16th to 24th, the second load peak appears. The optimal operation of power subsystem increases the output of the fuel cell. Microgas turbines are scheduled for output, and the energy storage equipment is in full discharge state. However, with higher TOU price, the system prefers to meet the power demand through other forms, instead of purchasing power from the external grid. As the heat load reaches the peak, the system increases heat purchasing from the external thermal network. To sum up, in scenario 1, the expected comprehensive operating cost of the IES is 1209040 yuan, and the expected consumption rate of wind power and PV reaches 98.235%.

| Collaborative optimization results of scenario 2
According to the bilevel optimization model, the collaborative optimization results of scenario 2 are obtained, as shown in Figure 6.
It can be seen from Figure 6 that compared with scenario 1, to cope with the uncertainty of distributed generation output, the controllable micropower output in scenario 2 is higher than that in scenario 1. For the upper layer model, considering the uncertainty of distributed power output, relative strategy is given out to increase the comprehensive operation cost and reduce the wind power and PV curtailment. Meanwhile, to avoid the worst impact, the IES increases the output of controllable microgenerators and power purchasing from the external network, and curtail partial wind and PV power to meet the peak load regulation capacity requirements of the system. In conclusion, in scenario 2, the expected comprehensive operating cost of the IES is 12880.910 yuan, and the expected comprehensive consumption rate of wind and PV power reaches 97.314%.
As shown in Figure 7, compared with scenario 1, the expected operating cost of the IES is higher in scenario 3 due to the load uncertainty, which is mainly due to the system sacrificing economics to reduce the negative impact of uncertainty. However, through the uncertainty simulation, the actual average operating cost of scenario 3 is significantly lower than that of scenario 1, about 6.30% from 13611.080 yuan to 12753.010 yuan. In conclusion, in scenario 3, the expected comprehensive operating cost of the IES is 13342.560 yuan, and the actual comprehensive operating cost of wind and PV is 12753.010 yuan, and the comprehensive consumption rate reaches 96.658%.

| Collaborative optimization results of scenario 4
According to the bilevel optimization model, the collaborative optimization results of scenario 4 are obtained, as shown in Figure 8.
As shown in Figure 8, the expected operation cost of the system in scenario 4 is 1384660 yuan, higher than that of other three scenarios due to the consideration of the dual uncertainty on both power supply side and demand side. However, considering the dual uncertainty, the actual operating cost of the IES is 12464.330 yuan. The result of scenario 4 is better than that of scenario 2 and scenario 3. Compared with the result of scenario 1, without considering the uncertainty, scenario 4 decreases by 8.43%, which effectively verifies the correctness and effectiveness of the proposed model. In conclusion, under scenario 4, the expected comprehensive operation cost is 1 384 660 yuan, the actual average comprehensive operating cost is 12 464.330 yuan, and the expected comprehensive consumption rate of wind power and PV reaches 92.486%. Therefore, the expected comprehensive operation cost, the actual comprehensive operation cost, and the consumption rate of wind power and PV can be obtained from the above four different scenarios. Among them, the expected operation cost and the wind power and PV consumption rate of the IES are obtained through the model. The actual comprehensive operation cost and actual wind power and PV consumption rate are obtained by simulating for 20 times with uncertainty and then taking the average value based on the system operation plan. The specific comparison results are shown in Table 4.
It can be seen from Table 4 that considering the system uncertainty, the IES operation plan adapts more for the actual operation conditions, which can significantly reduce the operation cost and improve the wind and PV consumption rate.

| Uncertainty analysis
Considering the source-load uncertainty in scenario 4, the uncertainty and index sensitivity of the IES is further analyzed. The game equilibrium analysis is carried out for the upper layer robust optimization model, and the confidence level analysis is conducted for the lower layer stochastic chanceconstrained programming model.

| Model uncertainty analysis
The game process between decision makers of the system and renewable energy in the upper robust optimization model is shown in Figure 9.
It can be seen from Figure 9 that the decision makers of the system and renewable energy gradually converge to the Nash equilibrium state by iterating. In Nash equilibrium, neither system decision makers nor renewable energy decision makers can make better game payment by changing strategies. Instead, the decision maker of the system formulates the operation strategy by considering the worst-case uncertainty of the system, making the optimal operation scheme of the upper layer model more suitable for the actual operation conditions.
The relationship between the comprehensive operation cost and the model confidence in the lower layer stochastic chance-constrained programming model is shown in Figure 10.
As shown in Figure 10, with the improvement of the confidence level of the objective function and constraint conditions, the optimistic value of the comprehensive operation cost of the lower layer model is closely related to the confidence degree of the objective function and constraint conditions. With the gradual improvement of the confidence level, the optimistic value of the system comprehensive operation cost also rises, which shows that when the confidence level increases, the conditions for the optimal operation scheme to meet the objective function value and constraint conditions become more stringent, resulting in the system optimization operation plan tend to be conservative and sacrifice certain economy.

| Sensitivity analysis of uncertain index
(1) Impact of output uncertainty on the source-load cooperative system.
In the upper layer robust optimization model, the uncertainty of the distributed wind power and PV output is described through intervals, which is directly related to the game strategy set range and the economic index of the upper level model. Therefore, considering the dual uncertainty on both source side and load side, the relationship between the comprehensive renewable power consumption rate and the operation cost is analyzed according to the different ranges of distributed wind power and PV output, as shown in Figure 11.  It can be seen from Figure 11 that the range of distributed wind and PV output affects the comprehensive consumption rate index and the comprehensive operation cost of the upper layer. With the increase of the uncertainty interval range, the comprehensive wind and PV power consumption rate decreases, and the comprehensive operation cost of the upper layer increases. In fact, by fitting the two curves in the figure, the sensitivity of the comprehensive wind and PV power consumption rate to the interval range is −0.0869/kw, and the sensitivity of the comprehensive operation cost in the upper layer to the interval range is 13.136 yuan / kW.
(2) Influence of load uncertainty on source-load cooperative system.
Considering the dual uncertainty of the source and load sides in scenario 4, the relationship between the comprehensive wind and PV power absorption rate index of the system and the comprehensive operation cost and the interval range of the upper layer is analyzed according to the different interval ranges of the system load level, as shown in Figure 12.
As shown in Figure 12, the load range also affects the comprehensive wind and PV power consumption rate index and the upper comprehensive operation cost. With the increase of the uncertainty interval, the comprehensive wind and PV power consumption rate decreases, and the comprehensive operation cost of the upper layer increases. In fact, through fitting, the sensitivity of the renewable power consumption rate to the load range is −0.1274/kw, and the sensitivity of the comprehensive operation cost in the upper layer to the interval range is 16.592 yuan / kW. Compared with the distributed power output interval, the sensitivity of the comprehensive operation cost is higher, mainly due to the large load capacity in the IES.

| CONCLUSION
Considering the impact of source-load uncertainty on the collaborative optimization operation of the IES, the operation cost can be effectively controlled and the utilization rate of renewable energy can be increased. Based on the establishment of the bilevel collaborative optimization model and the combination of the firefly algorithm, the following conclusions are obtained.
1. Comparing with the scenario without the uncertainty consideration, the uncertainty-involved collaborative operation plan is more suitable for the actual operation conditions and can significantly reduce the actual operation cost of the system. 2. During actual operation of the IES, on considering the source-load uncertainty, appropriate models and methods should be chosen based on the specific. By dealing with the generation uncertainty and load randomness, the introduction of the robust optimization and the stochastic rule method makes the results more instructive. 3. The model confidence analysis shows that the improvement of the confidence level of the objective function and constraints will enhance the optimistic value of the comprehensive operation cost of the system, that is, higher reliability can be obtained through partial economic expense. 4. The sensitivity analysis shows that larger ranges of the distributed wind power, photovoltaic, and load levels always lead to lower renewable energy consumption rate and the higher comprehensive operation cost of the IES.