An Energy Optimal Dispatching Model of an Integrated Energy System Based on Uncertain Bilevel Programming

An integrated energy system (IES) involving a large number of decision-makers causes problems of bad coordination between energy sub-networks and the IES and it is not able to fully consider the multi-energy complementarity among multiple decision-makers. In this context, firstly, this paper constructs an energy optimal dispatching model of an IES based on uncertain bilevel programming. The upper model takes the transformation matrix of energy hubs as the upper decision-maker, taking the minimum operation cost of the IES in the form of confidence as the objective function; the lower model takes each optimal operation plan of the electric power sub-network, the thermal energy sub-network, and the gas energy sub-network as the lower decision-makers, aiming at the operation economy of each sub-network and considering their operation as necessary constraints. Secondly, a firefly algorithm with chaotic search and an improved light intensity coefficient is designed to improve the proposed model. An empirical analysis was conducted on a pilot area of an integrated energy system in Hebei Province. The results show the following: (1) The typical daily operating cost of the integrated energy system in winter is lower than that in summer; (2) under the same load level, the typical winter and summer running costs of the integrated energy system are lower than that of the traditional microgrid; (3) compared with the particle swarm optimization algorithm, the improved firefly algorithm proposed in the paper has obvious advantages both in terms of running cost and solution time; and (4) when the confidence of the objective function and the constraints increases, the operating cost of various schemes also increase. The above results validate the effectiveness of the energy optimal dispatching model of the IES and the economy of the system operation under the multiple decision-maker hierarchy.


Motivation
With the rapid development of the national economy, the energy consumption system mainly based on fossil energy has brought about increasingly serious environmental issues [1]. Traditional energy systems subordinate to different departments for management and operation, such as cooling, heat, and power, cannot play a role in multi-energy coordination and optimization, which reduce the

Literature Review
Many scholars have studied the multi-energy cooperative and complementary optimization of IESs from different perspectives, as can been shown in Table 1. Wang et al. [8] proposed a cooperative and complementary operation scheme for distributed cooling, heating, and power supply systems in multiple regions by analyzing the disadvantages of the distributed combined heat and power system (distributed energy system (DES)/combined cooling heating and power system (CCHP)). Jiang et al. [9], considering the IES of photovoltaic, wind, natural gas, and power coupling and complementarity, and aiming at the lowest operating cost, constructed a mathematical model of energy equipment and energy flow, and proposed a corresponding optimal dispatching strategy. Wang et al. [10] established physical and mathematical models of a regional energy system including a photovoltaic energy storage system with the aim of minimizing operating costs and proposed corresponding optimal dispatching strategies. Gao et al. [11] studied the integrated modeling of a power system, a thermodynamic system, and a natural gas system coupling, and proposed a control strategy for collaborative optimization. Jiang et al. [12] believed that partial renewable power fluctuation on a power network is transferred to a gas system and a cooling or heating system by the coordinated operation of multi-CCHP. Li et al. [13] proposed a Chance-constrained programming CCP-based dispatching model, and its solution approach was proposed for isolated microgrids (MGs). Li et al. [14] proposed a collaborative operation optimization model of the bilateral cooperation of a system and its users in an integrated multi-energy system (IMES) built by an IMES operator. Wang et al. [15]. proposed an original multi-objective optimization model for IES operation formulated to minimize the operating cost, primary energy consumption, and carbon dioxide emission of IESs and to optimally reduce the power loss and voltage magnitude deviation of the utility grid. Wu et al. [16], on the basis of insufficient consideration of coordination and the mutual benefit of multiple energy supply regions and multi-energy flow dispatching, proposed a user-side energy Internet planning method based on bilevel programming. Zhang et al. [17] proposed a bilevel optimal dispatching model for a power-natural gas IES considering the optimal operation of a natural gas system. However, existing studies have considered only single decisions, and the significance of using bilevel programming theory is only in capturing key decision variables. However, with the orderly advancement of China's electricity marketization reform, the IES involves multiple operating entities and is a strongly coupled system that integrates energy flow, information flow, and business flow.
For a multi-energy complementary collaborative optimization model of an IES, mathematical programming methods and intelligent optimization algorithms are usually used to improve the operation benefit of the whole system. Sun et al. [18], aiming at minimizing the daily operating costs of operation and maintenance costs, power purchasing costs, fuel costs, and energy storage depreciation costs, constructed a multi-energy collaborative optimization model considering various working modes of ice storage air conditioning, and the model was solved by using mixed integer linear programming (MILP). Zhao et al. [19] proposed a three-level collaborative optimization model of cooling-heat-electric power supply with optimal equipment selection, capacity, and operation parameters solved by a particle swarm optimization (PSO) algorithm. Xu et al. [20], aiming at the lowest daily operating cost, proposed a collaborative optimization model of plant IESs considering cascade utilization, which was independently modeled under the constraints of cooling, heating, power

Our Contributions
Based on the orderly progress of China's electricity market-oriented reform, an IES involves a variety of operators, which comprise a strong coupling system of integrated energy flow, information flow, and business flow. It is suitable for modeling research using uncertain bilevel programming theory. Therefore, with the aid of uncertain bilevel programming theory, this paper establishes an energy optimal dispatching model of an IES with the randomness of the system considered. The main innovations are as follows.
(1) The structure of an IES is designed, and an optimal energy dispatching model of the IES based on uncertain bilevel programming theory is proposed. (2) The energy optimal dispatching model of the IES based on uncertain bilevel programming is established. The upper model takes the transformation matrix of energy hubs as the upper decision-maker, and takes the minimum operating cost of the IES in the form of confidence as the objective function; the lower model takes the optimal operation plan of the power grid, heat network, and gas network as the lower decision-maker, and takes the operation economy of each subnet as the objective. Considering the necessary constraints of each energy subnet operation, the interaction between the upper model and the lower model is determined. (3) Based on chaotic search and adaptive light intensity coefficient, an improved firefly algorithm is proposed to construct the optimal IES energy scheduling model for uncertain bilevel programming. (4) A pilot area in Hebei Province is selected for the IES. Two operation schemes of the IES for typical winter and summer days are set up, and the operation cost of the IES under different confidence levels is analyzed, which helps to bring the economy of system operation into full play.  It can be seen in Figure 1 that there is a clear energy flow relationship between different energy subnets of the IES. A two-level planning model is a special case of a multi-level planning model. The upper level is a decision-maker, and the lower level is a subordinate. There is an interactive relationship between them, taking into account the reflection of the lower layers. At the same time, since the IES contains multiple subnets, different individual interests with different goals and decision variables need to be met during the energy scheduling process. Therefore, in the IES energy optimization scheduling model based on uncertain two-level planning, the upper decision-maker is the center for formulating the energy conversion plan of the IES, and the decision variable is the energy pivot conversion matrix of the IES (i.e., the difference between different energy transformation plans for each time period); the decision-makers at the lower levels are the energy subnets, namely the electric energy subnet, the thermal energy subnet, and the gas energy subnet. Specifically, the information flow relationship of the energy optimization scheduling model of the integrated energy system for the bilevel planning is uncertain, as shown in Figure 2. It can be seen in Figure 1 that there is a clear energy flow relationship between different energy subnets of the IES. A two-level planning model is a special case of a multi-level planning model. The upper level is a decision-maker, and the lower level is a subordinate. There is an interactive relationship between them, taking into account the reflection of the lower layers. At the same time, since the IES contains multiple subnets, different individual interests with different goals and decision variables need to be met during the energy scheduling process. Therefore, in the IES energy optimization scheduling model based on uncertain two-level planning, the upper decision-maker is the center for formulating the energy conversion plan of the IES, and the decision variable is the energy pivot conversion matrix of the IES (i.e., the difference between different energy transformation plans for each time period); the decision-makers at the lower levels are the energy subnets, namely the electric energy subnet, the thermal energy subnet, and the gas energy subnet. Specifically, the information flow relationship of the energy optimization scheduling model of the integrated energy system for the bilevel planning is uncertain, as shown in Figure 2.  As can be seen in Figure 2, in the energy optimal dispatching model of an IES based on uncertain bilevel programming, the control variable of the upper model is the energy conversion matrix of the IES. The lower decision-makers make their own decisions according to the decision of the upper model (i.e., the optimal operation plans of each energy subnet), feed their operation costs back to the upper model, and modify their own control variables according to the feedback.

The Upper Model
In the energy optimal dispatching model of an IES based on uncertain bilevel programming, the upper decision-maker's goal is to minimize the comprehensive operating cost, which is affected by the objective function of the lower model, and its cost includes the loss of energy hub transformation. Therefore, the objective function of the upper decision-maker is not directly determined by its own decision variables, but by its own decision variables, so as to influence the decision-making of the lower model, thus indirectly determining its own objective function.
The energy transfer equation in the form of an energy hub in the t-th dispatching period of the IES is where ( ) C t is the energy coupling matrix at time t. ( ) P t is the energy input vector at time t. ( ) L t is the energy load vector at time t. Practically, there is some energy loss when each energy subnet of the IES performs energy conversion and transmission, which can directly influence the economy of the energy interconnected As can be seen in Figure 2, in the energy optimal dispatching model of an IES based on uncertain bilevel programming, the control variable of the upper model is the energy conversion matrix of the IES. The lower decision-makers make their own decisions according to the decision of the upper model (i.e., the optimal operation plans of each energy subnet), feed their operation costs back to the upper model, and modify their own control variables according to the feedback.

The Upper Model
In the energy optimal dispatching model of an IES based on uncertain bilevel programming, the upper decision-maker's goal is to minimize the comprehensive operating cost, which is affected by the objective function of the lower model, and its cost includes the loss of energy hub transformation. Therefore, the objective function of the upper decision-maker is not directly determined by its own decision variables, but by its own decision variables, so as to influence the decision-making of the lower model, thus indirectly determining its own objective function.
The energy transfer equation in the form of an energy hub in the t-th dispatching period of the IES is where c ij (t) is the energy coupling factor at time t, representing the proportion of energy from the j-th kind of energy network converted and input into the i-th kind of energy network at time t, and i, j ∈ E = e, h, g , E is the energy subnet set of the IES. The second equation is similar to the first.
Variable subscripts e, h, and g are, respectively, the electric power subnet, the thermal energy subnet, and the gas energy subnet. P i (t) is the input power of the i-th kind of energy network at time t. L i (t) is the load power of the i-th kind of energy network at time t. We change Equation (1) into a matrix equation form, as Equation (2).
where C(t) is the energy coupling matrix at time t. P(t) is the energy input vector at time t. L(t) is the energy load vector at time t. Practically, there is some energy loss when each energy subnet of the IES performs energy conversion and transmission, which can directly influence the economy of the energy interconnected microgrid operation. Therefore, the upper decision-maker needs to consider the efficiency of energy conversion. We decompose c ij (t) into where v ij (t) is the allocation factor at time t, representing the proportion of energy from the j-th kind of energy network converted into the i-th kind of energy network at time t. η ij (t) is the efficiency factor, representing the efficiency of the j-th kind of energy being converted into the i-th kind of energy. When η ij (t) is known, the upper decision-maker's decision variables are considered as v ij (t). However, η ij (t) changes as the operation conditions of the IES and time change, so this paper still considers the upper decision variable as c ij (t), and takes the influence on η ij (t) by the lower decision into account without changing the essence of the model. In the upper model, the total daily operating cost of the IES is shown in Equation (4). The first item is the total daily operation cost of each energy subnet, and the second item is the loss cost of energy conversion of the IES as an energy hub.
where T is the number of dispatching periods per day. In this paper, the dispatching time window of the IES is one day, and ∆t is the length of the dispatching period. C i (t) is the operation cost of the i-th kind of energy network at time t, which is determined by the operation plan made by the lower decision-maker. q j (t) is the marginal cost per unit power of the j-th kind of energy at time t, as shown in Equation (5).
where p grid (t) is the time-of-use (TOU) price of the external network at time t. C NG is the unit price of natural gas. Q LHV is the low calorific value of natural gas. P he (t) is the heating power of the micro gas turbine at time t. η he (t) is the heating efficiency of micro gas turbine at time t. In this paper, micro gas turbines are considered as heat network equipment. As a heat (cooling) source, micro gas turbines mainly take on a gas-to-heat (cooling) role, while utilizing waste heat to generate electricity. Micro gas turbines operate in a heat-to-electricity mode. q j (t) is also affected by the decision-making of the lower model. In the energy optimal dispatching model of the IES based on uncertain bilevel programming, random chance constraints are used for the objective function of the upper model, because the smaller the system operation cost, the better. Therefore, the objective function of the upper model can be obtained in the form of min-min, as shown in Equation (6): The upper objective function needs to meet the probabilistic constraints shown in Equation (7): where α up is the confidence of the upper objective function, and f up is the optimistic value of the objective function when the confidence of the upper model is α up . The upper model needs to meet the constraints as shown in Equation (8).
Energies 2020, 13, 477 The allocation factor balance and the power balance of each energy subnet are included in this equation.

The Lower Model
The control variables of the lower decision-maker are the specific optimal operation plans of the power electronic subnet, the thermal energy subnet, and the gas energy subnet. The plan is based on the transformation equation of the energy hub formulated by the upper decision-maker, and its objective is to minimize the operation cost of the lower decision-maker.

Objective Function
The objective function of the lower model is the comprehensive operation cost of the energy subnet, as shown in Equation (9).
where C e , C h , and C g are the daily comprehensive operation cost of the electric power subnet, the thermal energy subnet, and the gas energy subnet, respectively. The equation is where ∆t is the dispatching period. C e is the comprehensive operation cost of the electric power subnet. P e (t) is the purchasing and sale power from the external power grid at time t. P e,h (t) and P e,g (t) are the power of the electric power subnet converted into thermal energy subnet and gas energy subnet, respectively. F is the government's subsidized price for distributed photovoltaic power generation. P PV (t) is the power of the photovoltaic output at time t. P WT (t) is the power of the wind power output. k WT and k PV are the operation and maintenance cost coefficients of wind power and photovoltaic power, respectively. P SB (t) is the charge and discharge power of energy storage. u ch (t) and u dis (t) are the charge and discharge states of energy storage at time t: when it is charging at time t, u ch (t) = 1, u dis (t) = 0; when discharging, u ch (t) = 0, u dis (t) = 1. k ch is the cost coefficient of energy storage life loss during charging. k dis is the cost coefficient of energy storage life loss during discharging. The equations of k ch and k dis are shown in Equations (11) and (12), in which λ is the interruption compensation coefficient, and P cut (t) is load interruption capacity at time t. (12) where N(x) is the maximum number of cycles, which is determined by the depth x of charge and discharge. These equations are in [25]. C init is the initial fixed investment cost of energy storage. E SB,start is the initial charged states of the charge and discharge processes. E SB,end is the terminal charged states of the charge and discharge processes. E min SB is the minimum discharging depth (taken as 0.1 times the energy storage capacity). E max SB is the maximum allowed electric quantity of energy storage (taken as 0.9 times the energy storage capacity). c ch and c dis are the charging and discharging influencing factors, respectively.
Energies 2020, 13, 477 where C h is the comprehensive operation cost of the thermal energy subnet. P MT (t) is the CCHP Micro-turbine(MT) heating power at time t, which is also the energy input into the dual-effect absorption unit. η MT (t) is the CCHP MT heating efficiency at time t. θ is the classification of pollutants (the total is N). λ θ MT is the emission coefficient of the θ-th kind of pollutant of MT. c θ is the unit emission control cost of the θ-th kind of pollutant.
where C g is the comprehensive operation cost of the gas energy subnet. λ θ FC is the intake volume of the gas source point at time t. λ θ FC is the emission coefficient of the θ-th kind of pollutant of FC. P FC (t) is the generation power of the fuel cell at time t.
The objective function of the lower model is in Equation (15).
The probabilistic constraints that the lower objective function needs to satisfy are shown in Equation (16): Pr where α down is the confidence of the lower objective function. f down is the optimistic value of the objective function when the confidence of the lower model is α down .

Constraints (1) Electric Power Subnet Operation Constraints
The constraints that the electric power subnet decision-maker need to meet include the power balance, as shown in Equations (17) and (18), the energy conversion constraints, as shown in Equation (19), and the basic operation constraints of the electric power subnet, as shown in Equation (20). (18) where ∆P is the unbalanced power deviation of the electric power subnet. δ is the maximum operation deviation of the power balance, which is considered as 25 kW. β down is the confidence that power balance is satisfied, usually considered as 0.9. The chance constraint formula is used for the power balance in the lower model. In Equation (19), P e,h (t) and P e,g (t) are, respectively, the electric load power for heating and the power-to-gas power at time t. P h,e (t) and P g,e (t) are, respectively, the generation power of the CCHP unit and the fuel cell at time t. The equation of P loss (t) is in [27].
where C e,h (t), C e,g (t), and C e,t (t) are, respectively, the power-to-heat and power-to-gas energy coupling factors, which are determined by the upper decision-maker at time t.
Energies 2020, 13, 477 10 of 24 where S SB (t) and S SB (t + 1) are, respectively, the residual capacity at the end of time t and time t + 1. η dis is the discharging power of storage battery(SB). S min SB and S max SB are, respectively, the minimum and maximum residual capacity of energy storage. D SB is the self-discharging coefficient of energy storage. Q SB is the capacity of energy storage. P min SB and P max SB are, respectively, the minimum and maximum of energy storage output, respectively taken as −350 and 350 kW. P min grid and P max grid are, respectively, the minimum and maximum power for external exchange, respectively, taken as −500 and 500 kW. P max cut is the interruptible load capacity signed by the microgrid and users, considered to be 100 in the model.
(2) Thermal Energy Subnet Constraints Thermal energy subnet constraints include the thermal power balance, the cooling energy power balance, the energy input constraint, and the operation constraints of the energy storage device.
where C he and C co are, respectively, the coefficients of the heating and cooling of the dual-effect absorption unit at time t, respectively taken as 0.36 and 0.28. L he (t) and L co (t) are, respectively, the power of the heat load and the cooling load of the thermal energy subnet at time t. Q x (t) is the heat (cooling) stored in the energy storage device at time t. A value greater than zero indicates storage; otherwise, release is indicated. a is the residual heat coefficient of the micro gas turbine. X(t) and X(t − 1) are, respectively, the residual heat (cooling) of the energy storage device at time t and t − 1. λ x is the self-loss coefficient of residual heat (cooling) in the energy storage device, considered as 0.1. Q XU is the capacity of the energy storage device. When the energy storage device runs in the heating mode, the first constraint is satisfied; when it runs in the cooling mode, the second constraint is satisfied.
(3) Gas Energy Subnet The gas energy subnet constraints include gas balance, the operation constraints of the gas storage tank and the gas pipeline, and the energy conversion constraints.
where G s (t) is the gas emission from the gas storage tank at time t. G e,s (t) is the gas volume produced by power-to-gas conversion. L g (t) is the heating load of the natural gas. Q s (t) and Q s (t − 1) are, respectively, the residual gas of the gas storage tank at time t and t − 1, and Q min s and Q max s are the minimum and the maximum. G s (t) is the gas emission from the gas storage tank at time t, and G min s and G max s are the minimum and the maximum. G l (t) is the gas delivery capacity of the l-th pipeline at time t, and G min l and G max l are the minimum and the maximum. C e,g (t) is the energy coupling factor of the gas-to-power conversion at time t. C g,e (t) is the energy coupling factor of the power-to-gas conversion at time t, which is determined by the upper decision-maker.

Comprehensive Model
Combined with the upper model constructed in Section 2.2 and the lower model constructed in Section 2.3, a comprehensive energy optimal dispatching model of an IES based on uncertain bilevel programming was obtained, as shown in Equation (23).
where C is the energy coupling matrix of the upper decision variables. u e , u h , and u g are, respectively, the decision variables of the lower decision-maker, that is, the optimal operation plan of the electric power subnet, the thermal energy subnet, and the gas energy subnet. In this model, the energy interconnection matrix formulated by the upper decision-maker directly affects a part of their own objective function, namely the energy interconnection loss; on the other hand, it indirectly affects their own objective function by influencing the decision made by the lower decision-maker. The lower decision-maker's decisions require the precondition and foundation of the upper.

An Optimization Model Solving Algorithm
The energy optimization dispatching model of the IES proposed in this paper is a typical mixed integer nonlinear programming problem. Based on the features of this model, we used an improved firefly algorithm (FA) to solve it.

The Firefly Algorithm
In 2008, by exploring the mutual attraction and movement process among fireflies, Cambridge University scholars proposed a new group intelligence algorithm. It has been widely used in distribution network problems and generator set problems. However, in the FA, the location information of each individual firefly represents an optimization scheme within the space. Therefore, a traditional FA is insufficient to find an optimal solution and easily becomes a local optimal solution. On this basis, we constructed a new algorithm with which individuals can introduce a local-wide chaotic search, which can fully explore the search ability of each individual, thus improving the performance of the algorithm.
The main modeling idea of the FA is that individual fireflies with high brightness can guide the behavior of individuals with low brightness. The movement of the j firefly (weak brightness) to the i firefly (the strongest brightness) is mathematically expressed as shown in Equations (24) and (25). r 0 is the range of perception of each firefly.
X j (t+1) and X j (t), respectively, represent the position vector of firefly j at the t and t + 1 time periods; X j (t) represents the position vector of the brightest firefly in the sensing range of j in the t time period; β (0) is the attraction force when the distance between two fireflies is zero, and attraction decreases exponentially with distance when the distance is not zero; r is the distance between fireflies; γ is the light intensity coefficient, which is a constant, and represents the degree of inertia of each firefly's individual attraction; m is the number of fireflies. The parameter α in the formula is a random parameter between 0 and 1. u j is the random search volume adjusted by α. Therefore, the specific steps to improve the firefly algorithm are as follows: (1) Global Optimal Chaotic Search For the current globally optimal individual firefly, a chaotic search is performed within its sensing range. The search process is as shown in Equations (26) and (27). Logistic mapping is used as the chaotic motion formula. Firstly, each dimension position variable corresponding to the current firefly is mapped to [0,1]. Each iteration and the corresponding fitness of the firefly is recorded. When the iteration is completed, the firefly position variable with the highest brightness in the process is selected as the result of the chaotic search, which is finally mapped back to the original range of position variables.
In the formula, χ u+1 and χ u are the values of the chaotic variables at the u + 1 and u iterations, respectively; χ 0 represents the initial values of the chaotic variables, and there are fixed points including 0.25, 0.5, and 0.75 in the equation, so the initial values of the chaotic variables should be avoided.
(2) Adaptive Light Intensity Coefficient In the basic firefly algorithm, the light intensity coefficient is used as a constant value to characterize the degree of inertia of the firefly's individual attraction to the rest of the individual. In the optimization process, if a self-adaptive adjustment process can be introduced to the coefficient, it can give the current superior individual greater appeal. In order to further improve the efficiency of the firefly algorithm, an adaptive strategy is introduced for the light intensity coefficient in the speed update formula as shown in Equation (28) to form an adaptive light intensity coefficient firefly algorithm.
, f s,i ≤ f avg,s γ max , f s,i > f avg,s (28) γ s,i is the inertia weight of firefly i for iteration s, γ max is the inertia weight's upper limit, γ min is the inertia weight's lower limit. f s,i is the fitness value i for iteration s, and f avg,s is the average of all firefly fitness values at s. f min is the minimum fitness and f max the maximum fitness for all current fireflies.

The Solving Process
In this paper, the improved firefly algorithm is adopted to construct the energy optimization scheduling model of the established comprehensive energy system. The specific calculation is divided into two parts, namely creating the upper model and creating the lower model. The specific process is as follows: The process of the upper model is as follows: (1) Set the firefly population size and population parameters of the upper model, including the firefly number, the initial light intensity coefficient, the light intensity coefficient limit, the maximum attraction, the algorithm iteration number, and the chaos search algebra. corresponding to each firefly according to the cost. The lower the operating cost of the system, the higher the brightness. (4) For each individual firefly, select the individual firefly with the highest brightness in its perception range, calculate the distance and attraction between them, and update the position of the individual firefly. If the individual firefly is brightest within the range of perception, a chaotic search is conducted. (5) Determine whether the solution result of the algorithm converges to the set accuracy, indicating that the comprehensive cost index under the multi-energy complementary coordination plan corresponding to all fireflies is optimal, or whether the number of iterations has reached the set number of iterations. If "yes," then the current energy conversion matrix is output and the algorithm ends. Otherwise, return to Step (3).
The process of the lower model is as follows: (1) Set the firefly population size and population parameters, including the firefly number, the initial light intensity coefficient, the light intensity coefficient limit, the maximum attraction, the algorithm iteration number, and the chaos search algebra. At the same time, input the multi-energy complementary coordination plan provided by the upper model. (2) Within the scope of control variables, use the operation plan of each energy subnet equipment under the multi-energy complementary coordination plan as the location to randomly initialize the location of each firefly. Among them, the position information of each firefly corresponds to a device output plan. (3) Calculate the operating cost of the system corresponding to each firefly according to the objective function and constraint conditions of the underlying model, and calculate the brightness corresponding to each firefly according to the cost. The lower the operating cost of the system, the higher the brightness. (4) For each individual firefly, select the individual firefly with the highest brightness in its perception range, calculate the distance and attraction between them, and update the position of the individual firefly. If the individual firefly within the range of perception itself is the highest brightness, the chaos search has been conducted. (5) Determine whether the solution result of the algorithm converges to the set accuracy, that is the operating cost index of the equipment output plan corresponding to all fireflies has reached the optimal level, or whether the number of iterations has reached the set number of iterations. If so, the optimal cost will be fed back to the upper level solution model, and the algorithm will be finished. Otherwise, return to Step (3).
Based on this, the solving process of the energy optimal dispatching model for a comprehensive energy system can be obtained, using this improved firefly algorithm to perform the uncertain bilevel programming, as shown in Figure 3. so, the optimal cost will be fed back to the upper level solution model, and the algorithm will be finished. Otherwise, return to Step (3).
Based on this, the solving process of the energy optimal dispatching model for a comprehensive energy system can be obtained, using this improved firefly algorithm to perform the uncertain bilevel programming, as shown in Figure 3.

Simulation Example
In order to verify the validity and feasibility of the model, a pilot base in Hebei Province was selected as a research object for analysis of an IES example. The IES mainly consisted of a wind turbine, a photovoltaic generator, a cogeneration unit, an energy storage device, a methane reflection unit, a dual-effect absorption unit, and a fuel cell, and the system had obvious complementary characteristics. Its topological structure is shown in Figure 1.

Basic Data
The energy optimal dispatching model of an IES established in this paper takes into account the uncertainty of wind power generation, photovoltaic power generation, and load level, and the power curve data are the benchmark data. Simultaneously, the optimal dispatching time span is 1 day, divided into 24 runtime periods of 1 h. The system wind power output, photovoltaic power output,

Simulation Example
In order to verify the validity and feasibility of the model, a pilot base in Hebei Province was selected as a research object for analysis of an IES example. The IES mainly consisted of a wind turbine, a photovoltaic generator, a cogeneration unit, an energy storage device, a methane reflection unit, a dual-effect absorption unit, and a fuel cell, and the system had obvious complementary characteristics. Its topological structure is shown in Figure 1.

Basic Data
The energy optimal dispatching model of an IES established in this paper takes into account the uncertainty of wind power generation, photovoltaic power generation, and load level, and the power curve data are the benchmark data. Simultaneously, the optimal dispatching time span is 1 day, divided into 24 runtime periods of 1 h. The system wind power output, photovoltaic power output, power load, heat load, and cooling load curves within typical days in winter and summer are as shown in Figures 4 and 5. The above energy power levels follow the uncertain distribution; that is the wind power output adopts the Weill distribution [28], the photovoltaic output uses a distribution [29], and the load level uses a positive distribution [30]. power load, heat load, and cooling load curves within typical days in winter and summer are as shown in Figures 4 and 5. The above energy power levels follow the uncertain distribution; that is the wind power output adopts the Weill distribution [28], the photovoltaic output uses a distribution [29], and the load level uses a positive distribution [30].  Performance parameters of each piece of grid-connected equipment in the IES is shown in Table  2. Environmental conversion costs for cooling-heat-electric co-generation micro gas turbines and fuel cells are shown in Table 3. Whether purchasing or selling electricity to the outside network, the IES follows the TOU price level of the outside network, as shown in Table 4. the wind power output adopts the Weill distribution [28], the photovoltaic output uses a distribution [29], and the load level uses a positive distribution [30].  Performance parameters of each piece of grid-connected equipment in the IES is shown in Table  2. Environmental conversion costs for cooling-heat-electric co-generation micro gas turbines and fuel cells are shown in Table 3. Whether purchasing or selling electricity to the outside network, the IES follows the TOU price level of the outside network, as shown in Table 4. Performance parameters of each piece of grid-connected equipment in the IES is shown in Table 2. Environmental conversion costs for cooling-heat-electric co-generation micro gas turbines and fuel cells are shown in Table 3. Whether purchasing or selling electricity to the outside network, the IES follows the TOU price level of the outside network, as shown in Table 4.   In the process of IES energy optical dispatching, the gas price C NG is 2.06 ¥/m 3 ; the subsidy price of distributed photovoltaic generation F is 0.42 ¥/kWh; load interruption compensation coefficient β is 3; SB discharging efficiency η dis and SB charging efficiency η ch are both 0.88; the self-discharge coefficient of energy storage D SB is 0.1; the capacity of energy storage Q SB is 1000 kWh; the self-loss coefficient of the remaining heat (cooling) of energy storage device λ x is 0.1; the capacity of gas storage tank Q max s is 400 m 3 . In addition, in the improved firefly algorithm, it was assumed that the quantity of the initial state is 100, the initial firefly number is 20, the maximum number of iterations is 200 generations, and D is 100.

Scenario Setting
For the simulation analysis of the model proposed in this paper, two main scenarios were set up for discussion. One was the IES on a typical winter day taking into account the electrical and heat loads of the system, where the dual-effect absorption unit is in the heating state; the other was the IES on a typical summer day taking into account the electrical and cooling loads, where the dual-effect absorption unit is in the refrigeration state. From the given confidence level, various confidence levels, and different algorithms, the system operation schemes of typical days in winter and summer were further analyzed.

Simulation Results in Different Scenarios
In the energy optimization dispatching model, the confidence degree of the objective function was 0.9 in the upper model and the lower model, and the confidence degree was 0.95 in the constraint condition. The comparison discussion is as follows.

Results for a Typical Winter Day
For a typical winter day's operation, the IES energy optimization dispatching model established in this paper obtained an optimal operation scheme as shown in Figure 6.
As can be seen in Figure 6, the operation scheme of the IES fully invokes various power generation equipment and energy supply equipment, which satisfies the electric load and heat load demand of the system. From the 1st to the 5th period of the day, the external network TOU power price level is low, the energy storage is in the charging state, and the methane reaction unit is basically operating in a power-to-gas state. However, due to the high wind power output in the system and the low level of the electric load, the IES does not purchase electricity on a large scale from the external network. In the energy storage device charging and methane units, the electricity for power to gas (P2G) is basically from the distributed wind power output. From the 6th to the 8th period, the first peak of electric load happens, and the distributed new energy output is unable to meet the needs of system operation. As a result, the fuel cell is turned on, and the system purchases electricity from the external network. However, since the heat load level is still low, the economics of the operation of the CCHP unit is not superior; thus, the MT output is not large. From the 9th to the 13th period, the system is still at the peak of electric load, but also at the peak of heat load. The IES makes full use of the energy storage device discharge, and the CCHP unit provides the heat and power with the fuel cell output, while fully utilizing the methane reaction unit for gas-to-power activities. For the heat load, large-scale heat purchased from the thermal energy external network is also carried out to meet the system operation requirements. Since the external network TOU power price is relatively high, and the internal power supply cost of the IES is lower, the system sells electricity to the external network to obtain revenue. In the 14th period, the electric load enters a valley period. The methane reaction unit is transferred to the power-to-gas mode due to the relatively high operating cost, and the natural gas produced through P2G is input into the gas storage tank. From the 15th to the 23rd period, the system enters the second electric load peak, while the heat load level is still high. The operating conditions of the system are similar to those from the 9th to the 13th periods, but the difference is that, during the 20th period, the electric load reaches a peak of one day. Therefore, the system needs to buy electricity from the external network and even carry out a more expensive interruptible load to ensure the reliability of system operation.
Energies 2020, 13, x FOR PEER REVIEW 18 of 26 levels, and different algorithms, the system operation schemes of typical days in winter and summer were further analyzed.

Simulation Results in Different Scenarios
In the energy optimization dispatching model, the confidence degree of the objective function was 0.9 in the upper model and the lower model, and the confidence degree was 0.95 in the constraint condition. The comparison discussion is as follows.

Results for a Typical Winter Day
For a typical winter day's operation, the IES energy optimization dispatching model established in this paper obtained an optimal operation scheme as shown in Figure 6. As can be seen in Figure 6, the operation scheme of the IES fully invokes various power generation equipment and energy supply equipment, which satisfies the electric load and heat load demand of the system. From the 1st to the 5th period of the day, the external network TOU power price level is low, the energy storage is in the charging state, and the methane reaction unit is basically operating in a power-to-gas state. However, due to the high wind power output in the system and the low level of the electric load, the IES does not purchase electricity on a large scale from the external network. In the energy storage device charging and methane units, the electricity for power to gas (P2G) is basically from the distributed wind power output. From the 6th to the 8th period, the first peak of electric load happens, and the distributed new energy output is unable to meet the needs of system operation. As a result, the fuel cell is turned on, and the system purchases electricity from the external network. However, since the heat load level is still low, the economics of the operation of the CCHP unit is not superior; thus, the MT output is not large. From the 9th to the 13th period, the system is still at the peak of electric load, but also at the peak of heat load. The IES makes full use of the energy storage device discharge, and the CCHP unit provides the heat and power with the fuel cell output, while fully utilizing the methane reaction unit for gas-to-power activities. For the heat load, large-scale heat purchased from the thermal energy external network is also carried out to meet the system operation requirements. Since the external network TOU power price is relatively high, and the internal power supply cost of the IES is lower, the system sells electricity to the external network to obtain revenue. In the 14th period, the electric load enters a valley period. The methane reaction unit is transferred to the power-to-gas mode due to the relatively high operating cost, and the natural gas produced through P2G is input into the gas storage tank. From the 15th to the 23rd period, the system enters the second electric load peak, while the heat load level is still high. The operating conditions of the system are similar to those from the 9th to the 13th periods, but the difference is that, during the 20th period, the electric load reaches a peak of one day. Therefore, the According to the above analysis, on a typical winter day, the IES adopts multiple power supply and heating modes to perform a multi-energy complementary operation, which satisfies the user's demand and improves the social and economic benefits of IES operation. The comprehensive operating cost of the system is ¥7890. 25.
In fact, the operating cost curve of the IES can be obtained for each period of a typical winter day, as shown in Figure 7.
Energies 2020, 13, x FOR PEER REVIEW 19 of 26 system needs to buy electricity from the external network and even carry out a more expensive interruptible load to ensure the reliability of system operation. According to the above analysis, on a typical winter day, the IES adopts multiple power supply and heating modes to perform a multi-energy complementary operation, which satisfies the user's demand and improves the social and economic benefits of IES operation. The comprehensive operating cost of the system is ¥7890. 25.
In fact, the operating cost curve of the IES can be obtained for each period of a typical winter day, as shown in Figure 7.

Results for a Typical Summer Day
For a typical summer day's operation, the IES energy optimization dispatching model established in this paper obtained an optimal operation scheme as shown in Figure 8.

Results for a Typical Summer Day
For a typical summer day's operation, the IES energy optimization dispatching model established in this paper obtained an optimal operation scheme as shown in Figure 8. It can be seen in Figure 8 that the principle of the IES optimization operation scheme is basically the same as that for a typical winter day. In each peak and valley period of the system, the energy storage equipment of the power subnet, the storage device of the thermal energy subnet, and the methane reaction unit of the gas energy subnet and the gas storage tank all play a role in the peak shaving shift. Simultaneously, comparing the operation of a typical summer day and that of a typical winter day, the main differences lie in the wind power output, the level and distribution characteristics of the photovoltaic output, the load level, and the curve shape. Therefore, the minimum operating cost of a typical summer day is ¥8658.77 at a given confidence degree. It can be seen in Figure 8 that the principle of the IES optimization operation scheme is basically the same as that for a typical winter day. In each peak and valley period of the system, the energy storage equipment of the power subnet, the storage device of the thermal energy subnet, and the methane reaction unit of the gas energy subnet and the gas storage tank all play a role in the peak shaving shift. Simultaneously, comparing the operation of a typical summer day and that of a typical winter day, the main differences lie in the wind power output, the level and distribution characteristics of the photovoltaic output, the load level, and the curve shape. Therefore, the minimum operating cost of a typical summer day is ¥8658.77 at a given confidence degree.
In fact, the operating cost curve of the IES can be gained for each period of a typical summer day, as shown in Figure 9. As we can see in the figure, compared with the environment conversion cost, the IES has a larger proportion of expenditure in terms of economic cost. The operating cost curve of each period is basically positively correlated with the comprehensive load level of the system. In fact, the operating cost curve of the IES can be gained for each period of a typical summer day, as shown in Figure 9. As we can see in the figure, compared with the environment conversion cost, the IES has a larger proportion of expenditure in terms of economic cost. The operating cost curve of each period is basically positively correlated with the comprehensive load level of the system.

Result Comparison
We can compare the operation results between the upper decision-maker and the lower, subdecision-maker under the two operation modes of typical winter and summer days, as shown in

Result Comparison
We can compare the operation results between the upper decision-maker and the lower, sub-decision-maker under the two operation modes of typical winter and summer days, as shown in Table 5. Obviously, subnet operation shows that the operating cost of the power subnet is the largest, followed by that of the thermal subnet and subsequently that of the gas energy subnet. The cost of the operation of the IES in each period of a day can be seen visually. As can be seen in Table 5, under a given confidence level, the operating cost of the IES on a typical winter day is lower than that on a typical summer day. This is mainly due to the adoption of a multi-energy complementary and coordinated IES. In view of the heat load of the typical winter day, the IES can fully formulate a coordination plan among the electric power subnet, the thermal energy subnet, and the gas energy subnet, which further reduces the heating cost of the system. In addition, in order to show the operational advantages of the IES, considering the same load level, the operating cost of the traditional microgrid with only one form of energy supply is ¥8863.07 in winter and ¥10,774.88 in summer. Compared with the operation mode proposed in this paper, the operating cost of the IES in winter and summer is obviously reduced by 12.33% and 33.7%, respectively.

Multiple Algorithm Comparison
In order to verify the effectiveness of the improved firefly algorithm proposed in this paper, the basic particle swarm optimization algorithm, the basic firefly algorithm, and the improved firefly algorithm in this paper were used to solve the same IES optimization operation model, and the performance indexes under various circumstances were compared, as shown in Table 6. Compared with the particle swarm optimization and the basic firefly algorithm, the algorithm proposed in this paper has obvious advantages not only in running cost but also in solving time.

Sensitivity Analysis of Multi-Confidence
In order to compare the operating results under different confidence values, we selected other confidence values to obtain the minimum operating cost of the IES under the corresponding conditions, as described below. (1) For a typical winter day's operation, we selected the optimal operating cost distribution under different objective function confidence values and different constraint conditions, as shown in Figure 10.
proposed in this paper has obvious advantages not only in running cost but also in solving time.

Sensitivity Analysis of Multi-Confidence
In order to compare the operating results under different confidence values, we selected other confidence values to obtain the minimum operating cost of the IES under the corresponding conditions, as described below.
(1) For a typical winter day's operation, we selected the optimal operating cost distribution under different objective function confidence values and different constraint conditions, as shown in Figure 10. (2) For a typical summer day, we selected the optimal operating cost distribution under different objective function confidence values and different constraint conditions, as shown in Figure 11. (2) For a typical summer day, we selected the optimal operating cost distribution under different objective function confidence values and different constraint conditions, as shown in Figure 11. As can be seen in Figures 10 and 11, on typical winter and summer days, as the confidence of the objective function and the confidence of the constraint increase, the optimal operating cost of the system continuously increases. The results show that, with the increase in the confidence of the objective function, the confidence degree of the minimum value of the objective function is greater, and the operating plan developed by the system is more conservative, but the operating cost of the system increases. When the confidence degree of the constraint condition increases, the confidence of the system with respect to the power balance constraint requirement will also increase, which increases the operating cost of the system.
Based on this, in the development of an energy optimization dispatching plan for an IES, main decisions depend on the decision-maker's own risk preferences and realistic requirements. With increasing confidence, reliability requirements for the establishment of the objective functions and constraints of IESs will become increasingly stringent, planned operating costs of the system will increase, and reliability should be compensated by sacrificing economy. If the decision-maker is risktaking, or the actual operation requires less reliability, they will tend toward a lower objective function confidence and fewer constraint conditions to achieve a lower cost; if the decision-maker is risk-averse, or the actual operation requires higher reliability, they will tend toward a higher target function confidence and more constraint conditions, and they will compensate for the reliability with a higher planned operating cost. Figure 11. The relationship between the planned operating cost and the confidence degree of the IES for a typical summer day. like to express our gratitude to them for their help and guidance. The optimistic value of the objective function when the confidence of the lower model is α down c i j (t) Energy coupling factor in t period ∆P Power subnet unbalanced power deviation P i (t) Input power of the i-th energy grid in t period δ Power balance maximum running deviation L i (t) Load power of the i-th energy grid in t period β down Confidence in meeting power balance conditions C(t)

Conflicts of
Energy coupling matrix of IES in t period P e,h (t) Heating load power in t period

P(t)
Energy input vector in t period P e,g (t) Power to gas conversion power in t period L(t) Energy load vector in t period P h,e (t) Cogeneration unit power generation in t period v ij (t) Allocation factor c in t period P g,e (t) Fuel cell power generation in t period η i j (t) Efficiency factor in t period P loss (t) Loss of power in t period T Number of scheduling days per day C e,h (t) Energy coupling factor of power to heat in t period ∆t Length of scheduling time C e,g (t) Energy coupling factor of power to gas in t period C i (t) Operating cost of the i-th energy network in t period S SB (t) SB remaining power at the end of t period q j (t) The marginal cost per unit power of the j-th energy source in t period S SB (t + 1) SB remaining power at the end of t + 1 period p grid (t) External network TOU price level in t period η dis SB discharge efficiency C NG Natural gas unit price η ch SB charging efficiency Exchange power maximum with external network P e (t) The power of power subnetwork to buy and sell power to network in t period P max cut Interruptible load capacity signed by the microgrid and the user P e,h (t) the power of power subnetwork convert to heat subnetwork in t period C he Heating coefficient of double-effect absorption unit P e,g (t) The power of power subnetwork convert to gas subnetwork in t period C co Refrigeration coefficient of double-effect absorption unit F Government subsidized price for distributed photovoltaic power generation L he (t) Heat load power of thermal energy subnetwork in t period P PV (t) Photovoltaic output power in t period L co (t) Cooling load power of thermal energy subnetwork in t period P WT (t) Wind power output in t period Q x (t) Heat stored in the energy storage device (cooling capacity) in t period k WT Wind power operation and maintenance cost coefficient a Residual heat coefficient of micro gas turbine k PV Photovoltaic operation and maintenance cost coefficient Residual heat (cooling capacity) of the energy storage device during the t period P SB (t) Charge and discharge power of stored energy in t period X(t − 1) Residual heat (cooling capacity) of the energy storage device during the t − 1 period u ch (t) The state of being charged with stored energy in t period λ x Self-loss coefficient of residual heat (cooling capacity) of energy storage device u dis (t) The discharge state of stored energy in t periodQ XU Capacity of energy storage device k ch Energy storage life loss cost factor during charging G s (t) The amount of gas released from the gas storage tank in t period k dis Energy storage life loss cost coefficient during discharge G e,s (t) The amount of gas produced by power to gas in t period λ Interruption compensation coefficient L g (t) Natural gas heating load Energies 2020, 13, 477 23 of 24 P cut (t) Load interruption capacity in t period Q s (t) Gas tank gas remaining amount in t period N(x) Maximum number of cycles Q s (t − 1) Gas tank gas remaining amount in t − 1 period C init Energy storage initial fixed investment cost Q min