Operating Strategy for Local-Area Energy Systems Integration Considering Uncertainty of Supply-Side and Demand-Side under Conditional Value-At-Risk Assessment

To alleviate environmental pollution and improve the energy usage efficiency of terminals, energy systems integration (ESI) has become an important paradigm in the energy structure evolution. Power, gas and heat systems are becoming tightly interlinked with each other in ESI. The dispatching strategy of local-area ESI has significant impact on its operation. In this paper, a local-area ESI operational scheduling model based on conditional value-at-risk (CVaR) is proposed to minimize expected operational cost, which considers the uncertainty of energy supply-side and demand-side as well as multi-energy network constraints, including electrical network, thermal network and gas network. The risk cost is analyzed comprehensively under the condition of underor overestimated cost. On this basis, a hybrid method combining particle swarm optimization with interior point algorithm is executed to compute the optimal solutions of two-stage multi-period mixed-integer convex model. Finally, a case study is performed on ESI to demonstrate the effectiveness of the proposed method.


Introduction
The energy industry has gone through great changes as a result of the rapid development of technology and economy since twenty-first century.In recent years, the development of sustainable, affordable, and clean sources of energy is considered a prerequisite for today's economic strength and will benefit tomorrow's society.Under the impetus of competition in the energy industry, the unbundling of the electricity sector has introduced new technologies for the generation and the delivery of electricity, which signify less pollutant, highly efficient, and less costly ways of supplying the electricity.The large-scale integration of variable renewable units in power systems would require a fast response generating capacity.In this case, energy systems integration (ESI) as an important form of a new energy system is proposed; it covers electric power system, heat system and natural gas system with the purpose of coordinating the whole system in the processes of generation, transmission and consumption, which is considered to be the core technology of the "third industrial revolution" [1].A recent survey of the references shows that considering their characteristics range of application, the different mainstream ESI approaches adopted by researchers can be classified into two categories [2,3]: wide-area ESI and local-area ESI.Local-area ESI has many forms and includes city ESI, community ESI, industrial park ESI, etc.The paper focuses on the local-area ESI which lays solid foundations and provides valid support for entire energy systems; the local-area ESI can be seen as an upgraded level of a microgrid.The relationship between wide-area ESI and local-area ESI is shown in Figure 1.
forms and includes city ESI, community ESI, industrial park ESI, etc.The paper focuses on the local-area ESI which lays solid foundations and provides valid support for entire energy systems; the local-area ESI can be seen as an upgraded level of a microgrid.The relationship between wide-area ESI and local-area ESI is shown in Figure 1.There are a large number of references on ESI economic dispatch and many works have been conducted with regard to the decision making model.Parisio et al. [4] propose a robust optimization method for energy hub management which can be defined as an interface between generating units and the energy consumers, bounded uncertainties on energy hub parameters are taken into account and calculation methods are exploited to gain robust solutions which are feasible for all values.In [5], An integrated model of gas and electricity networks was developed to analyze the operational strategy using deterministic method, two stage stochastic method and multi stage stochastic method, which takes into account the uncertainty in wind power forecast and fuel availability of gas-fired generators.In [6], the coordination of constrained electricity and natural gas infrastructures is considered for firming the variability of wind energy in electric power systems, and the stochastic security-constrained unit commitment is applied for minimizing the expected operation cost in the day-ahead scheduling of power grid.In [7], an interval optimization based coordinated operating strategy for the gas-electricity integrated energy system is proposed considering demand response and wind power uncertainty.Zhang et al. [8] present a co-optimization planning model that considers the long-term interdependency of natural gas and electricity infrastructures, whose model is decomposed into a least-cost master investment problem for natural gas and electricity systems and interacts with two operation subproblems representing the feasibility (security) and the optimality (economic) of the proposed co-optimization.
The approaches mentioned above have the limitation that no mechanism was presented to deal with the gap between forecast values and real-time values for both energy supply and demand sides in ESI.Specifically, in energy supply side, high penetration distributed generations (DG) bring about more and more randomness and fuzziness, which is considerably difficult to make potentially costly operational decisions for energy management system (EMS).Meanwhile, in energy demand side, although several sources of loads have strong coupling relationships with each other that may There are a large number of references on ESI economic dispatch and many works have been conducted with regard to the decision making model.Parisio et al. [4] propose a robust optimization method for energy hub management which can be defined as an interface between generating units and the energy consumers, bounded uncertainties on energy hub parameters are taken into account and calculation methods are exploited to gain robust solutions which are feasible for all values.In [5], An integrated model of gas and electricity networks was developed to analyze the operational strategy using deterministic method, two stage stochastic method and multi stage stochastic method, which takes into account the uncertainty in wind power forecast and fuel availability of gas-fired generators.In [6], the coordination of constrained electricity and natural gas infrastructures is considered for firming the variability of wind energy in electric power systems, and the stochastic security-constrained unit commitment is applied for minimizing the expected operation cost in the day-ahead scheduling of power grid.In [7], an interval optimization based coordinated operating strategy for the gas-electricity integrated energy system is proposed considering demand response and wind power uncertainty.Zhang et al. [8] present a co-optimization planning model that considers the long-term interdependency of natural gas and electricity infrastructures, whose model is decomposed into a least-cost master investment problem for natural gas and electricity systems and interacts with two operation subproblems representing the feasibility (security) and the optimality (economic) of the proposed co-optimization.
The approaches mentioned above have the limitation that no mechanism was presented to deal with the gap between forecast values and real-time values for both energy supply and demand sides in ESI.Specifically, in energy supply side, high penetration distributed generations (DG) bring about more and more randomness and fuzziness, which is considerably difficult to make potentially costly operational decisions for energy management system (EMS).Meanwhile, in energy demand side, although several sources of loads have strong coupling relationships with each other that may improve the forecasting accuracy, uncertain energy consumption fashion in local-area still poses an enormous dilemma regarding the appropriate implementation of prediction algorithms.Therefore, as the risks of different energy systems are gradually aggregated, the interaction among various energy sources becomes more relevant, as illustrated in Figure 2. A number of energy service issues will be more complex to solve without assessing system risk, such as the economic scheduling of distributed energy system, scheduling of energy purchase, system security assessment, and planning for energy transactions and demand response.Thus, how to effectively quantify the risk of ESI is still an urgent problem to be solved.Value-at-risk (VaR) and conditional value-at-risk (CVaR), as well-known risk measures, have been widely used in various energy management problems for different entities in the electricity market such as retailers [9], producers [10][11][12], distribution companies [13], and in generation operation planning [14], optimal power flow [15] as well as coordinated energy trading problems [16,17].improve the forecasting accuracy, uncertain energy consumption fashion in local-area still poses an enormous dilemma regarding the appropriate implementation of prediction algorithms.Therefore, as the risks of different energy systems are gradually aggregated, the interaction among various energy sources becomes more relevant, as illustrated in Figure 2. A number of energy service issues will be more complex to solve without assessing system risk, such as the economic scheduling of distributed energy system, scheduling of energy purchase, system security assessment, and planning for energy transactions and demand response.Thus, how to effectively quantify the risk of ESI is still an urgent problem to be solved.Value-at-risk (VaR) and conditional value-at-risk (CVaR), as well-known risk measures, have been widely used in various energy management problems for different entities in the electricity market such as retailers [9], producers [10][11][12], distribution companies [13], and in generation operation planning [14], optimal power flow [15] as well as coordinated energy trading problems [16,17].The main contributions of this paper can be summarized as follows.
(1) To the best of our knowledge, the CVaR as risk evaluation approach has rarely been applied in local-area ESI operating strategy research.In this paper, risk assessment of local-area ESI is proposed to quantify the impact occurred by uncertainty of energy supply-side and demand-side.(2) In the modeling process, risk cost, which includes the fuel cost, maintenance cost and facility adjustment cost, is especially analyzed in the case of overestimation and underestimation as multiple sources of energy interact, and constraints of multi-energy systems network are fully taken into account.(3) By considering the advantage of the algorithm and the complexity of the ESI model, a bi-level optimization algorithm which combines particle swarm optimization (PSO) approach with interior point (IP) method is used to solve the stochastic multi-period mixed-integer model for dispatching issue.
The content of this paper is organized as follows.Section 2 presents the fundamental mathematical model related to energy supply equipment and uncertainty analysis on energy supply and demand side.Section 3 introduces the CVaR based economic dispatch model for integrated energy systems.The results of our numerical experiment and sensitivity analysis are presented in Section 4. Finally, the conclusion and plan for future work avenues are presented in Section 5.The main contributions of this paper can be summarized as follows.
(1) To the best of our knowledge, the CVaR as risk evaluation approach has rarely been applied in local-area ESI operating strategy research.In this paper, risk assessment of local-area ESI is proposed to quantify the impact occurred by uncertainty of energy supply-side and demand-side.(2) In the modeling process, risk cost, which includes the fuel cost, maintenance cost and facility adjustment cost, is especially analyzed in the case of overestimation and underestimation as multiple sources of energy interact, and constraints of multi-energy systems network are fully taken into account.(3) By considering the advantage of the algorithm and the complexity of the ESI model, a bi-level optimization algorithm which combines particle swarm optimization (PSO) approach with interior point (IP) method is used to solve the stochastic multi-period mixed-integer model for dispatching issue.
The content of this paper is organized as follows.Section 2 presents the fundamental mathematical model related to energy supply equipment and uncertainty analysis on energy supply and demand side.Section 3 introduces the CVaR based economic dispatch model for integrated energy systems.The results of our numerical experiment and sensitivity analysis are presented in Section 4. Finally, the conclusion and plan for future work avenues are presented in Section 5. Micro-turbine is a type of small gas turbine, whose rated power generally ranges from 25 kW to 300 kW.The efficiency of micro-turbine, which supplies energy source in the form of heat and electricity, can be up to 70%.Considering that the Capstone C65 micro is widely used in the practical application [18], power generation efficiency η MT can be approximated by Equation (1):

Mathematical Model Formulation
where P CHP denotes the electrical power of micro turbine; when the fluctuation of P CHP is limited, η MT can be regarded as a fixed constant.
LHV ng represents natural gas calorific value.The amount of nature gas consumed by CHP unit can be described as follow: C ng represents the price of nature gas.The fuel cost of CHP C CHP can be described as: η l denotes the heat loss rate, η rec denotes waste heat recovery efficiency, and COP h represents the heating coefficient of bromine cooler.The heating power of CHP can be noted as:

Fuel Cell Model
Fuel cell is a device that converts the chemical energy from a fuel into electricity through a chemical reaction of positively charged hydrogen ions with oxygen or another oxidizing agent [19].The amount of nature gas consumed by FC per unit time can be calculated as: where P FC presents generation power and η FC presents generation efficiency.The power efficiency of a fuel cell η FC and fuel cost of FC C FC can be described as follows: 2.1.3.Photovoltaic Power Generation Model Photovoltaics (PV) are a type of equipment that convert light into electricity using semiconducting materials that exhibit the photovoltaic effect.Considering that the output power of PV is closely related to temperature, radiance, pressure humidity, etc., the output of the photovoltaic cell can be expressed by Equation ( 8) where P PV (t) is the output of PV.G STC , T STC and P STC are radiance, photovoltaic cell temperature and maximum output power, respectively.G(t) and T c (t) are radiance and photovoltaic cell temperature at t-th period, respectively.k is the power temperature coefficient.

Wind Power Generation Model
A wind turbine is a device that converts the wind's kinetic energy into electrical power.Specifically, cut-in wind speed v ci rated wind speed v r and cut-out wind speed v co as three typical wind speed parameters and can comprehensively describe wind resource utilized by facility.The relationship between the WT output power and the wind speed are described in Equation ( 9): where P WT is actual output power of WT, P r is rated power, and k 1 and k 2 are fitting parameters.

Ground Source Heat Pump Model
Ground source heat pump model GSHP is a central heating or cooling system that transfers heat from the ground.It uses the earth as a heat source (in the winter) or a heat sink (in the summer) [20].The design takes advantage of the moderate temperatures in the ground to boost efficiency and reduce the operational cost of heating and cooling systems.In this paper, we employ a quadratic function as mathematical model to simplify the relationship of input and output.
where a, b and c are the coefficients of the quadratic function.

Gas Boiler Model
Gas boiler (GB) is the most commonly used thermal facility; it directly burns natural gas to produce heat [21].The amount of nature gas GB consumed per unit time can be calculated as follow: where Q GB is heating power and is heating coefficient.Fuel cost of GB C GB can be calculated as follow:

Uncertainty of Energy Supply Side
Given that renewable energy sources including PV and WT usually possess the characteristic of randomness and uncertainty in the ESI energy supply side, the Weibull distribution and the β distribution are commonly used to describe the uncertainty of wind speed and radiance for long-term power system planning and operation scheduling issues in conventional fashion, which ignores coupling between adjacent periods.These strategies might not be suitable for ESI economic scheduling in short-term or very short-term range to simulate the actual operation scenario.Thus, we employ short-term forecasts of wind speed and radiance to predict WT and PV output the next day.
The forecasting error of the wind speed ∆v can be obtained from the deviation between the wind speed predicted value _ v and the actual value v, which can be calculated as ∆v = v − _ v.Some studies have shown that the prediction error of wind speed is subject to normal distribution [22] whose mean and standard deviation are 0 and σ v , respectively.Thus, the probability of the actual wind speed is: The distribution function is: where Φ is standard normal distribution function.
According to Equation ( 9), for the reason that the relationship between WT output and wind speed is nonlinear, the output of the WT is not subject to the normal distribution.The output power probability distributions of WT can be calculated by Equations ( 9) and ( 14): The probability density function of WT output in interval [v ci , v r ] is: As per the analysis above, we can obtain the distribution of the WT output for the next 24 h after estimating the standard deviation σ v .
Similarly, according to Reference [23], PV output can be computed in a similar way.

Uncertainty of Energy Demand Side
A large body of references suggests that the probability density distribution of electrical and heat load approximately obeys normal distribution.Therefore, the normal distribution model with boundary is used to describe the uncertainty of electrical and heat loads [24,25] because the maximum and minimum points always appear during on-peak and off-peak periods.
Suppose that the random variable µ is subject to the normal distribution in the interval [a, b], the probability density function of µ can be described as: where µ 0 is the mean of random variables µ, and σ 0 is the standard deviation.Suppose that the predicted value of electrical load is P 0 , the actual value is P d .The minimum and maximum value are P d,min and P d,max , respectively.In addition, P d is subject to the normal distribution with standard deviation, σ d1 , thus the probability density function of the load actual value P d in the interval [P d,min , P d,max ] can be described as: Similarly, the probability density function f d (Φ d ) of the heat load actual value Φ d can also be obtained by method mentioned above.

Reduction of Uncertainty Scenarios
In this paper, we employ the Monte Carlo based Latin Hypercube Sampling (LHS) technique to generate scenarios by means of randomly sampling and comparing the probability distance.More details are provided in Reference [26].

Risk Modelling via CVaR
Conditional value at risk (CVaR) assessment is used to quantify the risk potential beyond VaR and represents an appropriate approach to integrate the inherent risk management problem in the commitment and dispatch procedures, which is an improved risk analysis method proposed by Rockafeller and Uryasev [27,28].
By definition, with respect to a specified probability level β, VaR of a portfolio is the lowest dispatch cost such that, with probability β, the dispatch cost will not exceed a considered amount, whereas the CVaR is the conditional expectation of dispatch cost above that amount.
Assuming that the probability density function of the random variable ξ is p(ξ), the distribution function of the loss function is: where f (X, ξ) is the loss function of portfolio X, and ξ represents a continuous random variable that may affect the loss function.
For a given confidence level β, the VaR of Portfolio X can be described as: The continuous CVaR can be illustrated as: Accordingly, the discrete CVaR can be calculated as: Here where n is number of period, and Pr s is the probability of occurrence for the s-th scenario.
The corresponding optimization formulation can be obtained according to Reference [22]: where f (X, z) is CVaR β , and z is VaR β .

Risk Cost of ESI
Since the predicted values of renewable energy source output and various types of loads may not able to track fluctuation of actual operational values comprehensively, which will lead to the unbalance of electrical flow, thermal flow, and natural gas flow, some measures must be taken to maintain the equilibrium state of ESI, which may bring about extra charge such as extra fuel cost, maintenance cost and additional adjustment expenses, besides the cost of purchasing natural gas and electricity.Thus, these sorts of extra charge can be defined as the risk cost of ESI economic scheduling, which can be divided into two categories: overestimated risk cost and underestimated risk cost.

Overestimated Cost
On condition that the forecasting value of PV and WT output are higher than the actual value or the forecasting value of electrical and thermal load are lower than those of actual value, the extra expense incurred by action, such as increasing output of controllable generation or cutting part of electrical and thermal load, can be defined as overestimated cost, which includes fuel cost for the controllable equipment, maintenance cost, adjustment cost and compensation for interruption of load.The overestimated cost can be described separately in terms of electricity and heat.

Overestimated Cost for Electricity
Overestimated cost for electrical subsystem can be described as follows: Here, superscript E is electrical power; Subscript t and i represent the t-th period and the i-th electrical power generation equipment, respectively; C E h,t is the overestimated cost; C E i,t is the fuel cost function, whose calculation detail has mentioned in Section 2; P i,t is planning electrical power output; N E is the number of controllable power generation facility; δ E i,t is the adjustment electrical power; r E u,i is the ramping up rate; K E i is the maintenance cost; c E T,i is the adjustment cost for per unit power; δ E h,t is the overestimated electrical power; c E loss is the compensation cost for interruptible load; and ∆P loss,t is the loss power for interruptible load.The second part of Equation (25) indicates that, even when all the controllable electrical power supply equipment operates at the maximum ramping up rate, the equilibrium point of ESI cannot be guaranteed, which will interrupt part of the load.The overestimated power can be described as follow: where P L,t represents electrical load forecasting power; P line,t is the exchange electrical power between the ESI and utility grid; P PV,i,t and P WT,j,t are PV and WT forecasting power output, respectively; δ E L,t , δ PV,i,t and δ WT,i,t are forecasting error of electrical load, PV and WT, respectively; and N PV and N WT are the number of PV and WT facilities, respectively.
The interrupt electrical load can be calculated as follow: Overestimated Cost for Heat The overestimated cost for heat is demonstrated as follows: where superscript H is heat power; Subscript t and i represent the t-th period and i-th thermal power generation equipment, respectively; C H h,t is the overestimated cost to guarantee thermal power balance; C H i,t is the fuel cost function and its detail can be found in Section 2; Φ i,t is the planning thermal power output; N H is the number of thermal equipment; δ H i,t is the adjustment power; r H u,i is the climbing rate; K H i is the maintenance cost for per unit thermal power; c H T,i is the adjustment cost for per unit thermal power; δ H h,t is the overestimated heat; c H loss is the compensation cost for interrupting thermal load; and ∆Φ loss,t is the interruption amount of the thermal load.The second part of Equation (28) indicates that, even when all the controllable thermal power supply equipment operates at the maximum ramping up rate, the equilibrium point of ESI cannot be guaranteed, which will interrupt part of the thermal load.The overestimated thermal power is noted as follow: where Φ L,t is the forecasting thermal load power, and δ H L,t is the prediction error of the thermal load.The interrupt thermal load can be calculated as follow:

Underestimated Cost
On condition that the forecasting value of PV and WT output are lower than the actual value or the forecasting value of electrical and thermal load are higher than those of actual value, some measures should be taken to keep power balance, such as controllable generation power reduction, and renewable energy curtailment.In this case, the related cost is defined as underestimate costs, which include the fuel cost for adjusting the controllable equipment, the maintenance cost, the punishment cost of renewable energy curtailment, and the adjustment cost.The underestimated cost can be described separately in terms of electricity and heat.

Underestimated Cost for Electricity
Underestimated cost for electrical subsystem can be described as follows: where C E d,t is the underestimated cost for electricity; r E d,i is the ramping down rate; δ E d,t is the underestimated electrical power; and c E pun is penalty cost for wasting electrical power.The second part of Equation (31) indicates that, even when all the controllable electrical power supply equipment operates at the maximum ramping down rate, the equilibrium point of ESI cannot be guaranteed, which will waste part of the electrical power.
The underestimated electrical power is noted as follow: The wasted electrical power is described as follow: Underestimated Cost for Heat Underestimated cost for thermal subsystem can be described as follows: where C H d,t is the underestimated cost; r H d,i is the climbing down rate; δ H h,t is the underestimated heat in the t period; and c E pun is penalty cost for wasting thermal power.The second part of Equation (34) indicates that, even when all the controllable thermal power supply equipment operates at the maximum ramping down rate, the equilibrium point of ESI cannot be guaranteed, which will waste part of the thermal power.
The underestimated thermal power is demonstrated as follow: The wasting heat is noted as follow:

Total Cost
The total risk cost of the ESI economic dispatching in the t-th period is the sum of overestimated cost and underestimated cost: where u E t and u H t are 0-1 variables; u E t = 1 and u H t = 1 indicate that overestimated cost of electrical power and thermal power supply will be incurred by the ESI operation, while u E t = 0 and u H t = 0 indicates that underestimated cost of electrical and thermal power supply will be incurred by the ESI operation.
Based on the analysis above, the total operating cost of the ESI in the s-th scenario for the entire T period can be calculated as follows where the first and second terms of the equation are the operation cost of power supply equipment and heating equipment in ESI in the case of the planned output, respectively; and M(P i,t ) = K E i P i,t and M(P i,t ) = K E i Φ i,t are maintenance cost for the corresponding device.The third term of Equation represents risk cost of ESI.The fourth item is the cost for electrical flow interaction between ESI and utility grid.P i,t , Φ i,t and P line,t without subscript variable s are the planned power output of controlled energy supply equipment.δ E i,t,s , δ H i,t,s and δ line,t,s with the subscript variable s are the actual adjustment output of the controllable energy supply equipment in the s-th scenario.

Objective Function
Total ESI operational cost based on CVaR is selected as the optimization objective to set up the economic dispatch model.Based on the theory of CVaR mentioned in Section 3.1, the formulation of the objective function can be presented as follow:

Constraint Conditions
The Upper and Lower Bound of Energy Supply Equipment Output where P i,min and P i,max are lower bound and upper bounds of output of the i-th electrical power supply equipment, respectively; Φ i,min and Φ i,max are lower bound and upper bounds of output of thermal power supply equipment, respectively; and U E i,t and U H i,t are 0-1 variables which represent the starting or stopping state of the thermal and electrical power supply equipment, respectively.

Power Balance Constraints
ESI operation should meet the requirement of balance on electrical power, thermal power and natural gas flow at any time: where ∆p t,s and ∆φ t,s are power loss of electrical network and thermal network in the t-th period, respectively; f source,t is the total amount of natural gas consumption of ESI in the t-th period; N e is the number of power supply device using natural gas as primary energy; f i,t is the amount of natural gas consumed by the i-th equipment in the t-th period; and δ g i,t is the variation of natural consumption due to the output adjustment of the energy supply equipment.

Tie-Line Power Constraints
The exchange power of tie-line between the ESI and utility grid should meet the upper and lower bounds: where P line,max and P line,min are the power upper and lower bounds for exchange power of tie-line between the ESI and utility grid.

Ramping Rate Constraints
The output variation of electrical power and thermal supply device needs to meet the ramping up rate and ramping down rate constraints: where r E u,i and r E d,i are ramping up rate and ramping down rate of the electrical power supply device, respectively; and r H u,i and r H d,i are ramping up rate and ramping down rate of the thermal supply device, respectively.

Multi-Energy Flow Constraints
The state variables and control variables of each line, pipeline and node which belong to the electrical power subsystem, thermal subsystem and natural gas subsystem should meet the constraints of multi-energy flow.Details can be found in Reference [29].

System State Variables Constraints
To ensure the ESI operates in the safe zone, state variables x = [θ, V, m, T s,load , T r,load , Π] T in the power subsystem, thermal subsystem and natural gas subsystem should be in the specified range, where x t,s is the state vector of the system in the s-th scenarios at the t-th (t = 1, 2, • • • , T) period; and x min and x max are the column vectors for each state lower and upper bounds.Details can be found in Reference [29].

Solving Algorithm
The CVaR based ESI economic scheduling model is concerned to a hybrid optimization design problem with a large number of disperse and continuous variables.This nonlinear mixed integer programming problem requires considerable calculations and storage resources, which is complex to solve without comprehensively analyzing the relationship between the variables and the model.Therefore, given that characteristic of different algorithms, we deploy a bi-level optimization model [28] that combines particle swarm optimization (PSO) approach [30] with interior point (IP) method [31].
In this paper, the particle swarm optimization algorithm is used to solve the sub problem of upper level, and the interior point method is used to solve the sub problems of lower level, the optimal solution can be obtained by the constant iteration of the upper level and the lower level until the variation of error is less than the mismatch tolerance.Furthermore, for the system control variables and the state variables in the optimization model, when the control variables are fixed in iterative process, the state variables can be determined by multi-energy flow calculation algorithm mentioned in Reference [31], which analyzes the flow of various energy systems jointly based on extended Newton-Raphson algorithm.This method can significantly reduce the number of unknown variables and simplify model complexity.
X represents the variables that are independent of the uncertain scenarios, including the planning output of each device and z, X s represents the variables related to the uncertain scenarios, which is the adjusted output of controllable energy supply equipment in each scenario.The optimization model based on CVaR is shown as Equation (47).min F(X, X s ) = e(X) where e(X) is a function of X; h(X) = 0 and g(X) ≤ 0 are the equality and inequality constraints regarding planning output mentioned in Section 3.3, respectively; and h s (X, X s ) = 0 and g s (X, X s ) ≤ 0 are the equality and inequality constraints regarding adjusted output of the equipment mentioned in Section 3.3, respectively.
In the simplified form of optimization model, the value of X s only has impact on the corresponding f s (X, X s ) solely, but X has impacts on f s (X, X s ) and F(X, X s ).In addition, we can see that f s (X, X s ) − e(X) and f s (X, X s ) have positive correlation, both Pr s and 1/(1 − β) are positive numbers.Thus, the objective function F(X, X s ) and f s (X, X s ) are positively correlated.Therefore, to get the minimum value of the objective function F(X, X s ), the minimum value of function f s (X, X s ) can be obtained by adjusting X s , when X is fixed.
According to the above analysis, the bi-level optimization model is described as follows: Upper Level Optimization Model The multi-energy flow in the upper level optimization model is calculated by the extended Newton-Raphson algorithm [31].

Lower Level Optimization Model
The lower level optimization model contains N s sub-problems, and the objective function and constraints of the s-th (s = 1, 2,• • • , N s ) sub-problem are shown as follows: 40), (41), ( 43), ( 44), ( 45), (45) Besides the constraints mentioned above, the inequalities and equality constraints in Reference [31] are also included in constraints.
The flowchart of hybrid algorithm is illustrated in Figure 3.

Case Study
The suitability of the proposed approach for conducting CVaR based economic dispatch studies was tested on the system network shown in Figure 4. EBi, HBi and GBi served as the nodes of the power, thermal, and gas systems, respectively.A gas compressor was installed in the gas subsystem to ensure pressure and gas supply.The electrical power subsystem is connected to the utility grid via the node EB12.PV with rated power of 160 kWp and 240 kWp is installed at nodes EB3 and EB8, and WT with rated power 1000 kW is installed at node EB4.CHP unit works in the following the electric load (FEL) model.The forecasting values of PV, WT, electrical load and heat load in each period are demonstrated in Figures 5 and 6.The parameters of power supply equipment in ESI are shown in Appendix A. The power factor of each electric load node is 0.9.Assuming that the forecasting errors of PV and WT are subject to the normal distribution with variance of σ μ =10% and σ μ =15% , respectively, the forecasting errors of the electrical load and thermal load are both subject to the normal distribution with variance of =5% σ μ , respectively.

Case Study
The suitability of the proposed approach for conducting CVaR based economic dispatch studies was tested on the system network shown in Figure 4. EBi, HBi and GBi served as the nodes of the power, thermal, and gas systems, respectively.A gas compressor was installed in the gas subsystem to ensure pressure and gas supply.The electrical power subsystem is connected to the utility grid via the node EB12.PV with rated power of 160 kWp and 240 kWp is installed at nodes EB3 and EB8, and WT with rated power 1000 kW is installed at node EB4.CHP unit works in the following the electric load (FEL) model.The forecasting values of PV, WT, electrical load and heat load in each period are demonstrated in Figures 5 and 6.The parameters of power supply equipment in ESI are shown in Appendix A. The power factor of each electric load node is 0.9.Assuming that the forecasting errors of PV and WT are subject to the normal distribution with variance of σ = 10%µ and σ = 15%µ, respectively, the forecasting errors of the electrical load and thermal load are both subject to the normal distribution with variance of σ = 5%µ, respectively.The time of use share price (TOU) for the electrical power exchange between the ESI and the external utility grid is demonstrated in Table 1.The price of natural gas is 2.51 Yuan/m 3 .The adjustment cost of the electrical power supply facility is 0.52 Yuan/kW.The adjustment cost of the thermal power supply facility is 0.42 Yuan/kW.The compensation cost for interrupting electrical loads and thermal loads are 0.612 Yuan/kWh and 0.452 Yuan/kWh, respectively.The penalty cost for wasting electrical power and thermal power are 0.32 Yuan/kWh and 0.22 Yuan/kWh, respectively.For line constants of the ESI, refer to Appendix B.  The time of use share price (TOU) for the electrical power exchange between the ESI and the external utility grid is demonstrated in Table 1.The price of natural gas is 2.51 Yuan/m 3 .The adjustment cost of the electrical power supply facility is 0.52 Yuan/kW.The adjustment cost of the thermal power supply facility is 0.42 Yuan/kW.The compensation cost for interrupting electrical loads and thermal loads are 0.612 Yuan/kWh and 0.452 Yuan/kWh, respectively.The penalty cost for wasting electrical power and thermal power are 0.32 Yuan/kWh and 0.22 Yuan/kWh, respectively.For line constants of the ESI, refer to Appendix B.  The time of use share price (TOU) for the electrical power exchange between the ESI and the external utility grid is demonstrated in Table 1.The price of natural gas is 2.51 Yuan/m 3 .The adjustment cost of the electrical power supply facility is 0.52 Yuan/kW.The adjustment cost of the thermal power supply facility is 0.42 Yuan/kW.The compensation cost for interrupting electrical loads and thermal loads are 0.612 Yuan/kWh and 0.452 Yuan/kWh, respectively.The penalty cost for wasting electrical power and thermal power are 0.32 Yuan/kWh and 0.22 Yuan/kWh, respectively.For line constants of the ESI, refer to Appendix B.    The dispatching period is . Most of the experiments could be completed on machine with a Core i7 central processing unit and 8 GB of memory using the 2014a MATLAB environment.

Simulation Based on Different Confidence Level
To analyze the influence of confidence level parameter on simulation results, the operating cost of VaR and CVaR under different confidence levels (β = 0.80, 0.85, 0.90, and 0.95) are presented in Figure 7.The output of different energy supply facilities under different confidence levels is shown in Table 1 and the output of energy supply facilities under different confidence levels is shown in Table 2.The dispatching period is T = 24 h, time interval is ∆t = 1 h.Most of the experiments could be completed on machine with a Core i7 central processing unit and 8 GB of memory using the 2014a MATLAB environment.

Simulation Based on Different Confidence Level
To analyze the influence of confidence level parameter on simulation results, the operating cost of VaR and CVaR under different confidence levels (β = 0.80, 0.85, 0.90, and 0.95) are presented in Figure 7.The output of different energy supply facilities under different confidence levels is shown in Table 1 and the output of energy supply facilities under different confidence levels is shown in Table 2.  Conclusions can be drawn from the simulation results.As the confidence level β increases, VaR and CVaR improve simultaneously, which means the dispatching risk of the ESI is aggregated.When β is increased proportionally, the difference between the VaR and the CVaR is gradually diminished because forecasting error of supply side and demand side are subject to the normal distribution.The potential risk of the energy systems gradually moves to the tail of the normal distribution, which makes the expected value of the corresponding cost increase, while the growth rate slows down.On the other hand, the mathematical model of the energy supply facilities deployed in the risk based economic dispatch possesses nonlinear characteristics, and the output variation of energy supply equipment presents no obvious regularity.

The Influence of Multi-Energy Flow Constraints on Simulation Results
The multi-energy flow constraints, which include the line, branch or pipeline limits of each subsystem, have been taken into account.To further analyze the influence of multi-energy flow constraints on simulation results, in the selected Ns scenario, we track the dynamic variation of essential variables during the dispatching process.The maximum and minimum value of voltage, temperature, and pressure, which belong to different subsystems, under different confidence levels, are illustrated in Figures 8-10.Conclusions can be drawn from the simulation results.As the confidence level β increases, VaR and CVaR improve simultaneously, which means the dispatching risk of the ESI is aggregated.When β is increased proportionally, the difference between the VaR and the CVaR is gradually diminished because forecasting error of supply side and demand side are subject to the normal distribution.The potential risk of the energy systems gradually moves to the tail of the normal distribution, which makes the expected value of the corresponding cost increase, while the growth rate slows down.On the other hand, the mathematical model of the energy supply facilities deployed in the risk based economic dispatch possesses nonlinear characteristics, and the output variation of energy supply equipment presents no obvious regularity.

The Influence of Multi-Energy Flow Constraints on Simulation Results
The multi-energy flow constraints, which include the line, branch or pipeline limits of each subsystem, have been taken into account.To further analyze the influence of multi-energy flow constraints on simulation results, in the selected Ns scenario, we track the dynamic variation of essential variables during the dispatching process.The maximum and minimum value of voltage, temperature, and pressure, which belong to different subsystems, under different confidence levels, are illustrated in Figures 8-10.In Figures 8-10, we can see that the nodal voltage of electricity subsystem, the temperature of thermal subsystem and the nodal pressure of natural gas subsystem are within the allowable range no matter in which confidence level.Compared with the conventional dispatching model based on CVaR, which neglects constraints of multi-energy flow, the proposed approach enables that ESI operates in the safe zone considerably, which can effectively eliminate the violated limit phenomenon and provide a reference for operator to control the system state of each node.

The Influence of Equipment Adjustment Cost on Simulation Results
For the sake of testing the influence of adjustment on simulation results, we define T c 0 as the original adjustment cost of energy supply equipment for per unit power.When T c are set as 0.25 times, 0.5 times, 0.75 times, 1.25 times and 1.75 times T c 0 under the condition that β = 0.90, the variation of the operating cost expectation value E, VaR and CVaR is presented in Figure 11.
In Figure 11, we can see that the distribution of the total operating cost is affected by T c : as T c increases, the expected value of the total operating cost E varies very limited.On the contrary, both VaR and CVaR decrease sharply, which means the scheduling risk of energy systems decreases at this time.Given these simulation results, EMS can effectively make corresponding decisions to reduce the risk cost.In Figures 8-10, we can see that the nodal voltage of electricity subsystem, the temperature of thermal subsystem and the nodal pressure of natural gas subsystem are within the allowable range no matter in which confidence level.Compared with the conventional dispatching model based on CVaR, which neglects constraints of multi-energy flow, the proposed approach enables that ESI operates in the safe zone considerably, which can effectively eliminate the violated limit phenomenon and provide a reference for operator to control the system state of each node.

The Influence of Equipment Adjustment Cost on Simulation Results
For the sake of testing the influence of adjustment on simulation results, we define c T0 as the original adjustment cost of energy supply equipment for per unit power.When c T are set as 0.25 times, 0.5 times, 0.75 times, 1.25 times and 1.75 times c T0 under the condition that β = 0.90, the variation of the operating cost expectation value E, VaR and CVaR is presented in Figure 11.
In Figure 11, we can see that the distribution of the total operating cost is affected by c T : as c T increases, the expected value of the total operating cost E varies very limited.On the contrary, both VaR and CVaR decrease sharply, which means the scheduling risk of energy systems decreases at this time.Given these simulation results, EMS can effectively make corresponding decisions to reduce the risk cost.

Convergence Analysis
To further illustrate the convergence of hybrid optimization algorithm based on PSO and IP for solving the proposed model, the convergence curves under the confidence levels β = 0.80, 0.85, 0.90, and 0.95 are demonstrated in Figure 12.We can see that the four groups of curves can converge effectively after 200 iterations, and the optimal solutions are obtained.

Conclusions
A CVaR-based optimization algorithm for energy management problem of ESI in the presence of uncertainty on energy demand and supply side was proposed in the paper, which has fully considered the constraints of electrical network, thermal network and natural gas network.Risk cost, including overestimated risk cost and underestimated risk cost, has been clearly analyzed according to practical demand.A bi-level optimization method based on the fast particle swarm algorithm and interior point method is used to solve ESI model.

Convergence Analysis
To further illustrate the convergence of hybrid optimization algorithm based on PSO and IP for solving the proposed model, the convergence curves under the confidence levels β = 0.80, 0.85, 0.90, and 0.95 are demonstrated in Figure 12.We can see that the four groups of curves can converge effectively after iterations, and the optimal solutions are obtained.

Convergence Analysis
To further illustrate the convergence of hybrid optimization algorithm based on PSO and IP for solving the proposed model, the convergence curves under the confidence levels β = 0.80, 0.85, 0.90, and 0.95 are demonstrated in Figure 12.We can see that the four groups of curves can converge effectively after 200 iterations, and the optimal solutions are obtained.

Conclusions
A CVaR-based optimization algorithm for energy management problem of ESI in the presence of uncertainty on energy demand and supply side was proposed in the paper, which has fully considered the constraints of electrical network, thermal network and natural gas network.Risk cost, including overestimated risk cost and underestimated risk cost, has been clearly analyzed according to practical demand.A bi-level optimization method based on the fast particle swarm algorithm and interior point method is used to solve ESI model.

Conclusions
A CVaR-based optimization algorithm for energy management problem of ESI in the presence of uncertainty on energy demand and supply side was proposed in the paper, which has fully considered the constraints of electrical network, thermal network and natural gas network.Risk cost, including overestimated risk cost and underestimated risk cost, has been clearly analyzed according to practical demand.A bi-level optimization method based on the fast particle swarm algorithm and interior point method is used to solve ESI model.
Through a case study, we can see the CVaR method can well evaluate the risk level of the ESI dispatching problem, which effectively avoids economic losses caused by the uncertain factors.As the confidence level increases, the scheduling risk increases accordingly.Furthermore, multiple constraints considering electricity, heat and gas flow ensure the system operation safely, which has more significance on theory and engineering.The hybrid optimization algorithm we deployed has good convergence property.

Figure 1 .
Figure 1.The relationship between wide-area energy systems integration (ESI) and local-area ESI.

Figure 1 .
Figure 1.The relationship between wide-area energy systems integration (ESI) and local-area ESI.

Figure 2 .
Figure 2. Economic dispatching mode for local-area ESI.

Figure 2 .
Figure 2. Economic dispatching mode for local-area ESI.

Figure 3 .
Figure 3.The solving process based on particle swarm optimization (PSO) and interior point (IP) method.

Figure 4 .
Figure 4. Case of local-area energy systems integration.

Figure 5 .
Figure 5.The forecasting values of PV and WT output power.

Figure 4 .
Figure 4. Case of local-area energy systems integration.

Figure 4 .
Figure 4. Case of local-area energy systems integration.

Figure 5 .
Figure 5.The forecasting values of PV and WT output power.

Figure 5 .
Figure 5.The forecasting values of PV and WT output power.

Figure 6 .
Figure 6.The forecasting values of electrical and thermal load.

Figure 7 .
Figure 7.The results of VaR and CVaR under different confidence levels.

Figure 7 .
Figure 7.The results of VaR and CVaR under different confidence levels.

Sustainability 2017, 9 , 1655 19 of 24 Figure 8 .
Figure 8.The maximum and minimum voltage of each node in the power subsystem.

Figure 8 .
Figure 8.The maximum and minimum voltage of each node in the power subsystem.

Figure 8 .
Figure 8.The maximum and minimum voltage of each node in the power subsystem.

Figure 9 .
Figure 9.The maximum and the minimum temperature of each node in the thermal subsystem.

Figure 9 . 24 Figure 10 .
Figure 9.The maximum and the minimum temperature of each node in the thermal subsystem.Sustainability 2017, 9, 1655 20 of 24

Figure 10 .
Figure 10.The maximum and minimum pressure of each node in natural gas subsystem.

Figure 11 .
Figure 11.The influence of energy supply equipment unit power adjustment cost on simulation results.

Figure 12 .
Figure 12.The convergence of the objective functions under different confidence level.

Figure 11 .
Figure 11.The influence of energy supply equipment unit power adjustment cost on simulation results.

Sustainability 2017, 9 , 1655 21 of 24 Figure 11 .
Figure 11.The influence of energy supply equipment unit power adjustment cost on simulation results.

Figure 12 .
Figure 12.The convergence of the objective functions under different confidence level.

Figure 12 .
Figure 12.The convergence of the objective functions under different confidence level.

Table 2 .
Output of energy supply facilities under different confidence levels.Figure 6.The forecasting values of electrical and thermal load.

Table 2 .
Output of energy supply facilities under different confidence levels.

Table A4 .
Pipeline coefficients of the natural gas system.