Operation Optimization of Integrated Energy System under a Renewable Energy Dominated Future Scene Considering Both Independence and Beneﬁt: A Review

: An integrated energy system interconnects multiple energies and presents a potential for economics improvement and energy sustainability, which has attracted extensive attention. However, due to the obvious volatility of energy demands, most existing integrated energy systems cannot operate in a totally self-sufﬁcient way but interact with the upper grid frequently. With the increasingly urgent demand for energy saving and emissions reduction, renewable resources have occupied a larger and larger proportion in energy system, and at last they may be dominant in the future. Unlike conventional fossil fuel generation, the renewable resources are less controllable and ﬂexible. To ease the pressure and guarantee the upper grid security, a more independent integrated energy system is required. Driven by that, this paper ﬁrstly reviews the optimal strategies considering both independence and beneﬁt from perspectives of individual efforts and union efforts. Firstly, the general optimization process is summarized in terms of energy ﬂows modelling and optimization methods to coordinate supply–demand side and realize beneﬁt maximization. Based on that, handling with uncertainty of high-ratio renewable energy is reviewed from uncertainty modeling methods and multi-stage operation strategy perspectives to make the strategy accurate and reduce the adverse effects on the upper grid. Then, the hybrid timescale characteristics of different energy ﬂows are explored to enhance operation ﬂexibility of integrated energy systems. At last, the coordination among different participants is reviewed to reduce the whole adverse effect as a union. Remarks are conducted in the end of each part and further concluded in the ﬁnal part. Overall, this study summarizes the research directions in operation optimization of integrated energy systems to cater for a renewable energy dominated scene to inspire the latter research.


Introduction
At present, the energy structure worldwide is dominated by primary energy, which raises a dilemma between the increasing demand and sustainable supply of limited energy resources [1]. To face that, vigorous renewable energies (RE) have been explored as alternative resources of fossil fuels [2], for instance, wind and solar generation. Since 2010, the installed capacity of renewable energy has increased at an average annual rate of 8%. In 2015, the United Nations adopted two historically significant agreements including the Paris Agreement [3] and the 2030 Agenda for Sustainable Development [4] to reduce global pollutant emissions. In the meantime, a novel energy structure named integrated energy system (IES) is proposed to integrate and coordinate electricity, gas, heat and other energies in a comprehensive manner to enhance energy utilization efficiency [5]. In an IES, multiple energies are interconnected and converted by various coupling components to satisfy terminal demands. With reasonable operation strategy, the correlation and independence of a single IES and relieve the upper energy system burden. On the basis of the above discussion on effort directions, this paper aims to summarize the existing work of each aspect in IES operation optimization and give our insights to inspire future researches. To have a better illustration, the efforts directions in IES independence promotion are displayed in Figure 1. IES is a promising energy utilization form to realize energy saving and economic improvement. To cater for an RE dominated energy structure in future, the independence improvement for IES has great significance. To our best knowledge, there is no related review work conducted from this new perspective aiming at IES particularity and facing the future scene. To fill this gap, this paper firstly discusses the effort directions in making optimal strategy with consideration to both independence and benefit, and then reviews the latest work of each aspect separately. The contributions of this paper are as follows: (1) This paper firstly analyzes the causes for dependence and discusses the effort directions to improve the IES operation independence from the perspectives of individual efforts and union cooperation. On this basis, this paper aims to summarize the existing work of each aspect in IES operation optimization and give our insights. To our best knowledge, this is the first time the IES operation optimization considering independence improvement to cater for an RE dominated energy structure in future has been focused on. (2) To reduce the adverse impact from the RE on independence, the uncertainty handlings are reviewed from the perspectives of uncertainty modelling methods and multi-stage optimization strategy to improve the accuracy of optimization results and reduce interactions with the upper grid. (3) To reduce the adverse impact from the strong couplings among multiple energy flows in IES, this paper reviews the exploration of the potential buffer in heat and gas pipelines. To our best knowledge, this is the first time the considerations of hybrid timescale characteristics of multiple energy flows in IES operation optimization have been comprehensively reviewed. (4) To reduce the adverse impact from the volatile multi-energy demands, this paper reviews the coordination among different participants to smooth load profile with price signals. The concept, mechanisms, and applications of game theory are introduced under the scenes of single IES and multiple IESs in a future energy network. To our best knowledge, this is the first time a review on the corporation of different interest participants in these two scenes has been conducted. The rest of this paper is structured as follows: Section 2 reviews the unified IES optimal operation process to realize individual benefit maximization. Based on that, more considerations can be added. Sections 3 and 4 introduce the directions of individual efforts in independence promotion. Section 3 concludes the handling methods facing RE uncertainties from two perspectives (i.e., uncertainty modelling and multi-stage dispatching) to reduce the adverse impact of RE on operation optimization results. Section 4 introduces the multiple time-scale characteristics of hybrid energy flows and separately reviews their depiction approaches in operation optimization to explore the potential buffer and increase system operational flexibility. Section 5 introduces the directions of union efforts and reviews the coordination among different participants in resources and interest to adjust the relation between supply and demand. The remarks and future directions are respectively discussed at the end of each section, and the study is moreover concluded with final remarks in Section 6.

General Formulation of IES Optimal Operation
The IES operation optimization problem is formulated based on mathematical modeling and optimization formulation. Firstly, the models of energy carriers and energy flows are built to represent device physics characteristics, energy transformation, and coupling relationship. Then, the objectives are developed as a baseline estimation in accord with requirements and realized with some optimization techniques. To make the results feasible and practical, constraints such as operation limits need to be considered. This part aims to give an intuitive show of the basic process of IES operation optimization considering individual benefit. Based on that, more considerations can be incorporated to promote independence. The device modelling methods have already been summarized in Ref. [12], and the other parts are summarized as follows.

Hybrid Energy Flows Modelling of IES
An accurate mathematical model is essential to depict and analyze the hybrid power flows among different devices. In present, there are two mainstream methods in IES energy flows modelling, i.e., bus based model and energy hub (EH) based model. In the operation optimization of most basic IES, such as a local IES, the different characteristics of energy flows are neglected, and the supply-demand relation is regarded as real-time balanced.

Bus Based Model
Bus based model [39][40][41] is the most commonly utilized hybrid energy flow modelling method. It can explicitly reflect the entire process of energy flows and energy coupling relationships. Its basic idea is to classify the bus according to energy transfer mediums (for instance, electricity bus, heat bus, cooling bus, gas bus) and then list energy flow balance equation in each bus node. Specifically, the main steps are the following: Step 1: analyze the energy flow relationship of IES in detail.
Step 2: define variables and list energy flow expressions of each component.
Step 3: write energy flow balance equation in each bus node in accord with the principle that the sum energy flows should be zero. Figure 2 shows a typical IES framework based on bus structure. In this system, electricity demand is satisfied by gas turbine unit and wind turbine jointly, while power grid is as a supplement. The waste heat from gas turbine is recovered and supply for heat load as well as absorption chiller driven. The output of absorption chiller is used to satisfy users' cool demand. Gas turbine, waste heat boiler, and absorption chiller realize energy cascade utilization in IES and improve primary energy efficiency. Additionally, gas boiler and electric cooler are set as auxiliary devices in case the heat or cooling load is greater than the supply of trigeneration unit. Three kinds of storage devices including electricity storage, cooler storage, and heat storage are used for peak load shifting and reducing excess energy. In the above IES, P WT , P in represent wind turbine output and input power from the power grid. G GT , G GB denote input natural gas of gas turbine and gas boiler, respectively. η e and η h represent electrical and thermal efficiency accordingly, while η GB represents gas boiler efficiency. COP ec , COP ab are the performance coefficient of electric cooler and absorption chiller. H ab represents the consumed heat power by absorption chiller. P ES , P TS , P CS , respectively, denote power exchange with energy bus of these three kinds of storage devices, where the subscripts c and d indicate energy absorption mode and energy release mode, respectively. L e , L h , L c denote power demand of different load. As aforementioned, there are three buses including electricity bus, heat bus, and cooling bus. The energy balance equation in each bus is listed in Equation (1).    electricity : P WT + P in + η e G GT + P ES,c = L e + P ec + P ES,d heat : η h G GT + η GB G GB + P TS,c = L h + H ab + P TS,d cooling : COP ec P ec + COP ab H ab + P CS,c = L c + P CS,d

Energy Hub Based Model
The concept of energy hub (EH) was put forward by Geidl et al. in 2007, and it was defined as a unit that provided the features of input, output, conversion, and storage of multiple energy carriers [42]. An EH consumes at the input ports and provides certain energy to users at the output ports through internal various energy carriers. The proposal of EH provides a new means of IES modeling and analysis [43]. The EH modeling method emphasizes the output and input characteristics from the perspective of efficiency analysis, while the internal complex structure is regarded as a black box. As shown in Figure 3 is a universal model of EH, wherein the L i (i = 1, 2 . . . m) and P j (j = 1, 2 . . . n) denote the i-th output energy power flow path and the j-th input energy power flow path, respectively. A coupling matrix C is defined to reveal the conversion relation among output and input ports wherein c ij is the relation coefficient. Specifically, the energy flow modelling can be described as Equation (2).  The EH based modelling method is from the angle of entirety to depict external characteristics of IES, and the core point is to obtain the coupling matrix C. In general, there are four steps for it: Step 1: define input and output power variables, namely L vector and P vector.
Step 2: define dispatch factors ν at input connection point of energy conversion device.
Step 3: denote the output of energy conversion device with its respective function.
Step 4: list the power balance equation at the output connection point and then arrange it in standard matrix form.
An example is given in Figure 4, wherein ν GT represents dispatch factor of gas turbine and other variables are as the same as Figure 2. The corresponding energy flow model is established in Equation (3). The high integrity and compactness of EH-based modeling attract much attention, and the coupling matrix has been expanded to combine more components, including demand response [44], storage devices [45,46], plug-in electrical vehicles [47], and renewable resources [48].

Remarks on These Two Modelling Methods
The bus based modelling method is conducted from the angle of locality to separately model each device and energy system. It is a quasi-steady state model and is easy to describe the system topology. The advantage is that it can explicitly reflect the entire process of energy flows and energy coupling relationships. However, the modelling process is based on the detailed analysis of energy flows to build balance equations in each bus. For a simple system, it is easy to carry out and express clearly. As for a complicated IES with more devices, energy types, and more intricate coupling relationships, the conventional bus base modelling method will have low efficiency and may complex equation formulation, reduce intuitiveness, and aggravate solving difficulty. In addition, some critical information (for instance, the transmission characteristics and network losses) is simplified or ignored, Energies 2021, 14, 1103 8 of 36 and the results accuracy is influenced. Hence, it is more applicable for an IES with simple structure.
As for EH based energy flow modelling method, it is conducted from the angle of entirety and concentrates more on the relation between input and output instead of the specified internal interconnections. It has better scalability and extensibility, which facilitate the calculation process and are effective in system operation optimization. Due to these advantages, it is especially applicable for the modelling of increasing complicated IESs. The EH base modelling method has attracted many researchers' attention and has grown rapidly in recent years. However, the coupling matrix is not easy to obtain, especially with more elements such as storage devices and demand response. Additionally, the standardized modeling method contains dispatch factors and products of decision variables in coupling matrix, and the operation problem becomes nonlinear with a high computation burden. Hence, it proposes higher requirement for optimization techniques.

Operation Optimization of IES
As denoted in Equation (4), the IES optimal operation problem is essentially an optimization problem with variables, objectives, operation constraints, and optimization techniques. f i (x) denotes i-th objective function, x denotes n-dimensional decision variables, m denotes the number of objective functions, and g j (x)(j = 1, 2, . . . , p) and h k (x)(k = 1, 2, . . . , q) represent the j-th equality and k-th inequality constraint function, respectively. The decision variables include a number of operating devices, on-off signals, output of each generation device, storage state at different time slots, etc. More detailed description about optimization objectives, major constraints, and solving methods are concluded separately as follows.

Objectives
From the whole system perspective, multiple objectives are presented in the optimization model with different purposes. Aiming at different objectives, the optimization results provide theoretical guidance for the system scheduling and transactions. The optimization results change dramatically when the optimization objectives are different. Normally, the operation economy is the most concerning issue to have a better market competitiveness and benefit. The economy evaluation can be subdivided into cost reduction and benefit improvement, and they are represented with corresponding functions. For the former, it aims to minimize the total cost in the operation. In general, the main cost includes operation cost (such as the fuel cost, start/off cost and maintenance cost) and purchasing cost from external energy markets. To relieve the pollution, it is noted that the concern for the environmental impact is growing. The excess carbon cost and penalty cost of wind curtailment are also included for some scenes [49]. From the benefit perspective, it aims to increase the potential income through energy transactions. The cost-benefit evaluation can also be considered with the concept of net profit, which indicates the difference between the total benefit and the total cost. With the pursuit of these economic objectives, the cost and benefit balance can be achieved, and the optimal economic operation strategy is obtained. The common economic indicators include the total operation cost, profit, social welfare [50,51], and so on.
On the basis of the economy pursuit, the concern of energy savings is also introduced, sometimes faced with the energy crisis. The energy saving is evaluated from the perspectives of energy generation [52,53], energy consumption [54], and energy utilization efficiency [55,56]. Some typical energy saving indicators, including primary energy consumption (PEC), primary energy ratio (PER), and primary energy efficiency (PEE), are utilized to increase energy savings. Additionally, the environmental protection also raises attention to reduce adverse impact on the environment during operation, such as to minimize the NOx emissions [57], the carbon emissions [58,59], and the total emissions [60]. With all these considerations from different perspectives, the objective may be multiple. Common indicators are concluded as follows: Faced with different demands and preferences, the optimization objectives are various, and the optimal operation strategy varies greatly. The objectives may be multiple and sometimes mutually exclusive. For instance, the improvement of environmental impact will bring reduction of economic benefit. Thus, it necessitates an efficient coordination approach to promote the development of multi-objective optimization. Additionally, the existing indicators take the different grade energies equally, which is not reasonable from an energy utilization perspective. Some exergy analysis methods can express the difference from the second law of thermodynamics, yet are normally adopted in system thermodynamic analysis [61,62] and seldom utilized in the IES operation optimization process. In summary, a comprehensive and systematic evaluation system to better guide optimization directions is lacking.

Constraints
Each component in IES has its operation characteristics and technical constraints, which must be considered in dispatching to guarantee practicality and authenticity. When associated with other parts in IES, additional constraints are imposed to satisfy device physical conditions, and operators' or users' specific demand at a system level. The main constraints in a single IES are concluded and listed as below.

Optimization Techniques
The optimal IES operation has been studied extensively. Among them, the major solution methods can be divided into mathematical optimization methods and heuristic methods.
The conventional mathematical methods are designed with appropriate iteration formulas and get the approximate solution finally after multiple iterations. They usually start with a given initial point and calculate the next iteration point according to the descending information such as gradient. Thus, they have fast convergence speed and high searching efficiency. Many typical mathematical methods have been used in IES operation optimization problems, such as the Lagrangian Relaxation (LR) method [69], Generalized Reduced Gradient (GRG) method [70], Interior Point (IP) method [71] and Benders Decomposition (BD) method [72]. The GRG is easy to implement but has poor convergence, and hence it is not applicable for a complicated problem. The IP method performs well in many aspects, but is not applicable for all problems. Likewise, the LR method is also faced with the problem of computational complexity. Compared with others, the BD method can reduce the computational complexity and is more effective for a large-scale and complicated problem. However, it has slow convergence rate and low efficiency. In general, the mathematical optimization methods have a long history of developments and are relatively mature in theoretical research. Hence, they are easier to conduct and implement, even though the computing cost of these methods are expensive, especially when solving a complex problem. Additionally, they require more information, which is not always known beforehand. Hence, they are more applicable to a small-scale IES operation optimization problem with fewer devices and interconnection relationships.
Intelligent methods are a kind of experience-based approach, which have been rapidly developed and implemented on various optimization problems. They are inspired by behaviors, reactions, and communication mechanisms in nature and imitate the process to seek the optimal solution in a more effective group search way. The most extensively utilized intelligent methods include Genetic Algorithm (GA) [73], Particle Swarm Optimization (PSO) [74], Bacteria Foraging Algorithm (BFA) [75], Dance Bee Colony Algorithm [76], etc. The GA method has good diversity of solutions. Since it is easy to implement, it is extensively utilized. However, it is probable to lose the optimal solutions and needs long time for convergence. To solve more complex problems and improve calculation efficiency, many algorithms derive into many improved forms, for instance, Non-dominated Sorting Genetic Algorithm (NSGA) and NSGA-II [77]. Compared with GA, the PSO method is easier to find the global optimal solutions, yet it is not sure to be strictly convergent and is probable to be immersed in local optima in multi-modal problems. Some improved variants also appeared like multi-objective PSO(MOPSO) [78], and Adaptive Multi-Objective Particle Swarm Optimization (AMOPSO) [79]. The SA method is more probably to find the global optimal solution, but it is time consuming and may be constrained by initial conditions. Compared with classical mathematical methods, intelligent methods need neither extra initial points nor gradient information of objective functions, and additionally have dominant superiority in universality, parallel computing ability, and possibilities of global optimization solution obtained. Thus, they are more suitable and efficient for an optimization problem with little information/large-scale/high-dimension/multiple objectives. However, the conventional intelligent methods are easy to fall into local optimum and the convergence speed is expected to be improved. Detailed mechanism description can be found in Ref. [80], and the comparison of different optimization methods are concluded in Table 1. Poor adaptability to complex problems, sensitive to external environment change From the comparison of different methods, it is found that there is no one that performs the best in all aspects. Although there are already many studies and continuous developments, there still exist difficulties to be solved or performance to be improved. With the fast development of IES, it tends to combine more energies and devices, which will complex the operation optimization problem in both calculation burden and difficulty. Hence, the computation complexity remains to be reduced, and computation time shortened, and computation efficiency needs to be improved. Decomposing the original complicated problem into some small sub-problems or using some dimensionality reduction techniques is a promising direction, yet the implementation in mathematics is the most difficult part. Additionally, from the comparison in Table 1, it can be seen that each method is unique with distinct features. The calculation performance relies on the characteristics of optimization problems, i.e., objectives, variables types (continuous, discrete, or integer), and constraints forms (linear or nonlinear), as well as convexity of problem. At present, the selection of the optimization algorithm for a concrete operation optimization problem often relies on expert knowledge or production experiences and existing previous research. To assess the performance of application optimization results, a standard or benchmark is necessitated. Even though some indexes such as convergence, computation time, and uniformity are normally adopted, a comprehensive and systematic evaluation method is still lacking.

Handling Uncertainty of High Accommodation Renewable Energy in IES
To reduce pollutant emission and faced with the fossil energy crisis, the installed capacity of renewable energy is rapidly growing worldwide, and an RE-dominated system is proposed for the future. In this context, renewable resources will also account for an increasingly high ratio in IES. However, the power generation of renewable resources such as wind and solar is weather-dependent and with strong uncertainties. Normally, the output of RE is estimated based on prediction means and then incorporated into the IES operation optimization problem. Inevitably, errors exist between predicted and real value, which reduces the accuracy of results and causes adverse impact on the interaction with the upper grid. Compared with uncertainties from other issues, such as loads and device technical performance deviation, the RE has the greatest uncertainty, especially for a high RE proportion IES. Different from single power system, the uncertainty is easy to expand through strong coupling and interconnections among different energy flows and various energy conversion devices in IES. The uncertainty cascade will influence the energy supply stability and cause serious safety problems, especially for an IES with less controllable units. In addition, the IES will require more support from the upper grid to maintain normal operation. Therefore, the consideration of RE uncertainty is particularly vital in making IES operation strategy, and the effective handling method will be a direction of individual efforts to promote IES independence and benefit. Some measures have been taken to tackle with that, and it is categorized in this paper into three terms including (1) increasing prediction accuracy, (2) modelling RE uncertainty, and (3) step-wise eliminating the adverse influence by multi-stages optimization. As the first category is restricted to technical level, the last two are respectively reviewed in this paper.

Uncertainty Modelling Methods
To obtain enough information, many uncertainty modelling methods have been applied to depict RE output with considerations of uncertainty in the operation problem. The most common approaches are probabilistic approaches, wherein uncertain parameters (for example, output of wind unit and photovoltaic) are modeled by probability density functions (PDFs) and dealt with probabilistic strategies such as Monte Carlo simulation (MCS) [81,82], scenario-based analysis (SBA) [83], and point estimate method (PEM) [84]. The core ideas of these three common methods are described as follows.
• MCS utilizes PDFs of input uncertain variables (i.e., wind speed) to calculate objective functions on the basis of input-output relationship and select the best individual at different iterations as optimal solution. The more the number of trials is, the more accurate result is achieved, but computational time increases. Overall, MCS is simple to implement but computationally expensive. • SBA method subdivides the PDF curve of uncertain input variables into several regions. The probability of falling each region can be calculated and expressed based on PDF and then the expected value of output is represented with the total sum of the product of probability and corresponding value in each region. Likewise, higher number of scenarios increases the accuracy of the achieved results. Hence, SBA method is also a technique with high computational burden. • PEM works based on the statistical information provided by the first few central moments of an input random variable on K points for each variable. With these points and the function, which reveals the relation between input and output variables, the information about the uncertainty associated with problem output random variables can be obtained [85]. Similar to MCS, PEM method uses deterministic routines for solving optimization problems with uncertainty. As it provides an approximate description of the statistical properties of random variables, it can significantly reduce the calculation times. However, the PEM methods cannot describe the correlation between different uncertain inputs.
These conventional probabilistic approaches are conducted based on the assumption that the random variables accord with the given PDFs (for example, wind speed prediction obeys the Weibull distribution [86], solar irradiation obeys beta distribution [87]).
With the increasing penetration of renewable resources in IES, more renewable energy generation units with different power generation characteristics in time and space scale may be integrated in the system, and these existing PDFs are not fairly applicable and sufficiently accurate. To solve that, some data-driven improved probabilistic approaches are imposed. In this paper, they are subdivided into three categories as follows.
(1) Non-parametric probability distribution models Non-parametric probability distribution models are with no restricting assumption on the distribution of random parameters but abstract statistics information from samples, and hence they can obtain more accurate probability distribution of random variables. Ref. [88] adopted non-parametric optimization method to operate and plan energy storage units in power distribution grids. The results were compared with conventional deterministic and parametric stochastic approaches based on Gaussian approximation, and it was proved to significantly release computation burden. Kernel density method was proposed in Ref. [89] for wind speed distribution model, and it was demonstrated to be more accurate and more superior in adaptability using three years' actual wind speed data at ten wind farm sites. Based on that, Ref. [90] presented an improved diffusion-based kernel density method (DKDM) to estimate wind speed probability distributions with consideration on both band width selection and boundary correction. The statistics of goodness of fit tests of DKDM were smaller than the conventional methods in most stations, and its practicality was demonstrated.
(2) Stochastic process model Stochastic process model is normally adopted in prediction of renewable resources output, and recently it is extended to model uncertain parameters in dispatch optimization by some researchers. Ref. [91] quantified the uncertainty of uncertain wind power in stochastic economic dispatch with the proposed method, which was on the basis of Karhunen-Loeve expansions and historical data at each wind farm. The experiments showed that the proposed method was one to two orders more efficient than Monte Carlo-based estimates. However, it faced a curse-of-dimensionality challenge.
(3) Scenario generation methods based on artificial intelligence technologies Scenario generation methods based on artificial intelligence (AI) technologies take the advantage of strong data mining ability to directly generate uncertain scenarios with actual data instead of fitted probability distribution density functions. Hence, they are more accordant with the real characteristics of random variables. Ref. [92] utilized generative adversarial networks to capture renewable energy production patterns in both temporal and spatial dimensions for a large number of correlated resources. In Ref. [93], artificial neural networks (ANNs) were proposed to generate scenarios to account for various stochastic variables including electric load, PV, and wind production, and it applied this method to an optimal participation problem of a PV agent in day-ahead market. The percentage difference of ANN-based model with respect to perfect forecast could be up to −0.78%, −1.64%, and −4.47% in three scenarios, and it showed its accuracy compared with other methods.
Compared with conventional probability distribution density functions, these datadriven methods depict the variation characteristics of renewable resources more accurately and hence improve system operation strategy. Nonetheless, the performance of these methods strongly relies on the external manifestation of samples, which may influence the results of dispatch. Additionally, these methods are investigated in the conceptual stage and most are applied in simple power system to validate their effectiveness. As for an IES with plenty of different variable types and complex coupling relations, the size of optimization problem is much larger, and their shortness in computation efficiency problem will aggravate.
Another normal uncertainty modelling method is robust optimization method (RO). The conventional RO considers the worst-case condition and yields a conservative scheduling in which the uncertain variables allow change within a given range around their forecasted values. The dispatch problem based on RO can be expressed as Equation (5), where f (x) is the objective functions, x is the robust feasible solutions of the decision variables, A and b represent the constraint vector in robust model, and D is the uncertainty set.
It should be noted that the set of range D significantly influences the optimization results. At present, the most commonly utilized uncertainty set is polyhedral uncertainty sets, as shown in Equation (6). The subscript of ω denotes uncertain variable d, which constitute a set Ω, d ω represents the value of d, and D l ω and D u ω are the upper and lower limit of d ω .
It is easy to understand and implement, yet the simple consideration on upper and lower limits cannot truly reflect the variation range of uncertain variables and furthermore leads to a conservative solution from mechanism. To face that, some other uncertainty sets are introduced, represented by ellipsoidal uncertainty sets [94,95] and extreme scenarios based uncertainty sets [96][97][98]. Though these methods reduce the conservativeness and guarantee the robustness of optimal decision scheme, they are still on the basis of worst scenarios, and the over-conservative problem exists. Probability distribution of random variables is combined into the RO method to make up for the deficiency. The core idea of data driven robust optimization (DRO) methods is to establish an ambiguity set based on probability distribution with data and make the optimal decision under worst scenario in ambiguity set. By fully utilizing the statistical information of uncertain parameters, the conservativeness can be reduced. It should be noted that the key point of this method is to establish the ambiguity set. The two most commonly utilized methods are statistical moment based [99,100] and distance based [101] ambiguity sets. Ref. [99] demonstrated that the conventional RO method was too conservative, while the DRO method could greatly reduce the redundant cost by over-conservativeness. Ref. [101] designed a distancebased ambiguity set to capture uncertainty of wind power distribution in a two-stage UC problem. The results showed that the cost decreased from $4.9131 × 104 with data size 50 to $4.6850 × 104 with data size 5000, and that the conservativeness of the problem could be controlled by adjusting the data size and confidence level. However, the average calculation time was over 5700 s for a practical system. Overall, the statistical moment based DRO methods reflect the statistical information of random variables through data analysis, which is easy to implement and achieve. On the other hand, it cannot fully mine the statistical information and leads to a respectively conservative result, while the distance-based method builds empirical probability distribution based on data and defines the distance between the actual and established empirical one. It needs fewer data and thoroughly excavates information. Nevertheless, the optimization problem may have difficulties in solving, and even cannot be solved with the expansion of the scale. To have an intuitive comparison, the advantages and disadvantages of these above uncertainty modelling methods are summarized in Table 2. Additionally, some other stochastic methods including possibilistic methods, hybrid probabilistic-possibilistic methods, interval optimization (IO), and information gap decision theory (IGDT) are also extensively employed in handling uncertainties of RE generation output in IES operation optimization. The specific description of these methods can be found in Refs. [102,103], and they are not deeply reviewed in this paper. In general, all these methods are trying to seek an expression form to quantify the effect of uncertain inputs on the optimization results. Their main difference reflects on the uncertainty depiction techniques. Based on these methods, the uncertain variables are represented with a stochastic form in objective equations and constraints equations, and the non-deterministic dispatching problem is usually converted to a deterministic one so that the adverse effects brought by uncertainties can be reduced.

Multi-Stage Optimization Strategy
In uncertainty modelling methods, the information of RE output forecast values are estimated by probability distribution or other forms before the operation horizon begins and the uncertainty level is contained to be used in later system operation optimization. These methods are efficient on the condition that the RE penetration is relatively low, and hence the negative impact of deviation between predicted and real values is not great on the IES as well as upstream main grid. For a large-scale RE incorporated IES, the error between expected and actual values become larger and then may affect the results accuracy. Moreover, the impact of IES on the upstream main grid is expected to be minimized to reduce penetration and release power grid pressure. Thus, more efforts are required, especially for a renewable-rich IES.
It is obvious that the prediction error is getting smaller as time goes on. Multi-stage operation optimization is therefore proposed to successively mitigate uncertainty of RE at different time scales and then reduce adverse influence on the IES optimal operation. Generally, the process can be divided into three stages including day-ahead dispatch (DA), intra-day scheduling (ID), and real-time control (RT) corresponding to timescale as illustrated in Figure 5. DA is the first stage to determine the operation strategy of each device in IES with 24 h forecasting data of loads and REs output obtained at time t − ∆t, wherein t presents the DA dispatch start time and usually sets as 0:00, and ∆t denotes the time delay due to communication or calculation. The dispatch horizon of DA stage is 24 h (∆T da ) and corresponding resolution is 1 h (∆T 0 ). Model predictive control (MPC) based rolling optimization has been extensively applied in intra-day scheduling [104,105]. In MPC-based ID stage, the time horizon is ∆T id , the start time is t 0 , the corresponding resolution is ∆T 1 , and the rolling window is ∆T r . In the rolling process, only the first interval scheduling plan is implemented, and the optimization process will be repeated at time t 0 − δt + ∆T 1 . The time-scale of ID is always 15 min level, while RT is always in the time scale of second/millisecond. The control resolution is denoted as ∆T 2 , and dt is the time delay of measurement devices. Each optimization stage has different functions and objectives, as illustrated in Figure 6. In the DA stage, based on day-ahead input information and devices' physical models as well as the energy flows model, the optimal decision is made. Since it is a relatively long-term scheduling, it always pursues economic or environmental performance maximization. The DA stage provides a base operation decision of variables for the following ID stage, and the uncertainty of RE can be considered, utilizing uncertainty modelling method. After that, ID strategy is optimized with the almost same process as DA but based on short-term prediction information. It should be noted that the ID stage aims to adjust the optimization strategy to eliminate the adverse impact from RE generation output prediction error between long-term DA stage and short-term ID stage. Hence, the optimization objective is always to maintain the IES, a controllable individual from the perspective of upper power grid. With the time approaching, the main characteristics of RE generation output are captured, and the operation strategy obtained can mostly satisfy the actual condition. However, deviation between predicted and real values still exists. As the electricity power is real-time matching, some fast fluctuations may sometimes cause safety problems. Hence, the RT stage optimization is usually integrated from the perspective of reactive power adjustment to mitigate the voltage fluctuation caused by RE variability utilizing real-time data to guarantee the network security. Through coordinating optimization objectives in multiple operation optimization stages, the renewable resources fluctuation and uncertainty can be eliminated step by step.
The multi-stage optimization strategy is originated from the power system and was applied to IES operation optimization recently. Ref. [107] proposed a two-stage operation optimization scheme, wherein the first stage was based on probability analysis before the renewable output is not realized, and the second stage was set after the uncertainty has been observed in a high-penetrated distribution system. The results demonstrated that the energy loss could be reduced from 110.3 to 54.9 kWh, and the voltage violation could be eliminated. Based on that, a three-stage operation optimization scheme consisting of DA, ID, and RT was presented for a renewable-rich system to promote RE consumption and mitigate the system uncertainty in Ref. [106], wherein the dispatch horizons for multiple  Faced with system integrated multiple heterogeneous energy resources, Ref. [108] conducted a two-stage operation to smooth out the RE fluctuations and follow the load variations in an electricity-cooling coupled IES system. In the day-ahead stage, hourly generation schedules of each MG component based on the forecasted electricity and cooling demands of the next day were determined with expectation of minimizing operation cost. The uncertainty of RE was characterized to provide reference for the following rolling process. In the real-time stage, necessary generation was adjusted on the basis of wind and solar short-term power predictive values to reduce the discrepancy of exchanged power with the main grid and minimize possible penalty to IES. Ref. [109] formulated a day-ahead and real-time coupled energy optimization model for a heat-electricity IES and proposed an event-triggered-based distributed algorithm to solve the distributed optimization problem. The real-time dispatch was divided into hourly timescale dispatching (1 h) and minute timescale dispatching (15 min) according to real-time dynamic variations of heat and electricity power features. The former was to make the intra-day dispatching results track the DA results and to accommodate variations of heat load along with solar water heater, while the latter was aiming at accommodation fluctuations of RE and electricity loads. In Ref. [110], a multi-time scale framework of coordinated optimization strategy was presented including multi-objective DA dispatch, two-layer ID dispatch, and emergency control aiming at a typical IES. To consider dynamic characteristics of heterogeneous energy resources under different time scales, the ID dispatch was divided into intra-day cooling-heat dispatch and intra-day electricity dispatch with time span of 2 h and 1 h, time interval of 1 h and 5 min, respectively, in the rolling process. The case results showed that the multi-stage operation optimization could effectively stabilize fluctuations of wind and solar generation and had the least cost compared with conventional strategies.
From the existing work review, it can be concluded that the application of multi-stage operation optimization in IES is still in the relative initial stage. Most literatures adopt the same research routine as power system and subdivide the stage into two or three stages, as aforementioned.

Challenges and Prospects
In IES, renewable resources are a significant part, and the ratio is increasing. The uncertainty and fluctuation characteristics of RE may cause a cascade effect through strong couplings among energy flows, and devices bring challenges to IES operation strategy making and increase the interactions with the upper grid to maintain normal operation. To relieve the upper grid burden and promote IES independence, effective handling methods with RE uncertainty are significant. The latest work is reviewed from the uncertainty modelling methods and multi-stage optimization strategy perspectives. Based on that, the existing challenges and prospects are given as follows: (1) Deep study in the novel uncertainty modelling methods. To solve the limitations in conventional methods, some novel data-driven approaches are proposed and are demonstrated to be effective. However, these methods are mostly applied in a simple case and seldom extended to the IES scenario. Accounting for the complex coupling relations and various variable types in IES, some work needs to be further conducted.

•
Novel uncertainty modelling methods with consideration of temporal-spatial correlation of different RE generation output.

•
Novel uncertainty modelling methods with consideration of coupling of various types of variables, such as continuous and discrete types.

•
Deep sensitivity analysis of the impact of some crucial parameters on the dispatch results and how to select them or optimize them when facing different demands. For instance, the uncertainty set D is an important parameter that influences the conservativeness of DRO method directly. Most work sets it according to experiences and there is a lack of theoretical study.
(2) Stage subdivision accounting for multiple energies flows characteristics in IES. The majority of literatures adopt the same dispatching interval as the power system. More concretely, 1 h resolution for the day-ahead, 15 min interval for the intra-day, and 15 s resolution for the real-time stage. However, different energy systems have great difference in dynamic characteristic, equipment control, network response process to the dispatching command, and load pattern. For instance, power generation units need to carry on intraday schedule on a 5 to 15 min interval to accommodate intermittent resources, while heat load is respectively stable in the short time-scale, and consequently it is unnecessary for heat units to regulate as frequently as power generation units. Frequent movements of devices have adverse impact on both unit life and operation economy. Therefore, further studies are essential on the optimum dispatch period selection for IES optimal operation with considerations of distinguishing dynamic characteristics of energy flows.

The Time-Scale Characteristics of Different Energy Flows
Electricity, heat, natural gas, and other energies are integrated and comprehensively utilized in IES to enhance efficiency on the whole. They are tightly coupled as they share the same transport paths and are transferred to others from production to utilization process. Remarkably, the dynamic characteristics and topologies of multiple energies are significantly different from each other. Electricity travels at light speed with little time delay and energy loss, and hence it is possible to be transported over long distances, while gas and heat transport at a slow speed through pipelines and are applicable for short-distance transportation to reduce friction loss. Additionally, though the dynamic characteristics of gas and heat are similar to some degree, there remain some unneglected differences in practical process. Different energy flow characteristics are summarized in Figure 7, and they can be attributed to response timescale differences in essence. The different dynamic characteristics of hybrid energy flows are non-negligible. The transmission delay, losses, and virtual storage characteristics under complex topology of heat/gas networks change heat/gas balance between supply and demand, and consequently affect the electricity balance even further to the whole IES operation through coupling devices. Compared with a single power system, the multi-timescale characteristics of hybrid energy flows make them nontraditional and challenging in IES operation optimization. On the other side, the different dynamic characteristics of these energy flows can be regarded as a potential energy storage option to coordinate the real-time supply and demand relation and guarantee the feasibility of optimized solutions. These characteristics can be utilized to exploit potential in decoupling restriction among multiple energies. Compared with introducing additional storage devices, the potential storage capacity in existing pipelines is almost costless. Hence, the adjustment of fluctuating wind power, as well as volatile load, can be improved by concerning flexibility, and then the independence on the upper grid and benefit can be promoted.

Thermal Dynamic Characteristics Considerations
Heat and electricity are usually provided simultaneously by cogeneration or trigeneration units to improve energy utilization efficiency, and their coupling presents a typical IES. For the strong coordination between power generation and heating supply, many researchers investigate the combined operation optimization of integrated electricity and district heating systems (IEHS). District heating networks (DHNs) are the major component of a DHS, which deliver heat from generation sources to users through thermal mediums (usually flowing water).
The general formulation for the joint of power system and district heating system is similar to the process in Section 1, but with additional heating system constraints and modelling of IEHS. The additional constraints include heating network constraints, node balance constraints, heat loss constraint of the pipe section, pipe input and output constraints, heating power output constraints, etc. The modelling of IEHS is categorized based on whether dynamic characteristics of the heat transfer process in DHSs are considered and how to depict it.
In general, researchers have developed data-driven methods and model-driven methods. For the first category, the temperature quasi-dynamics of pipelines are described with a black box, which is estimated and validated by massive data. Ref. [111] proposed a black-box compact physical model to calculate IES power flow, and Ref. [112] established an individual model to depict the dynamic heat process. Data-driven methods are easy to implement, yet need numerous data in training and are not general in different networks structures or operating conditions. Model-driven methods are based on characteristics of heating system and physical knowledge, which are readily adaptive to changes in system and possess more accuracy. They can be further subdivided into three main classifications according to model complexity. That is, the steady-state method [113][114][115], node-method [31,[116][117][118], and delicate dynamic network modelling methods [99][100][101], ranked from low complexity to high.
The heat transfer process in DHN steady-state model is described holistically as Equation (7) shown. T in and T out signify the inlet and outlet temperature of a pipe; T a represents the ambient temperature; m is the mass flowing into the pipe; L is the pipe length; c p denotes the specific heat of water at constant pressure; λ denotes heat transfer coefficient of a pipe per unit length; m is the water mass flow rate in pipeline. In the steady-state model, the temperature drop along pipelines is described, but the temperature distribution in DHN and time delay effect are neglected. Ref. [113] exploited the flexibility of DHSs with a steady-state model in IES operation stage. Aiming at profit maximization in heat-power spot market, Ref. [114] established both the hydraulic and thermal model of DHN and proposed a two-step decomposition algorithm to solve the nonlinear and nonconvex optimal problem. Ref. [115] detailed the supply and return networks, and additionally considered the variation of mass flow rates as well as flow directions in a heat-electricity IES economic dispatch problem. The comparison showed that the cost was reduced, but the network limits were violated in several periods if the heat transfer process was ignored.
The second category, node-method, is a commonly utilized dynamic approach to describe temperature distribution in pipelines. Its principle is to estimate the outlet temperature preliminarily from the inlet temperature with regard to flow time from one node to another and then evaluate the heat loss of the pipe as well as outlet temperature. The schematic diagram is illustrated in Figure 8, and the pipeline outlet temperature of mass flow is calculated by Equation (8). The DHNs are finally embedded as a set of mixedinteger constraints into the combined dispatch model. In Ref. [116], a general dispatch method considering heat dynamics of DHNs under quality regulation mode was built up, wherein the transfer process constraints were emphasized and an iteration method was proposed to handle nonlinear constraints. A comparison with the case ignoring the DHN storage capacity, considering both storage capacity and detailed heat transfer process, and only considering storage capacity was conducted, and the coal use values were 3383.4, 3372.4, and 3348.7 ect. The results demonstrated that the slow characteristics of heating could reduce cost and provide additional flexibility for IES dispatch, whereas the flexibility might be overestimated, and the optimization results were not practical the heat transfer process was not fully considered. Ref. [117] established a fully nonlinear model of DHN under quantity regulation mode to exploit storage potential and solved by a deterministic method. Based on that, a linear model of DHN was proposed under quality regulation in Ref. [118] to relieve the calculation burden brought by various constraints and variables in heat transfer. Ref. [31] utilized the thermal inertia of DHN to improve system flexibility for wind power integration, and quantitatively evaluated the thermal state of DHN to present its dispatch potential.
These blocks represent a sequence of water masses that flow into the pipeline at consecutive periods. ms t represents the water mass flow rate, ∆t is the time step, and ms t ·∆t denotes the mass of water flowing into the pipe at period t. t − n and t − m, respectively, denote the index of the last period whose mass outflows the pipeline before period t and period t − 1, wherein n and m are corresponding time delay. W represents the total water mass contained in the pipeline, while X represents total mass flowing into the pipeline from period t − n to t, and Y represents that from t − m to t. With the above equation, the outlet temperature T out,t at current time step t can be estimated by the average temperature of the dark red part weighted by water mass. For the third category, delicate dynamic network modelling methods, they fully consider heat dynamics in DHN with multiple differential equations. These dynamic methods have been widely studied; nevertheless, the detailed information of DHNs, such as topology data and pipeline parameters, greatly aggravates computational complexity in operation optimization. Hence, most studies are conducted with many simplifications or assumptions and meanwhile devoted to exploiting a feasible solving method. Driven by this, some equivalent methods are derived to represent the complicated physical model and release computational stress in recent years. Ref. [119] proposed a water mass method for pipeline thermal inertias by removing the integer variables and differential equations to make the combined dispatch model tractable. In Ref. [120], the thermal inertias of DHNs were represented by an outer approximation from likewise the state-of-charge formulation of thermal storage devices or batteries. The additional savings of several thousand Euros per day were achieved using the thermal inertia of a heating grid as storage. Ref. [121] presented an equivalent representation of DHNs from electrical-analogue perspective, wherein the input information from DHS was regarded as a black-box. Comparison experiments were conducted in the cases of no consideration of heat transfer dynamics and with consideration of heat transfer dynamics (the pipeline transfer delays were 24 and 48 min). The wind curtailment was reduced considering slow dynamic characteristics of heating, and the slower the storage, the greater was the potential. These equivalent methods are efficient in the research case, but generality is not sufficiently strong.
Moreover, the thermal inertias in the end-user side are also studied to explore the potential of heat inertia in overall efficiency improvement. Ref. [122] proposed a virtual energy storage system model based on the characteristics of buildings for the heat load. In Ref. [123], a building indoor temperature calculation model was established, and the response ability of buildings in dispatch was relevant to the set indoor temperature boundary. With simultaneous consideration on transmission delay in DHN and the storage capacity of buildings, the cost could be reduced by 9.3% and the wind power consumption was improved by 44.6%. Ref. [124] detailed the buildings' heat load model with additional consideration of building characteristics' diversity as well as outdoor ambient temperature variation. To release calculation burden and protect privacy, Ref. [116] introduced an equivalent building model based on entransy dissipation-based thermal resistance theory, and Ref. [125] proposed synchronous response model of buildings on the basis of the first-order equivalent thermal parameter (ETP) model. Considering the influence among adjacent buildings under different regulation modes, Ref. [126] studied the different thermal capacity of buildings and their impact on IES operation optimization separately. Similar to DHNs, physical models with some simplifications of buildings' heat load are adopted. Generally, the concrete heat transfer between pipelines and indoor air is ignored, and only the indoor temperature variation is concerned to indicate actual heat load and potential for heat storage. Additionally, it is noted that most existing work assumes the buildings are with the same structure and response behavior. That is, the comfortable indoor temperature range is fixed among all buildings and all time periods, and different buildings as virtual energy storage show identical charge/discharge state. It may be a future research direction for higher practicability.

Gas Dynamic Characteristics Considerations
The natural gas is environmentally friendly and effective with no production of SO 2 and substantially less NO x compared with other fossil fuels. Hence, it promotes the utilization of some gas-related furnaces such as gas-fired generators and gas boiler, as well as the development of some new technologies like power-to-gas (P 2 G). Systems containing these devices and technologies impose a linkage between gas and electric power transmission system and hence can be regarded as an IES.
Conventionally, power systems and gas systems are separately considered, and it is assumed that the supply from the gas network is always well sufficient. On this condition, the IES operation optimization is conducted, neglecting gas flow characteristics as well as their impact on the whole system. It simplifies the optimization process; nevertheless, it is not practical with the increasingly tighter interconnection between two systems. Specifically, the power system flexibility is limited to gas networks' ability to meet requirements of gas-fired units and fast change of gas load. Facing that, many scholars focused on two systems' interdependence assessment. They researched the intrinsic characteristics of gas networks evaluation, and determined possible limits in terms of constraints on operation flexibility and security of power system during a period of time in the operation optimization process. Ref. [127] proposed an integrated model to assess the impact of natural gas networks on power system operation, and incorporated these security constraints into the unit commitment optimal problem with the objective of minimizing cost. The natural gas networks model was represented by various bound constraints on pipelines, plants, contracts, and so on. Ref. [128] characterized natural gas transmission and then studied the system operational flexibility under these security constraints. On this basis, Ref. [129] and Ref. [32] extended it to the economic dispatch of a system with more energy forms. The former targeted a CHP-based system and proposed a general model to schedule all energy sources with considerations of gas losses. As more factors were involved, which complexed the problem, an advanced PSO technique was adopted to solve it more efficiently. The latter studied a system with both gas-fired generation and heating sectors. For the strong coupling between these furnaces and gas networks, detailed discussion was conducted on system operation flexibility under different heating scenarios.
The above studies emphasized the limits brought by the interconnection between two systems and incorporated them into the optimization model to find the optimal solution that can meet security constraints. However, the slow response characteristic of the gas system is neglected. Figure 9 shows the framework of a natural gas system and transmission system. Gas is transported from the source to the load side through long pipelines with a slow speed. The large inertia decides that gas cannot realize real-time balance like electricity. With increasing RE incorporation, gas-fired units are utilized to fast adjust load conditions to accommodate more RE, and it in turn leads to a fluctuation of the gas load. If the slow dynamic characteristics of the gas system are disregarded, the dispatch results are not accurate and practical. Hence, some scholars tried to depict the gas network model and schedule two systems jointly from a holistic angle.
Attributed to the Weymouth function, the general steady-state equation for gas flow models can be illustrated as Equations (9) and (10), which is widely utilized in most researches. f mn represents the gas flow value from node m to node n. ω m and ω n are the gas pressure at these two nodes; sgn(ω m − ω n ) represents the gas flow direction; "1" means from node m to node n, while "−1" means from node n to node m; C mn is a constant, which is determined by the pipeline characteristics and natural gas compositions. Ref. [130] represented the natural gas reservoir storage capacity in a given bound and modelled the gas flow in a steady-state way in a co-optimized problem. Likewise, based on the steady-state relation of gas pressure and flow, Ref. [131] estimated the upward and downward line pack utilizing nodal pressure limits, and optimized the system operation with more energy forms (combined cooling, heat, and power) and more objectives (economic benefits, safety, and efficiency). Steady-state gas flow models are easier to implement and bring less challenge in calculation, but they cannot depict the dynamic process and are demonstrated to be less accurate in weighing the flexibility of the gas system to provide for gas-fired generation units in case of a contingency [132].
Driven by this, some dynamic network models are developed. Ref. [133] established a detailed dynamic gas flow model with multiple partial differential equations (PDEs) by applying the laws of conservation of mass and energy, and the laws of momentum, to characterize the important parameters tendency with time and position. It also compared the results with the steady gas flow model, and found that the results with the steadystate and transient-state gas models yielded two distinct results. Hence, it is necessary to consider gas dynamics to obtain practical results and suboptimal schedules. To handle the challenging computational tractability of numerous PDEs, Ref. [134] approximated the complex PDE constraints with the reduced network flow (RNF) and studied joint optimization problem for different coordination scenarios. Comparison experiments were conducted, and the results showed that the steady-state gas transmission model would cause pressure violations, while the consideration of gas slow dynamic characteristics could remove pressure violations, but the generation dispatch cost was increased.
Additionally, the gas slow dynamic characteristics make it possible to store gas in large quantities within pipelines, and the potential storage capacity can be explored in short term or some contingencies (i.e., insufficient production, congested networks, or price fluctuations) to guarantee the adequacy and improve system flexibility. The gas stored in pipelines in a period is known as line pack. Except for economics improvement, coordinated operation optimization is also explored with considerations of RE variability to fully utilize the line pack as a buffer and guarantee system reliability. In general, increasingly high penetration of intermittent RE inevitably leads to frequent fast condition change of gas-fired units and then affects gas demand from networks, which may cause security problems in both systems. Based on the steady-state gas flow model, Ref. [135] and Ref. [136] took the uncertainties of RE generation into the co-scheduling model to ensure the optimization results were within the security boundary, wherein the interval optimization and stochastic approach were adopted to depict the wind uncertainties. In Ref. [137], a two-stage coordinated dispatch problem was proposed based on distributional robust optimization model, and the line pack was utilized to improve system operation flexibility. Despite the capability of line pack in gas networks, Ref. [138] also took the P 2 G technology into account to alleviate the effects of RE in the operation of gas and electricity systems. The gas line pack was as a buffer, and the PtGs transformed the surplus energy of wind generations to gas network during peak times. The steady-state gas flow model is still the most extensively adopted at present for its better operability. Dynamic models are proven to be more accurate. However, these precise modelling methods based on multiple PDEs involve numerous variables and more dimensions. They are mainly solved by various finite element numerical methods with a huge calculation burden. Additionally, the analyses of gas and electricity networks are not unified and impose a barrier on the integrated operation optimization problem.

Challenges and Prospects
An IES incorporates many forms of energy flows with distinguishing characteristics via multiple energy carriers, which provide a potential buffer. It can be regarded as an effort direction to increase IES operation flexibility and promote individual independence as well as benefit for a future RE dominated energy system. The considerations of distinguishing characteristics of multiple energy flows bring new challenges and difficulties to the original IES operation optimization problem. From the literature review, it is noted that the majority focus on coordinated optimization of either heat-electricity IES or gaselectricity IES with considerations of time-scale characteristics of different flows. They are devoted to modelling the transmission in networks and incorporating them into the dispatching model.
Studies considering both heat and gas inertia have emerged, yet are still in their infancy. Ref. [139] optimized operation of an islanded electricity-heating-gas IES with the objective of minimizing operation cost, and the storage capability of gas pipelines was demonstrated to be effective in increasing system scheduling flexibility against RE uncertainties. It also verified that the results might not be optimal if neglecting the built-in storage capabilities of pipelines and the slow travelling of gas flows. However, the heat inertia was ignored. Ref. [140] established both the thermal and natural gas network model to consider practical constraints and their impact on results. The electricity-thermal-natural gas networks were steady-state based and the energy flow characteristics of heat and gas were disregarded. Ref. [141] considered both heat and gas dynamic characteristics with different inertia time constants in the optimal operation model establishment, and adopted a linearization method to make the problem tractable with the objective of minimizing the total operation costs.
On the whole, the most challenging issues in existing work are formulating hybrid energy flows and solving the nonlinear nonconvex optimization problem in a tractable and efficient way. A detailed dynamic model to depict gas or heat inertia is conducive to improving the rationality of optimization results and providing potential to increase system operational flexibility. However, the networks bring plenty of limits and multi-type variants into the optimization model, which increase computational stress and difficulty greatly. To handle that, the existing work concentrates on assuming simplifications like the steady-state model, or developing a linearization approximation, or seeking an equivalent model. Nonetheless, they are merely applicable to an identified problem and lack of sufficient generality and extendibility.
In future, a more detailed and unified dynamic network model for both gas and heat is necessary. Additionally, a more efficient solving method needs to be exploited simultaneously. An equivalent representation may be an alternative and feasible way to ease contradictions between model accuracy and computation stress, wherein its reliability and generality must be considered. Additionally, the existing work cannot realize unified analysis of different energy flows and leads to barriers in IES operation optimization. Intuitive and unified modelling method and theory analysis framework of hybrid energy networks can be further explored in future. In terms of joint optimization, it is critical to exploit a reasonable time resolution for each energy system in the dispatch modeling to properly coordinate the difference in dynamic characteristics.

Interest Coordination among Multiple Participants
One of the main reasons for IES' strong dependence on the upper grid is the volatile multi-energy loads. When the energy demand exceeds the supply capability of IES, the system has to ask for help from the upper grid to keep operation reliability. Since the IES is always close to users and the multiple energy is volatile, the interaction between the upper grid is frequent and increases its adjustment pressure. In a smart grid scenario, IES can manage the load demand effectively based on the formula of high consumption-high payment during peak hours and low payment for the valley periods. Since the load pattern is smoothed, the single IES can better satisfy the demand side without or with less external help. The energy price from the upper grid can also be formulated according to actual supply-demand state. In that case, IES can act as both a price-taker and a price-maker. From a union effort prospective, the coordination among these participants provides potential to enhance IES independence without reducing individual benefit.
With the development of information and communication technologies, a large-scale IES is conceived to consist of many lower-level sub-IESs under the interconnected energy internet framework in the future vision. As shown in Figure 10, an IES is equipped with energy generation and conversion devices to satisfy diverse end demand, and it can also obtain energy supply from the upper public network such as power grid and natural gas network. Beyond a certain IES, there are many parallel autonomous IESs with diverse size and configurations. Though these sub-IESs can act as a self-regarding entity in terms of operation strategy decision and price responsiveness, the decision of an individual IES is usually tightly coupled with others' choices. The reasons are as follows. The local market retail energy prices depend on both the time of energy consumption and the total demand of all IESs to provide better incentives for energy users to adjust their usage pattern and alleviate network pressure. Hence, the decisions of each IES can affect the market prices, and the purchase cost will affect other IES's operation strategy decisions. Additionally, the limited transmission capacities of networks lead to limits on the supply of upper gas/electricity energy systems to these sub-IESs. That is, the total dispatchable energy of all sub-IESs is constrained globally. For the aforementioned interconnections in the IES cluster, it is crucial to consider the impacts from other neighboring IESs when determining its own operation strategy. Since each pursues a better interest, it poses the emerging problem of how to coordinate all the participants who have effects on the final decision.
Game theory is introduced to provide resolution for interest conflicts among these interacting parties and seek an optimal result. Game theory is a formal analytical as well as conceptual framework with a set of mathematical tools enabling the study of complex interactions among independent rational players [142]. The main components involved in a game model include the set of players, the action sets, and the utility functions. Each player adjusts their own strategy to maximize their objectives in accord with others' choices. The equilibrium is achieved when no player can benefit more by unilaterally changing their strategy.
The non-cooperative game model represented by Stackelberg game has been extensively utilized in the interest coordination among multiple participants for designing operation strategies, especially with the trading process between energy provider and users. It stresses the autonomous decision of one individual behaving selfishly to seek their own interest maximization. The basic Stackelberg game model Ψ is displayed in Equations (11)- (13).
• ℵ ∪ is the player set. ℵ is the set of leaders, while denotes the followers. They represent all the participants associated with the problem. For instance, IES operator and lower energy users.
• s n and S k are, respectively, the strategies space of each leader and follower, which are feasible sets of actions that the players can take. For example, the scheduling strategy of generation units of IES and the user's energy usage plan. • W leader,n and U f ollower,k denote payoff functions, which are the objectives each player pursuing to evaluate the efficiency of selected strategy, such as profit.
and configurations. Though these sub-IESs can act as a self-regarding entity in terms of operation strategy decision and price responsiveness, the decision of an individual IES is usually tightly coupled with others' choices. The reasons are as follows. The local market retail energy prices depend on both the time of energy consumption and the total demand of all IESs to provide better incentives for energy users to adjust their usage pattern and alleviate network pressure. Hence, the decisions of each IES can affect the market prices, and the purchase cost will affect other IES's operation strategy decisions. Additionally, the limited transmission capacities of networks lead to limits on the supply of upper gas/electricity energy systems to these sub-IESs. That is, the total dispatchable energy of all sub-IESs is constrained globally. For the aforementioned interconnections in the IES cluster, it is crucial to consider the impacts from other neighboring IESs when determining its own operation strategy. Since each pursues a better interest, it poses the emerging problem of how to coordinate all the participants who have effects on the final decision.  Game theory is introduced to provide resolution for interest conflicts among these interacting parties and seek an optimal result. Game theory is a formal analytical as well as conceptual framework with a set of mathematical tools enabling the study of complex interactions among independent rational players [142]. The main components involved in a game model include the set of players, the action sets, and the utility functions. Each player adjusts their own strategy to maximize their objectives in accord with others' choices. The equilibrium is achieved when no player can benefit more by unilaterally changing their strategy.

IES
The non-cooperative game model represented by Stackelberg game has been extensively utilized in the interest coordination among multiple participants for designing operation strategies, especially with the trading process between energy provider and users. It stresses the autonomous decision of one individual behaving selfishly to seek their own interest maximization. The basic Stackelberg game model Ψ is displayed in Equations (11)- (13). When no player can increase his utility by choosing a different strategy other than s * n /S * k , the Nash equilibrium is obtained. s * n and S * k represent the optimal strategy of leader n and follower k. s * −n and S * −k represent the optimal operation strategy of the leaders and followers except for leader n and follower k. s * n = argmaxW leader,n (s n , s * −n , S * k ) Stackelberg game clearly differentiates all the participants into leaders and followers according to the sequence of their actions [143]. The leaders have the priority to make decisions and the followers adjust their strategies based on that. Several leader-follower models have been developed for analyzing the interactions between the energy suppliers and consumers, which can be categorized into single-leader multi-follower (SLMF) and multi-leader multi-follower (MLMF) models. Ref. [144] introduced a multi-party operation optimization framework based on SLMF Stackelberg game theory for an IES consisting of CHP operator and PV prosumers. The results demonstrated that the model could effectively determine the time-of-use (TOU) prices scheme and optimize the net load characteristics. Ref. [145] adopted SLMF model in a bi-level IES programming model including upper energy generation level and lower load level. The optimal dynamic pricing scheme was obtained to motivate consumers to participate in load profile adjustment. SLMF model is always adopted in single-energy provider with consumers. Compared to SLMF model, the MLMF models could better simulate the real energy market transactions. In Ref. [146], the MLMF model was utilized for analyzing the multiple energies trading problem in IESs. A number of distributed IESs lead the game and decide the unit prices of energies they generated, while multiple energy users perform as followers. The profits of all participants were optimized to realize maximization in a competitive energy market. The results showed that even one participant was provided more resources, its profit was almost unchanged. That is, the equilibrium strategies of all participants were obtained. Ref. [147] studied the optimal operation of all market participants in a three-level IES, which was composed of utility companies, IESs, and users. The optimal operation of each market participant was obtained under a dynamic pricing electricity mechanism. The results showed that the profit tends to decrease with the increment of energy hubs, and therefore IES should devote to raise its market competition so as not to be crowded out. Additionally, with the rational interest coordination among different participants, the load profile could be flattened, which was a benefit for intermittent renewables generation accommodation. In Refs. [148,149], a non-cooperative Stackelberg game was utilized to solve the duck curve to increase solar penetration. Likewise, Ref. [150] identified the best interactions between a grid and buildings with Stackelberg game. The results showed that the net profits increased by 8% and the demand fluctuation was reduced by about 40%.
Considering the diversity of configurations and demands in different IESs, it is possible in the future that adjacent IESs can exchange energy and jointly optimize their operation strategy to increase operational flexibility, reduce waste of surplus energy, and improve comprehensive utilization energy efficiency. The concept of energy network is therefore imposed, which aims to reduce the independence of a single IES and relieve the upper energy system burden [37,38]. With the communication technology advancement, it is possible for an IES to collaborate with neighboring IESs and complement spatio-temporal difference between supply and demand in a global win-win way. The upper retail energy price as well as operation cost may also be decreased by the total demand reduction in the region. Under this background, these IESs will have potential collaboration, and the above autonomous non-cooperative game theory is not well applicable. In consideration of the future energy network transaction mode, cooperative game theory has been introduced in some work recently. The cooperative theory stresses the maximization of the collective interest of coalition and that any individual in the coalition is increased or at least not damaged compaired with in non-cooperative mechanism. To this end, the distribution mechanism of reduced cost is particularly vital.
The Shapley value is a frequently used method for the profit allocation in the cooperative games [151]. It is based on the participants' marginal contribution to the alliance in the cooperation and can allocate profit in a fair way. The marginal contribution of player a i to the cooperative alliance S is the cost or benefit of the alliance increased by participants i joining the cooperative alliance. It can be calculated with Equation (14). The alliance A contains n players, and denoted as A = a 1 , · · · , a n . w(S i ) is the probability of forming a cooperative alliance S containing the subject a i . v(S i ) denotes the cooperative surplus of the alliance S with the participant a i . S i − a i denotes the remaining game alliance in S after removing a i , and the v(S i − a i ) is the cooperative surplus of the alliance without the participation of a i . The Shapley value ϕ i (v) represents the benefit allocating to a i for the alliance including n participants. With that, the surplus benefit can be distributed fairly according to the contribution of a certain participant in a quantitative way.
Ref. [152] adopted a cooperative game model and Shapley value to deal with multiple IESs cooperation optimal operation. The results showed that the pollutants emission and peak load were reduced by 45.59% and 27.06%, while all IESs had profit increment, and the total profit was increased by 53.54%. Similarly, Ref. [153] studied the optimal operation strategy of distributed PV systems with an externality-corrected mathematical model based on Shapley value and gave a benefit analysis. However, the Shapley value based allocation mechanism will be computationally complex and time-consuming when there are a number of participants. Some simplified methods are therefore introduced in view of that. Ref. [154] employed the bilateral Shapley value (BSV) model in the operation optimization of multiple prosumers integrated in a grid-connected IES. BSV solved combinatorial explosion via bilateral-marginalization of coalition, and the results showed that BSV could reduce calculation numbers effectively compared with the Shapley value model. Ref. [155] proposed a simplified profit allocation method based on the contribution of each prosumer's participation in additional profit compared with non-cooperative mode from the upper management perspective. Compared to the non-cooperative mode, the prosumers' costs were reduced by 4%. These simplified methods are effective in the research cases, yet their generality is not sufficient, and the mainstream is still the conventional Shapley value methods.
Game theory is extensively implemented into the collaborative operation optimization of multiple neighboring sub-IESs to handle the potential interconnections in both energy and market price. It is verified that game theory-based methods are effective in the coordination optimization problem and provide a promising way for sub-IES operators to realize the individual maximized interest considering others' effect.
It is noted that the market environment plays an essential role in the operational decisions. Energy companies can determine different prices for different periods to motivate consumers to adjust their consumption. Normally, the market-driven idea to coordinate resources is mainly employed in the power system. By now, the pricing schemes for electricity trading mainly focus on two aspects: dynamic (time-varying) pricing and static (pre-determined) pricing. The conventional static pricing scheme can be subdivided into the step wise pricing [156], time-of-use pricing [157], and critical peak pricing [158]. Today, the dynamic pricing becomes a hotpot, since it can effectively enhance participation motivation of consumers with market driven. It is proven that reasonable pricing mechanism can help in easing the strain on the electricity networks during peak periods and reducing the cost [159].
Currently, the electricity market, natural gas market, and heat energy market are operating separately. For heat and natural gas, they are conventionally considered as a natural monopoly product whose price is regulated by the government instead of load variation. That is, they are commonly fixed. However, some studies show that the marketdriven pricing policy is also beneficial for the gas market and heat market in peak shaving and cost reduction [160]. Moreover, due to the tight coupling of multiple energies in IES through various energy conversion devices, the original separate markets are not beneficial to realize optimal allocation of energy market sources and the greatest social benefits. For instance, when the electricity price is high in the wholesale market while the gas price is low, the IES operator can choose to generate power through the internal gas-driven generation units to satisfy users' electricity demand. In this case, the cost is reduced, and the peak burden from upper power system is released as well. Driven by that, the multi-energy trading market has recently drawn some researcher's attention.
For the gas-electricity market, Ref. [161] proposed a TOU gas pricing strategy for an IES. Gas storage, gas-fired absorption chiller, and gas-driven combustion driven engine in IES play an important role in the operation strategy adjustment according to price signals. The results showed that the production cost of the IES dropped by 20.54%, and meanwhile the utility went up to 5407.02 Yuan with TOU gas pricing compared with the fixed gas pricing. It indicated that the TOU gas pricing was more economical. In Ref. [162], an energy management model for IES based on dynamic energy pricing for electricity and natural gas was proposed.
Multiple trading of heat and electricity markets was also considered in some recent work, and the benefits are proven. Ref. [163] proposed an electricity-heat coordinated double auction retail energy market framework. By price signals guiding, the penetration ratio of wind power was more than 80%, and the energy expenditure could be reduced by 18.76% without sacrificing users' comfort. In Ref. [164], the conventional TOU electricity pricing was expanded to a multiple energies TOU pricing scheme (i.e., TOU heating pricing, TOU cooling pricing, and TOU electricity pricing) to explore the IES potential. In comparison with fixed pricing, the load curves were smoothed. The loads of cooling, heating, and power decreased by 0.719, 0.971, and 10.989 MW, respectively, during the peak period, and the difference between the peak and valley decreased by 0.121, 0.496, and 1.687 MW, respectively. Ref. [165] proposed a multi-energy market considering gas, heating, and electricity trading simultaneously and compared the results under four different market mechanisms (i.e., three independent energy markets, electricity-gas market, electricity-heat market, and the proposed multi-energy market). It demonstrated that the multi-energy market led to more profit, better end-user satisfaction, and flatter loads fluctuations.
In summary, the multiple-energy market has proven to be effective in mining potential of peak shaving, energy saving, and profit increment by price signals from the above work. As a medium that connects the upper energy systems and the lower end-users, IES plays a critical role through adjusting its operation strategy. However, the multiple-energy market concept is still in the preliminary stage of theoretical research. The various energy supply paths, strong couplings, and varying price signals dramatically complex the problem. There is a long way to go in transaction mechanism, pricing schemes for different energies, privacy protection, etc.

Challenges and Prospects
The volatile multi-energy loads largely increase the dependence of IES on the upper grid. The energy price signal is a promising method to motivate users to change their usage pattern. From a union perspective, it provides a potential to enhance IES independence. Since there are many participants who affect the strategy making and pursue a better individual benefit, such as the upper grid, the multiple IESs, and end-users, the key point lies in how to coordinate them and make an optimal strategy considering other's effect.
Game theory is introduced to solve interest coordination among multiple IESs. Noncooperative game theoretic methods are employed with SLMF and MLMF model and have been well developed. Recently, cooperative game is introduced to account for potential cooperation among sub-IESs in a future energy network. Some challenges and prospects are summed up as follows.

•
Models based on cooperative theory are still in their preliminary stage and some promotions are required. A rational and fair profit allocation mechanism is vital in cooperative games. The existing work is mostly based on the Shapley value method. However, it is not applicable to a complex problem. Some efficient and simplified allocation mechanism should be researched.

•
Market operation mechanism has a critical impact on the interconnection among multiple IESs through price scheming and trading mode. In present, the trading rules of energy trading among the internal sub-IESs are not formed. Trading mechanism in both internal energy trading and upper transactions with networks needs further study. The dynamic pricing scheme is a tendency to inspire more participation, and hence novel approaches such as the real-time game model may be in need to better satisfy the practical demand. • Except for electricity, the markets for heat and gas, which are conventionally regulated by the government instead of load variation, have been proven to be effective in mining potential of peak shaving, energy saving, and profit increment by price signals. However, they are still in the preliminary stage of theoretical research. More work is needed, especially in the transaction mechanism, pricing schemes for different energies, and privacy protection.

•
The above studies concentrate on the mutual influence of IESs on energy retail prices and assume the supply from upper energy networks are always sufficient. To obtain more practical and accurate results, the supply limits from upper input energy systems should be considered.

•
The users' characteristics can be comprehensively depicted. The existing work mostly focuses on the interest coordination among a group of IESs with identical configurations and load curves, while IESs group with various types is more in accord with practice. Diversity of sub-IESs may be considered in future research. Additionally, the preference and uncertainty of users and their impact on the final optimal IESs coordination strategy are always disregarded.

Conclusions
Integrated energy systems combine multiple energies to improve comprehensive efficiency and present a potential to face the energy crisis. However, due to volatile energy demands and restricted adjustment ability internal strong couplings, most of the existing integrated energy systems are dependent on the upper grid and cause a large pressure. The frequent interactions between integrated energy systems and the upper grid cause greater burden, since the renewable resources have occupied a larger and larger proportion in the energy system and at last may become dominant to reduce the fossil utilization. These weather-dependent generations are less controllable and flexible, and the system security and stability are hence more easily influenced by load fluctuations. To have better adaptability and relieve the grid burden, it requires the integrated energy system to increase independence and reconcile benefit.
This paper summarizes the existing work, which contributes to this aim from both directions of individual efforts and union efforts. More concretely, this study has systematically reviewed three aspects including the general formulation of integrated energy system optimal operation (individual benefit pursuit), individual efforts (i.e., handling methods with high proportional renewable energies to reduce interactions and exploiting the storage potential of multiple timescale characteristics of different energy flows to improve flexibility), and union efforts (coordination among multiple participants to manage load profile and exchange resources to realize energy complementary at the IES level).
(1) In terms of uncertainty handling, uncertainty modeling methods and multiple timescale dispatching strategy are utilized to jointly eliminate the adverse impact. At present, the data-driven uncertainty methods are emerged and fast developed to solve the problem of probability density distribution deficiency and the over conservativeness in conventional methods. Deep theoretical study is needed on the set of some critical parameters. For time-scale strategy, this is the optimum dispatch period selection for integrated energy system optimal operation with considerations of distinguishing dynamic characteristics of energy flows. (2) In terms of hybrid energy flows characteristics, gas-electricity integrated energy systems and heat-electricity integrated energy systems, as well as heat-gas-electricity integrated energy systems are reviewed. The main work focuses on the depiction of networks to reflect slow response characteristic of heat and gas. For both, the steady-state model is easy to implement and can be utilized when the accuracy is not concerned. Dynamic models will also be encouraged, as if it can be formulated in a tractable form. Some simplified or equivalent models seem to be efficient, which strike a balance between the tractability and accuracy, while the general form and verification need to be further explored. (3) In terms of multiple participants' efforts, game theory is adopted for the interest coordination problem. Cooperative optimization operation of different integrated energy systems is a trend for the urgent need under an energy network scenario in future. A rational and fair allocation mechanism is critical, and more work is in need. Additionally, the market operation mechanism plays an important role in the process. The dynamic pricing scheme is promising to encourage more flexible sources participants. More work can be done under this framework, especially with a multiple energies trading market.
Overall, the integration and interconnection of multiple energies in an integrated energy system greatly complicate the operation optimization problem, especially under a diversified complex trading market environment. This paper firstly aims at the uniqueness of integrated energy systems, concluding the effort directions considering independence and benefit to cater for a renewable energy dominated energy structure. This paper makes a thorough survey of each aspect to draw on our sights in possible solutions to existing