Design and Evaluation of Operational Scheduling Approaches for HCNG Penetrated Integrated Energy System

This paper proposes and assesses three different control approaches for the hydrocarbon natural gas (HCNG) penetrated integrated energy system (IES). The three control approaches adopt mixed integer linear programing, conditional value at risk (CVaR), and robust optimization (RO), respectively, aiming to mitigate the renewable generation uncertainties. By comparing the performance and efficiency, the most appropriate control approach for the HCNG penetrated IES is identified. The numerical analysis is conducted to evaluate the three control approaches in different scenarios, where the uncertainty level of renewable energy (within the HCNG penetrated IES) varies. The numerical results show that the CVaR-based approach outperforms the other two approaches when renewable uncertainty is high (approximately 30%). In terms of the cost to satisfy the energy demand, the operational cost of the CVaR-based method is 8.29% lower than the RO one, while the RO-based approach has a better performance when the renewable uncertainty is medium (approximately 5%) and it is operational is 0.62% lower than that of the CVaR model. In both evaluation cases, mixed integer linear programing approach cannot meet the energy demand. This paper also compares the operational performance of the IES with and without HCNG. It is shown that the IES with HCNG can significantly improve the capability to accommodate renewable energy with low upgrading cost.

INDEX TERMS Conditional value at risk, HCNG, integrated energy system, mixed integer linear programing, robust optimization. The Integrated Energy System (IES) is considered as an emerging energy system concept that has great potential on VOLUME 7, 2019 breaking the barriers and bridging different types of energy systems. Especially for the industrial park or area with considerable amount of heating/cooling demand and renewable energy resources, the IES has been proved as an efficient system for mitigating surplus renewable energy and providing high quality heating and cooling energy. The intension to convert a conventional energy system to an IES can be frustrated in specific areas due to the limited transportation capability of natural gas pipeline network and the congestion in the electricity power network. A more common way to upgrade an existing energy system is either to adopt energy storage system or expand the capability and scale of the energy network. However, such method can be very expensive and time consuming. With the proliferation of renewable energy resources and clean vehicles, the energy demand of end users can become harder to predict, which can lead to greater uncertainty level in the energy system. To give more credit when justifying the benefits of converting an existing energy system to an IES rather than using common system upgrading methods, it is important to show that the IES can accommodate the uncertainties in renewable resources and end user loads without the need to expand the existing gas pipeline network and the electric network. Power-to-gas (P2G) technology is a conventional method to convert surplus electricity to gas fuel [1]. The converted gas fuel can be natural gas (through methanation) or hydrogen (through electrolysis). The P2G technology using methanation is commonly used for large-scale natural gas production [2] whilst the P2G technology using electrolysis is more usual in medium or small scale hydrogen production aiming to satisfy the local hydrogen demand [3]. The produced hydrogen can be directly used as the fuel for fuelcell vehicle. It can also be injected into the gas pipeline to produce Hydrocarbon Natural Gas (HCNG) [4]. Due to the limitations of the compressor, the hydrogen content of HCNG cannot be too large [5]. Researches have proved that the HCNG can be used as the fuel for combustion engines, gas turbines and home appliances if the proportion of hydrogen is below certain percentages [6]- [10]. Reference [6] conducted the experiment of using HCNG to power the CHP, and the CHP still had satisfactory performance when the hydrogen proportion reached 8%. Reference [9] evaluated the effect of hydrogen percentage (i.e. from 0% to 30%) within the HCNG on CHP, and it was discovered that the CHP heat recovery efficiency was improved as the percentage of hydrogen raised and at the same time the CHP system could still work without any failures. CHP is the core equipment within the IES, and it is considered as the major alternatives to traditional electricity generators. Using HCNG as the fuel for CHP can potentially reduce the natural gas consumption. In addition, the hydrogen is easy to produce and store, which makes it easier for the IES to accommodate more surplus energy generated by renewable generators, in both short and long-term period. Considering the global commitment to reduce CO 2 emission, the HCNG penetrated IES can be the development direction of IES, and the investigation of how a HCNG penetrated IES can be to effectively operated is valuable. The use of HCNG as a fuel in the engine will significantly reduce emissions of THC, CO, CO 2 and CH 4 , but will increase NO X emissions [11]- [14]. And NO X can be processed to meet emission standards [12], [14].

Symbol Quantity
A number of previous studies were carried out on the operational optimization of IES. Reference [15] proposed an interval optimization approach for gas-electricity IES taking into account the demand response scheme and renewable uncertainties. Wei Gu and Jun Wang designed an integrated energy system optimization method which combines the thermal inertia of the regional heat network and the thermal inertia of the building to improve the utilization of wind power [16]. The interval based optimization approach mitigated the effects of worst scenario on the overall optimization, compared to the robust control. In [17], an operational optimization approach for heating-electricity IES based on heating inertia was presented. It has proven that the increased heating storage capability of a heating network would increase the operation flexibility of a CHP, which is a valuable lesson for the HCNG penetrated IES optimization (where the stored object is hydrogen instead of heat). Dolatabadi and Mohammadi-Ivatloo designed a conditional-value-atrisk (CVaR) based optimization approach for scheduling the operation of energy hub (a system coupling among different energy networks, similar to IES) [18]. Such methodology has shown good performance on mitigating the risk of load forecasting volatilities and energy price fluctuation. Reference [19] proposed a risk-averse optimal operation strategy for multiple-energy carrier system (similar concept to IES), using CVaR method to quantify the risk associated with loads uncertainties. With regard to the operation scheduling problem of energy system, the risk-accounted optimization approaches (e.g. CVaR, Robust Optimization, Interval Optimization) have been demonstrated as efficient ways to mitigate the risks caused by the uncertainties in renewable generation, load and price [20]- [24].
To the best knowledge of authors, there are little research that has been conducted on the operational optimization of the HCNG penetrated IES, as well as the mathematical modeling of the CHP using HCNG as the fuel. In the foreseeable future, increasing number of clean vehicles, renewable generators and real-time pricing schemes will be deployed and it will be necessary to identify methods to mitigate such uncertainties within the energy system. An expectation of using the HCNG penetrated IES is that it can provide appropriate risk mitigation capabilities through the use of flexible injection and storing of hydrogen. To demonstrate the advantages HCNG penetrated IES.
In this paper, the detailed model of HCNG penetrated IES is developed and such system can produce and use hydrogen to absorb extra wind and PV energy whilst at the same time reduce the carbon emissions. In order to mitigate the influence of uncertainties in renewable generation and loads, three planning strategies are developed and tested in this paper. Furthermore, two hydrogen fuel cells are applied in the HCNG penetrated IES model and their impact on the planning results are also investigated. After that, an operation planning method is proposed which can optimize the electrical and thermal efficiency of the CHP that is related to the change of the hydrogen concentration in the gas mixture. Finally, to validate the performance of the proposed model and algorithm, two kinds of deviations from the predicted values are used for numerical analysis and comparison. The main contributions of this paper can be concluded as below.
(1) In this paper, the quantitative mathematical model of CHP considering the influence of hydrogen concentration change in fuel is given, as well as the models of other refrigeration equipment, heating equipment and hot water equipment in the HCNG penetrated IES. The operation constraints of these facilities and the interaction with the power grid are also considered. (2) The relationship between the change of load electrothermal ratio and the hydrogen proportion of CHP fuel is studied. In addition, the performance of the CVaR model, robust optimization (RO) model and mixed integer linear programing is compared. (3) The optimization results of CVaR, RO and mixed integer linear programing based approaches are comprehensively analyzed from the perspective of expectation of running cost, resistance of uncertainty and running time of the program. The remainder of this paper is organized as follows: section 2 describes the mathematical model of the HCNG penetrated IES whilst section 3 presents three optimization methods for day-ahead economic dispatch. Section 4 discusses the advantages and disadvantages of the three optimization schemes under different fluctuation conditions.

II. SYSTEM MODEL A. HCNG PENETRATED IES
This paper presents and compares three day-ahead control approaches for the hydrogen penetrated IES. Compared to the traditional IES, a hydrogen penetrated IES employs hydrogen generation (e.g. electrolyzer), hydrogen storage (e.g. hydrogen tank) and hydrogen consumption (e.g. fuel cell or HCNG). Figure 1 presents a typical hydrogen penetrated IES. Apart from the aforementioned hydrogen related components, such system also has the gas related components, smoke recycling components, electricity based components and cooling components, and all of them within the system are connected by different energy-carrier busbars.
Most of CHPs can take natural gas mixed with hydrogen of certain percentage without any upgrading and this paper VOLUME 7, 2019 focuses on the hydrogen related components within the IES. The generated hydrogen is mixed with methane and used as the fuel for CHP. For the IES involving renewable energy resources such as solar and wind turbines, the hydrogen related devices can help the system operator to store the surplus electricity by transforming electricity to hydrogen. The stored hydrogen can be injected into the gas pipeline, which can contribute to the electricity and heat generation of CHP when there is a shortage in electricity and heat. As hydrogen is a carbon-free fuel and can be stored for long term (seasonally or even yearly) in a tank, it is considered as one of the potential replacements for the future fuel.
Note that all types of gas mentioned and discussed in this paper are measured in volume under normal condition and the unit is Nm 3 .

B. CONSTRAINS OF ENERGY-CARRIER BUSBARS
The hydrogen penetrated IES incorporates multi-type energy busbars in one single system. The system should satisfy the constraints of different energy-carrier busbars, which are defined by equations (1) -(4).
Equation (1) indicates the balancing constraint of the electrical busbar and equation (2) shows the balancing constraints of the hot water busbar. The balancing constrain of the smoke busbar is presented in equation (3) and the balancing constraint of the air busbar is illustrated in equation (4).

C. CHP MODEL USING HCNG AND ITS ASSOCIATED CONSTRAINTS
The fuel used by CHP in the hydrogen penetrated IES is the HCNG which is a mixture of methane and hydrogen. The volume proportion of hydrogen ε in the gas mixture is constrained by the volume of methane or Methane Number (MN) [6]. The detailed constraints of the hydrogen within the HCNG are given as below.
Equation (7) is a nonlinear equation, which will dramatically increase the optimization complexity for the dayahead control of hydrogen penetrated IES. Thus, equations (7) and (8) are rewritten as equation (10). The optimized proportion of hydrogen within the gas mixture in various time intervals can be derived from equation (7).
To illustrate the influence of the hydrogen injected into the methane, the electricity generation power, hot water output power and smoke output power generated from HCNG can be calculated based on equations (11) -(13), respectively. Assume the smoke output power is 10% of the total power produced from HCNG. Note that Lower Heating Value (LHV) is used when calculating power generated from HCNG.
According to [6], when the volume proportion of hydrogen in the mixed gas changed, the electricity-heat efficiency of the CHP would correspondingly change. Based on the experimental results given in [6], the efficiency change of the CHP using mixed gas could be reflected by equations (14) - (16).
Since equations (11) -(13) are nonlinear equations, it will be sophisticated to optimize the scheduling of the CHP operation. To reduce the calculation complexity, the output power model of the CHP is simplified and the relationship between the output power values and the flow rate of natural gas and hydrogen is shown in linear equations (17) - (19).
For the CHP, the electric power constraints and power ramping constraints are presented in equations (20) - (22). The binary variable µ CHP is used to indicate the operating status of the CHP where 1 means operating whilst 0 denotes not operating.

D. OTHER DEVICES MODEL AND OPERATION CONSTRAINTS 1) CONSTRAINTS OF HYDROGEN STORAGE TANK
In the proposed IES system, the hydrogen is stored in high pressure tank(s) under normal temperature. More specifically, the hydrogen will first be generated from the electrolyzer and subsequently compressed and injected into the hydrogen storage tank. Due to the negligible magnitude, the electrical power consumed by the air pump (of the gas pipeline) and compressor will not be considered in the modeling. The stored hydrogen will be consumed by the CHP and fuel cell to generate energy in different formats (e.g. electricity and heat). Equations (23) -(26) show the relationship between the hydrogen volume in the storage tank and the flow rate in different devices as well as the associated constraints.

2) CONSTRAINTS OF ELECTROLYZER AND FUEL CELL
Even though the electrolyzer consumes electricity to generate hydrogen whilst the fuel cell uses hydrogen to generate electricity, they cannot operate simultaneously to maintain a continuous loop where electrolyzer takes electricity from the fuel cell and the fuel cell takes the hydrogen from the electrolyzer. This is because there will be loss in both processes. Therefore, two binary variables µ EL and µ HFC are used to indicate the device operating status: 1 as operating and 0 as not operating.
No heat energy could be retrieved from a conventional fuel cell during the electricity generation process. However, if a Solid Oxide Fuel Cell (SOFC) is used, the resulted heat from electricity generation could be captured and reused. Assume the heat generated from SOFC is used to boil the water and the resulted hot water power is 80% of electricity power counterpart of SOFC. The power relationship and constraints of the electrolyzer and fuel cell are presented by equations (27) to (33).

3) CONSTRAINTS OF ELECTRIC BOILER AND HEAT EXCHANGER
Electric boilers and heat exchange devices are all devices that provide hot water power, and the constraints are presented by equations (34) to (37).

4) CONSTRAINTS OF COOLING AND HEATING DEVICES
The cooling and heating system consist of absorption chiller, heat pump and electric cooler. In terms of the absorption chiller, it will use the energy from smoke and hot water to drive the cooling process. For the heat pump, it can consume electricity to implement either cooling or heating process. In principle, the cooling process can be achieved by the absorption chiller, heat pump and electric cooler whilst the heating process can only be implemented by the heat pump. Two variables µ c and µ h are employed to denote the working status of the cooling and heating process; µ c becomes 1 when cooling process starts and µ h turns 1 when heating process begins. The power relationship among the absorption chiller, heat pump and electric cooler and the constraints are defined in equations (38) to (48).
x Absorption chiller y Heat pump z Electric cooler { Cooling and heating process variables Note: if the heating load L ta = 0, Note: if the cooling load L ca = 0, By default, the IES system either purchases or sell electricity from/to the power grid at any given period and there will be no simultaneous bi-directional electrical power flow. Two variables µ b and µ s can be used to represent the electricity buying and selling process: µ b = 1 means the IES system is buying electricity from the grid and µ s = 1 suggests selling electricity to the grid. The constraints of electrical power flow are shown in equations (49) to (52).

III. SYSTEM DAY-AHEAD ECONOMIC DISPATCH
The day-ahead economic dispatch of the system described in this paper is to predict the output and load data of renewable energy in each period of the next day in the previous day. The optimization algorithm aims to minimize the 24-hour system operation cost. In general, a single day is divided into 24 time slots and all devices of the IES system will be dispatched based on the forecast of system electrical load, hot water load, (air) heating and cooling load, wind generation and solar generation. In this section, six groups of forecast data will be used as a benchmark for system planning. It is worth noting that there will be error between the actual and forecast data and CVaR and RO approaches are proposed to reduce the impact of the error and improve the accurateness of the algorithm.

A. DETERMINISTIC METHOD: MIXED INTEGER LINEAR PROGRAMING
Without considering the uncertainty of the predicted value, that is, without considering the risk caused by the difference between the actual value and the predicted value. The system operation cost can be divided into three categories: cost of buying natural gas (C ng ), electricity consumption cost (C grid.b ) and revenue obtained from selling surplus electricity to grid (C grid.s ), where the unit price of purchasing natural gas and selling surplus electricity are fixed. On the other hand, there are three price bands when buying electricity from grid, depending on the time of buying. Equations (53) to (55) show how the cost value of each category can be derived.
Mixed integer linear programing needs to meet the equipment operation constraints and strictly meet the bus energy balance constraints. The objective function to minimize the system operation cost and subjected to various constraints can be expressed as In the system described in this paper, all predicted values are uncertain. The output of PV and wind power is affected by the weather and the load is affected by people. These uncertainties lead to the uncertainty of the operation of all equipment in the system. In the current planning, there will be a risk that the output cannot meet the load demand and the system operation cost is too high. To this end, the CVaR model and the robust optimization model are separately built to plan the system and reduce the risks.

1) UNCERTAINTY METHOD 1: CVaR MODEL
Since there is no historical data, a simple Monte Carlo simulation method is used here to calculate CVaR. It is assumed the number of Monte-Carlo simulations is N and the confidence level is α. All known values will follow the normal distribution with expected uncertainty of E P and standard deviation of σ t . The forecast values at a later time slot would be less accurate that those at an earlier slot. Assume Y is a real number greater than 24, σ t can be calculated as: When generating a random number based on the probability distribution of the uncertainty, it is necessary to satisfy that all the uncertainties ω are positive values, and that the power of the photovoltaic and wind power does not exceed the maximum value.
The calculation method of the system running cost C i Bill for each simulation is the same as the objective function of Deterministic method, and the expectation of the system running cost is the average value of each simulation running cost, which are expressed in equations (60) and (61), respectively.
In this paper, the objective function of CVaR is defined as the CVaR whose operating cost exceeds the expected operating cost. In conventional CVaR method, a prerequisite is the value of VaR which make it much more complex to solve the equation. According to [25], the CVaR could be derived without knowing the value of VaR. This approach is applied in this study to obtain the CVaR when the actual system operation cost exceeds the expected value. Two intermediate variables ξ and φ i are used to calculate CVaR and ξ represents VaR. The proposed objective function is equivalent to CVaR, shown in equations (62) to (64).
After combining CVaR, the objective function of the system planning is changed to the sum of expected operation cost and CVaR (when the operation cost outnumbers the expected value). At this point, it becomes a multi-objective plan. We need to define the weight of the CVaR. The weight coefficient β of the CVaR is in the range of 0 and 1.
The CVaR model also needs to meet the energy bus constraints and equipment operation constraints. Unlike the common optimization, the known quantities of each simulation in the CVaR model are different. Each energy bus cannot take the equal sign, and the inequality is used to satisfy the supply. At the same time, the constraints (60) and (61) for solving CVaR shall also be satisfied. The objective function and constraints of the CVaR model are: In this paper, a single-stage robust optimization model is used, without knowing the probability distribution of uncertain quantities. Instead, the robust optimization model is built by describing the fluctuation of parameters through an uncertain set, such as load, wind power and photovoltaic power generation. The uncertainty fluctuation is denoted as E P . The amplitude of the fluctuation is τ t , and the uncertainty of defining the uncertainty is r t , which is used to control the amplitude of the uncertainty fluctuation. Similar to Uncertainty method 1, the uncertainty in the robust optimization model must also be a positive number. The value of photovoltaic and wind power should be less than the maximum value. Then the uncertainty set of the uncertainty ω can be expressed as Because robust optimization with box set tends to be conservative, a robust control factor t is introduced. Equation (65) states that the control factor t is equivalent to the sum of uncertainty level r t for all uncertain variables and the robustness can be altered by manipulating the magnitude of t .
The decision variable is represented by x. f and g are functions of x where f is identical to the objective function of Deterministic method and g is a set of common planning constraints and uncertainty constraints as described in (66) ∼ (68). The busbar constraint is also changed so that the supply will be either greater than or equal to the demand. ω is the uncertainty whilst is the uncertainty set. The objective function and constraints of the robust optimization can be presented as: Constraints: g (x, ω) ≤ 0∀ω ∈   Table 1 lists the parameters of each device within a HCNG penetrated IES. The unit price of buying natural gas and selling electricity is 3.1/Nm 3 and 0.1/Nm 3 , respectively. Figure 2 illustrates the three price bands and the corresponding time periods, when electricity is purchased from grid. The forecast load, wind generation and solar generation are plotted in Figure 3 and Figure 4. For the hydrogen storage tank, the initial amount of hydrogen being stored is 10 Nm 3 . The fuel cell used in the case study is the conventional model which means the heat produced from the electricity generation process will not be reused. Analysis of all cases is encoded using YALMIP and the equations are solved using CPLEX 12.8.0. In particular, Uncertainty Method 2 uses the Robust optimization module in YALMIP to solve.

A. DETERMINISTIC METHOD
Using the methods and constraints defined in Chapter II and III, the simulation results indicates the system operation cost of a single day is 1129.5. The program running time is 0.5579s. Figure 3 to Figure 5 show the optimization results of three important busbars in the system respectively -the electric busbar, the hot water busbar and the air busbar.

1) OPTIMIZATION RESULT ANALYSIS
Electric busbar: The electrical energy within the system is mainly derived from the CHP, the photovoltaics and the wind power. From 6 to 22, the CHP keeps operating and provides the base electrical energy. When the photovoltaics (from 6 to 18) and the wind power (19 -22) are available, they can contribute to the total electricity generation and thus reduce the output from the CHP. In the period when CHP is not running, the system's electrical power is mainly generated from the wind power.
Hot water busbar: When the CHP is in operation from 6 to 22, the hot water load is high (larger than 30 kW) and most of the hot water power comes from the CHP. During this period, the hot water power provided by the heat exchanger and the electric boiler is very small. When the CHP is not running, the hot water demand is met by the electric boiler.
Air busbar: The cooling power is mainly provided by the energy recovered from the hot water and the flue gas through the absorption chiller. Note that the vast majority of the power contained in the hot water and all power within flue gas are supplied by the CHP. Only a small amount of cooling power is provided by the electric cooler.

2) UTILIZATION OF PHOTOVOLTAIC WIND POWER
As mentioned earlier, the power of photovoltaics is fully utilized during the day. During the night time, the power of wind is high and is greater than the load of the system. The electric hydrogen production equipment can make full use of such additional wind power by triggering the hydrogen production process to fill the hydrogen storage tank. In this way, the required amount of hydrogen produced from other energy sources during the day time can be reduced and thus the required energy within the system can be reduced. Figure 6 plots the hydrogen volume change in a hydrogen storage tank. In the period of 1 to 7, the system uses wind power to produce hydrogen which can provide energy during the day. From the simulation results, the wind power is sold to the grid for a total amount of 58.2 kWh from 22 to 24, and the utilization rate of wind power was 89.4%.
When the electrolysis cell, hydrogen storage tank and hydrogen fuel cell are removed for the simulation, the operating cost is 1168.9, which is higher than that of the system containing these devices ( 1129.5). Moreover, a larger amount of wind power is sold to the power grid and the utilization rate of wind power is only 28.9%. This results in the fact that more power needs to be purchased from the power grid during the peak hours during the day time. Figure 7 shows the influence of whether the system contains hydrogen production and storage equipment on wind energy utilization, power grid purchase and natural gas consumption. Obviously, after the employment of hydrogen and hydrogen storage equipment, the system operating costs are reduced and the utilization of wind energy within the IES is more efficient. As a result, the electricity demand from the grid is reduced.

3) HYDROGEN UTILIZATION
It is mentioned in [6] that the change of hydrogen content in the fuel of CHP unit will change the electric-heat ratio of CHP output power, and the increase of hydrogen content will raise the electric-heat ratio. The hydrogen content change of CHP fuel in the system is shown in Figure 8. In order to study the relationship between load electro-heating ratio and CHP output electro-heating ratio, the electro-heating ratio of the system load is defined as electric load divided by the sum of hot water load and air cold load. From previous analysis, the loading on electric busbar, hot water busbar and air busbar can reflect the operation of the system. It can be observed from Figure 8 that the CHP output electric-heating ratio and the system load electric-heating ratio follows the same pattern and the load electric-heating ratio is larger than the CHP output electric-heating ratio because the electric energy in the load also includes the photovoltaic output. During the period from 7 to 10, although the load electro-thermal ratio is changing, the value stays at high level and the hydrogen content of the CHP fuel reaches the maximum value which is 20%. In the period of 11 to 17, the photovoltaic output grows and thus the electric-heating ratio of CHP decreases. Between 20 pm and 21 pm, the energy ratio provided by CHP is not particularly large, and the load-to-heat ratio is difficult to reflect the CHP output electric-to-heat ratio. In the case where electric energy is mainly supplied by the CHP, when the hydrogen content in the CHP fuel does not reach the peak value, the hydrogen content of the CHP fuel changes as the load electro-thermal ratio fluctuates. In this way, the system can plan the hydrogen content of the CHP fuel according to the load electro-thermal ratio so that the hydrogen can be utilized more efficiently.

4) INFLUENCE OF SOFC ON IES
Among all types of fuel cell, the SOFC operates at the highest temperature. SOFC is known as high temperature fuel cell which can provide large amount of heat that can be reused for CHP system and thus improve the efficiency of fuel utilization. However, if the simulation case uses SOFC and optimization plan merely based on forecast data, the system operation cost is 1128.2 which is only 1.3 less than the case where conventional hydrogen fuel cell is used. It should be noted that the cost of SOFC is very high. In addition, the SOFC consumes more hydrogen than the CHP and only 23.7% of hydrogen will go to the CHP, which will dramatically reduce the proportion of hydrogen within the CHP, as illustrated in Figure 9.

1) UNCERTAINTY METHOD 1
In the Uncertainty method 1 described in section III, Y is configured as 100 and the confidence level α as 0.95. The number of Monte Carlo simulations is 500 and they are VOLUME 7, 2019 FIGURE 9. Volume proportion of hydrogen within CHP when using SOFC and conventional hydrogen fuel cell. run under different CVaR weights. Figure 10 illustrates the minimum operating cost expectations and the CVaR. When CVaR is not considered (i.e. CVaR weight = 0), the operating cost is expected to be 1214.9, but the CVaR is very large at 422.4. When the weight of CVaR gradually increases from 0, although the expectation of the operating expense only slightly increase by less than 20, the CVaR decreases to approximately 100. After the weight exceeds 0, the values of the cost expectation and the CVaR do not obviously change. From the optimization results, the power provided by the electrical busbar, hot water busbar and air busbar is greater than the load, and the later the time period is, the greater the gap between the provided power and the load will be.
In order to balance the electric busbar, the purchased power P grid.b and the sold power P grid.s are set to be change freely in this section and in the robust optimization section (Because the system does not have access to the heating network and the air network, the balance of the hot water busbar and the balance of the air busbar are not considered at present. The only requirement is to satisfy the demand that is more than or equal to the electricity supply). Because the operating cost of the system is only related to the purchase and sale of electric power as well as buying of natural gas (other related costs are negligible), the loading of the electrical bus can reflect the change of the operating cost of the entire system. Take the case where the CVaR weight is 1 as an example. Figure 11 presents the optimization result of the electrical busbar when the interaction with the grid is excluded. In the period when the output of the photovoltaic and wind power is small and the CHP operation is required, the electric power provided by the system is less than the load. This insufficient power should have been supplemented by purchasing electricity from the grid. When the wind power is relatively large, the excess electricity should have been sold to the grid. When the electrical load fluctuates within appropriate range, the scheduled operating plan can meet the operational requirements.

2) UNCERTAINTY METHOD 2
Uncertainty method 2 described in section III is used for day-ahead planning. In order to facilitate the comparison with Uncertainty method 1, the uncertainty set of the uncertainty is set as the confidence interval of the uncertainty in the upper section of Uncertainty method 1. When the confidence level α is 0.95, the confidence range is [E P − 1.96σ t , E P + 1.96σ t ] with resulted uncertainty fluctuation τ t of 1.96σ t . The system is simulated by taking a number of different t values between 0 to 6. The expected values of the running cost varying with t are shown in Figure 12. When t is 0, Uncertainty method 2 is converted to Deterministic method, and the running cost is also 1129.5. The operating cost expectation increases when t grows, and the slope in the interval (1, 2) becomes smaller. When t is 2, the operating cost expectation reaches the maximum value of 1687.2, after which it does not increase when t keep growing beyond 2. Similar to Uncertainty method 1, the optimization result of Uncertainty method 2 is that the power provided by the electric busbar, the hot water busbar  and the air busbar is greater than the load, and the gap between the provided power and the load will become larger in later period. Figure 13 shows the electrical bus optimization results using Uncertainty method 2 when the grid interaction is removed ( t = 0.25). Similar to Uncertainty method 1, when the CHP is not working, the wind power is larger than the load, and the system shall sell electricity to the grid. However, the CHP operating period is different from Uncertainty method 1. In the period of 10 to 17, the power provided by the IES is greater than the load. If the load fluctuates downward or the wind power and PV output fluctuates upwards, the excess power of the system needs to be sold to the grid, which is not efficient from the perspective of energy and economy. Therefore, this optimization scheme tends to be conservative during the period of CHP operation.

3) ANALYSIS OF WIND POWER UTILIZATION AND HYDROGEN UTILIZATION
In this section, Uncertainty method 1 with a CVaR weight of 1 and Uncertainty method 2 with a robustness control factor of 0.25 are used for wind power utilization and hydrogen utilization analysis.
x Analysis of optimization results of each busbar Figure 12 and Figure 13 show the results of the electrical busbar optimization using Uncertainty method 1 and Uncertainty method 2 after removing the interaction with the grid. If Deterministic method is used, the electrical energy of the system is mainly derived from the CHP, photovoltaics and wind power. In all three methods, the period of CHP operation is 6 -22. The power of the system is mostly derived from CHP during the period of CHP operation. Whilst in the period when CHP is not running, the power of the system is mainly produced from wind power. The hot water busbar and the air busbar optimization results of Uncertainty method 1 and Uncertainty method 2 are similar to Deterministic method.
y Utilization of wind power All three methods generate hydrogen from the electrolysis cell at night, which is stored in the hydrogen storage tank. The total amount of hydrogen produced is shown in Table 2 and it can be observed that all three methods generate similar amount of hydrogen. Referring to the operation plan, the period of hydrogen generation is from 23 pm to 8 am the next day.  z Hydrogen utilization In Deterministic method, the hydrogen content in the CHP fuel varies with the load electro-thermal ratio, and this relationship is also identified in Uncertainty method 1 and Uncertainty method 2. As shown in Figure 14, the trend of VOLUME 7, 2019 Uncertainty method 2 and Deterministic method is almost the same in the period of 8 to 17, because the uncertain variables in Uncertainty method 2 are centered on the predicted values and there is a close relationship between them. The trend of Uncertainty method 1 is similar to Deterministic method only in the period of 12 to 15. This is because the planning of Uncertainty method 1 needs to meet all conditions of the 500 sets of Monte Carlo simulations and eventually it only reflects the statistical characteristics of the 500 simulations to some extent. Because various energy busbars must strictly abide by the characteristic where the supply is greater than demand, the changing pattern of the CHP fuel hydrogen content in Uncertainty method 1 is different from Deterministic method. If the 500 Monte Carlo simulations are independently optimized, the statistical trend should be the same as Deterministic method for most of the time.

C. COMPARISON OF OPTIMIZATION RESULTS
In order to validate the resistance of proposed optimization approaches (i.e. mixed integer linear programing, CVaR model, robust optimization model) to the fluctuation in the forecast data, two groups of actual data are established. Out of these data, one group is built by introducing massive variation to the forecast data whilst the other group is constructed by adding tiny fluctuation. The values are shown in Figure 15. From Figure 10 and Figure 11, the expected system operation cost is similar when CVaR weight is 1 (in approach: Uncertainty method 1) and when t is 0.25 (in approach: Uncertainty method 2). Hence, these two optimization scenarios are selected for the investigation of fluctuation resistance. There are three main points of comparison: -whether the optimization approach can meet the demand of each energy busbar given that the supply is greater than or equal to the demand -the running cost -the speed of the program running.

1) VERIFICATION THAT THE SUPPLY ON THE BUSBAR IS GREATER THAN OR EQUAL TO THE DEMAND
The optimization results/plans need to be validated whether they can satisfy the operation requirements. More specifically, the balancing equations of the energy-carrier busbars (1) to (4) are transformed to equations (69) to (71). If these equations can be satisfied for all time intervals, the optimization results/plans can satisfy the operation requirements. Note that equation (66) of electrical busbar is slightly different because it is assumed the buying and selling of electricity can flexibly change without restriction.
−P max grid.b ≤ P PV + P W + P HFC + P CHP − (L E + P HP + P EC + P EL + P EB ) ≤ P max grid.s (69) Q HFC + Q HE.water + Q CHP.water + Q EB − (L HW + Q AC.water ) ≥ 0 (70) The operating costs of the system are mainly composed of electricity purchasing cost, natural gas purchasing cost and electricity selling revenue. The amount of electricity being purchased and sold can be derived from equation (69). Furthermore, the cost of buying natural gas and electricity as well as the revenue of selling electricity can be calculated using equations (53) to (55). The resulted system operation cost is illustrated by equation (73).
The comparison results of section (1) and (2) are shown in Table 3. The * indicates that the supply cannot meet the demand. Regardless of the fluctuations, the running cost of Deterministic method is always the lowest. However, the associated planning scheme cannot meet the operational needs. By simply analyzing the uncertainty in load and output, if the load increases or the output decreases, both by a specific amount, Deterministic method scheme will not be able to meet the demand.
When encountering the small fluctuations described in this paper, Uncertainty method 1 with weight 1 and Uncertainty method 2 with t of 0.25 can meet the operational requirements. The cost of Uncertainty method 2 is 0.62% less than that of Uncertainty method 1. When large fluctuations exist, the operating cost of Uncertainty method 2 is 1.52% less than that of Uncertainty method 2, but Uncertainty method 2 with t of 0.25 cannot meet the operational requirements of the system. Although Uncertainty method 1 runs at a higher cost, it can operate the system regardless of the fluctuation size.
In order to find a robust optimization scheme that can resist the aforementioned large fluctuations, the data in Figure 14 are verified from low t to high t , and it is discovered that when t = 1, the planning scheme can satisfy the demand in the case of large fluctuations. However, the operating cost at this time is relatively high, and the operating costs in the cases of no fluctuation, small fluctuations and large fluctuations are 19.75%, 16.54% and 9.04% higher than the operating expenses of Uncertainty method 1, respectively. In order to analyze why Uncertainty method 2 with t of 0.25 cannot meet the system operation requirements during large fluctuations, the air busbar is selected as the investigation object. The reason for choosing the air busbar is that this bus is the simplest busbars that is with uncertainties. Only the load is uncertain and this bus is enough to reflect the characteristics of the system planning. Figure 16 plots the cooling power and the cooling load corresponding to a mixed integer linear programing, a CVaR model with weight 1, a robust optimization model with t of 0.25 and a robust optimization model with t of 1 when large fluctuations are encountered. In three intervals near 12, 16, and 18, the cooling power of Deterministic method and Uncertainty method 2 with t of 0.25 is less than the load. Consequently, Deterministic method and Uncertainty method 2 with t of 0.25 cannot satisfy the system operating requirements. For the other two optimization methods satisfying the system operation requirements, Uncertainty method 1 with weight 1 has a higher output than Uncertainty method 2 with t of 1 at any time and can thus resist greater fluctuations.

3) PROGRAM RUNNING TIME COMPARISON
During the planning process, the speed of the program should also be considered, in addition to the performance of the method. Table 4 summarizes the average running time of the program in which all three methods run 10 times under different conditions. Deterministic method run time is less than 1 second. The average running time of Uncertainty method 2 is only 0.2% of Uncertainty method 1 when Monte Carlo simulation is required to run 500 times. Even if only 10 times of Monte Carlo simulations are performed, the time used is still more than twice that of Uncertainty method 2. Moreover, too few Monte Carlo simulations will result in uneven distribution of random numbers, which will affect the reliability of the simulation results. The recommended number of Monte Carlo simulations is between 200 to 500 times.

4) SUMMARY OF COMPARISON RESULTS
According to the comparison results, although Deterministic method is the fastest, its performance is the worst and it is not suitable for economic dispatch. When the prediction accuracy is high and the external influence on the system is small, Uncertainty method 2 with appropriate robustness control factor is more appropriate because it has the advantages of relatively fast computing speed and relatively low operating cost of the system. When the environment is volatile and the accuracy of prediction is difficult to guarantee, the speed of the program should be de-prioritized, and Uncertainty method 1 with better risk resistance should be selected.

V. CONCLUSION
In order to absorb wind energy at night and reduce the amount of electricity purchased from the grid during the daytime, this paper deploys hydrogen generation and storage equipment to the integrated energy system, and supplies hydrogen to the CHP unit and the hydrogen fuel cells. The hydrogen content of the fuel can affect the output electro-thermal ratio of the CHP which can be used to design control strategy to achieve the most efficient use of hydrogen. Taking into account the uncertainty of wind power, photovoltaics and load, the CVaR and the robust optimization can be used to reduce the risk. By simulating a set of predicted values for wind power, photovoltaics, and loads, and then constructing two sets of data based on the predicted values, the following conclusions can be made.
1) Compared with the integrated energy system without hydrogen equipment, the HCNG penetrated integrated energy system described in this paper significantly improves the night time wind power utilization rate by using wind power to generate hydrogen at night, and consumes the corresponding stored hydrogen during the day time to reduce the electricity demand from the grid during the day.
2) The content of hydrogen in the CHP fuel will change with the variation of the load electro-thermal ratio.
If the CHP output ratio is extremely high and the hydrogen content does not reach the maximum value, the increase of the load electro-thermal ratio will increase the hydrogen content. Using this relationship, the efficiency of using hydrogen can be maximized and the operating cost can be reduced. 3) When the actual value significantly deviates from the predicted value, Uncertainty method 1 has lower system operating cost and higher risk resistance capability, while Uncertainty method 2 has lower operating cost when the fluctuation is small. Considering the performance and the computing speed, it is ideal to use Uncertainty method 2 when the prediction is accurate. By contrast, Uncertainty method 1 performs better when it is difficult to guarantee the prediction accuracy. CAMPBELL BOOTH received the B.Eng. and Ph.D. degrees in electrical and electronic engineering from the University of Strathclyde, Glasgow, U.K., in 1991 and 1996, respectively, where he is currently a Professor and the Head of the Department for Electronic and Electrical Engineering. His research interests include power system protection, plant condition monitoring and intelligent asset management, applications of intelligent system techniques to power system monitoring, protection, and control, knowledge management, and decision support systems. VOLUME 7, 2019