# Robust Stochastic Dynamic Optimal Power Flow Model of Electricity-Gas Integrated Energy System considering Wind Power Uncertainty

Authors: Zhengfeng Qin; Xiaoqing Bai; Xiangyang Su
Journal: Complexity
 (2020)

---

## Abstract

The application of gas turbines and power to gas equipment deepens the coupling relationship between power systems and natural gas systems and provides a new way to absorb the uncertain wind power as well. The traditional stochastic optimization and robust optimization algorithms have some limitations and deficiencies in dealing with the uncertainty of wind power output. Therefore, we propose a robust stochastic optimization (RSO) model to solve the dynamic optimal power flow model for electricity-gas integrated energy systems (IES) considering wind power uncertainty, where the ambiguity set of wind power output is constructed based on Wasserstein distance. Then, the Wasserstein ambiguity set is affined to the eventwise ambiguity set, and the proposed RSO model is transformed into a mixed-integer programming model, which can be solved rapidly and accurately using commercial solvers. Numerical results for EG-4 and EG-118 systems verify the rationality and effectiveness of the proposed model.


---

## Body

## 1. Introduction

With the development and practical application of the theory of the “integrated energy system (IES),” it has become a research upsurge to realize the mutual transformation of multienergy flows over the worldwide [1–8]. Compared with other primary energies, natural gas and wind power have the characteristics of clean and environmental protection. Thus, it is of practical significance to study the interaction among natural gas, wind power, and the power system.Some research studies about the coordinated planning and optimal operation of the electricity-gas IES with gas turbine and power to gas equipment as the coupling relationship [9, 10] are in literature. Fan et al. and Dai et al. [11, 12] studied the optimal power flow (OPF) problem of the electricity-gas IES with natural gas network constraints. Still, the impact of the uncertainty of wind power output on the optimal operation of the integrated energy system is not considered. It is necessary to take effective uncertainty optimization techniques and solve approaches with risk consideration for the OPF problems of the electricity-gas IES.How to deal with the uncertainty of wind power output is one of the important difficulties in the collaborative optimization of electricity-gas IES. Stochastic optimization (SO) [13] and robust optimization (RO) [14, 15] are two effective methods. The SO method is described by probability distribution with a heavy computational burden. At the same time, it is difficult to obtain accurate probability distribution in practice, which may lead to inaccurate models [16]. RO algorithm is characterized by boundary parameters of uncertain factors. The computational cost of RO is low, but it is hard to select a suitable robust set, and the decision results are conservative in most cases [17]. In recent years, scholars have tried to use a distributionally robust optimization method (DRO) [18, 19] to deal with uncertain problems. This method combines stochastic programming with robust optimization to find the worst probability distribution function under the condition of satisfying uncertain parameter information. However, the solution process for this method is very complex.In the DRO method, the construction of the ambiguity set is very important. The form of the ambiguity set will not only affect the transformation and solution of the distributed robust optimization model but also affect the conservative degree of the optimization results. The construction of ambiguity sets for DRO can be divided into the following categories: constructing an ambiguity set based on moment information of uncertainty [20]; constructing an ambiguity set based on probability distribution functions of uncertainties; and constructing an ambiguity set in other ways [21]. Although the ambiguity set based on the first-order and second-order moment information has the advantages of equivalent model processing and simple set construction, the real probability distribution information of uncertainty is limited, and the results are conservative. To make full use of sample data and reduce the conservatism of model decision-making, the ambiguity set based on probability distribution function can provide more uncertain probability distribution information, thus making the scheduling decision-making process more accurate. However, when the support sets of the two probability distribution functions do not overlap or overlap very little, the ambiguity set based on Kullback–Leibler (KL) divergence may not be meaningful [22], while the ambiguity set based on Wasserstein distance (WD) is feasible. The ambiguity set constructed by WD can still reflect the distance between two probability distributions under the above circumstances and can measure the “distance” between any two probability distributions more accurately [23, 24]. Owing to the current research on constructing ambiguity sets based on WD focuses on the Wasserstein distance of 1-norm, the other forms of Wasserstein ambiguity sets are not involved in. Therefore, this paper uses RSO algorithm to calculate the influence of the 1-norm ambiguity set and ∞−norm Wasserstein distance on the test system, in order to find a better ambiguity set.For different ambiguity sets, DRO has different decision models and corresponding processing methods, which may lead to the model and solution only applicable to specific problems in the electricity-gas IES and lacking generality. In order to find a general model expression, the RSO model is proposed [25]. The author proposes to construct DRO with an eventwise ambiguity set and proves that the eventwise ambiguity set has good generality. Besides, the eventwise ambiguity set can describe the discrete distribution, moment-based distributed ambiguity set, data-driven K-means ambiguity set, and Wasserstein ambiguity set. The RSO model is applied to the economic dispatch problem of the power system considering wind power uncertainty for the first time [26], and the validity and accuracy of the RSO model are verified by simulation.This paper proposed an RSO model to solve the dynamic optimal power flow model for electricity-gas IES considering wind power uncertainty, where the ambiguity set of wind power output is constructed based on Wasserstein distance and then affined to the eventwise ambiguity set to solve. The main contributions are as follows:(1)
The ambiguity set of wind power output is constructed based on Wasserstein distance and then affined to the eventwise ambiguity set of the RSO model, which can be solved by commercial solvers such as CPLEX, GUROBI, and MOSEK.(2)
Aiming at the uncertainty of wind power output, Wasserstein ambiguity sets based on 1-norm and∞−norm are constructed, respectively, to study the influence of wind power penetration, adjustment cost coefficient of wind units, and sample size on the operation of the electricity-gas IES. The simulation results verify that the Wasserstein ambiguity set based on ∞−norm has a more economic optimal solution.(3)
The AC power flow model is approximated as a decoupled linear power flow model to improve the feasibility of the generation plan obtained from the model.The rest of this paper is organized as follows. Section2 presents a robust stochastic OPF model for the electrical-gas IES connecting wind power uncertainty where the natural gas flow is modeled as an enhanced second-order cone (SOC) relaxation formulations. In Section 3, considering the uncertainty of wind power output, the model is established based on RSO. In Section 4, the case studies are conducted to demonstrate the effectiveness of the proposed model. The conclusion is presented in Section 5.


## 2. Dynamic Optimal Power Flow Model of Electricity-Gas IES

Electricity-gas IES is a new integrated multielement network system, which integrates the data of the power network, natural gas network, and new energy power generation. The power system and natural gas system are coupled by energy conversion equipment, such as gas turbines and the power to gas (PtG) equipment to realize bilateral energy conversion and heterogeneous energy complementary.With the connecting of the natural gas system and wind generators to the electricity system, the coordinated and optimized operation of the electricity-gas IES is much more complex than considering the power system and gas system separately. The objective function is to minimize the total generation costs, which include the operation cost of thermal power units, the natural gas consumption cost of the natural gas subsystem, and the worst-case expected adjusting cost. It can be formulated as(1)min∑t=1T∑i=1NpaiPGi,t2+biPGi,t+ci+∑j=1NgasCgjwgj,t+maxp∈Ω∑s=1Spsmin∑k=1Nwck,twindPk,t,sw0−Pk,t,sw,where T is the total number of time periods in the dispatching cycle; Np, Ngas, and Nw are the sets of the thermal power units, the gas supply point, and the all distributed wind, respectively; PGi,t is the active power of the thermal power generator of bus i in period t; Cgj is the price of gas supply point i; wgj,t is the injection flow at the ith gas source point in period t; ai, bi, and ci are the consumption characteristic parameters of the thermal power unit at bus i; ck,twind is the adjusting cost coefficient of wind unit k in period t; ps is the probability of each discrete scenario s; and Pk,tw0 and Pk,tw are the predicted output and actual output of wind power in period t of the wind turbine k in the scenario s.


### 2.1. Constraints of the Electricity System

#### 2.1.1. Power Balance Equation Constraint

(2)Pi,t=PGi,t+Pgj,t+Pk,tw−Pi,tPtG−Pdi,t,=Vi,t∑j=1nVj,tGijcosθij,t+Bijsinθij,t,(3)Qi,t=QGi,t−Qdi,t,=Vi,t∑j=1nVj,tGijsinθij,t−Bijcosθij,t.Note that formulae (2) and (3) are a nonconvex quadratic optimization problem. In order to improve the efficiency of the model, the nonconvex nonlinear AC power flow equation is approximated to a decoupled linear power flow equation [27]. First of all, expand equations (2) and (3) to obtain the following results:(4)Pi,t=Vi,t∑j=1nVj,tGijcosθij,t+Bijsinθij,t=giiVi,t2+∑j=1,j≠ingijVi,tVi,t−Vj,tcosθij,t−bijVi,tVj,tsinθij,t,(5)Qi,t=Vi,t∑j=1nVj,tGijsinθij,t−Bijcosθij,t=−biiVi,t2+∑j=1,j≠inbijVi,tVi,t−Vj,tcosθij,t+gijVi,tVj,tsinθij,t.The nonlinear term in the above formula is approximated as follows:(6)Vi,t2≈Vi,t,Vi,tVj,tsinθij,t≈θi,t−θj,t,cosθij,t≈1,Vi,tVi,t−Vj,tcosθij,t≈Vi,t−Vj,t.Then, (4) and (5) can be simplified as(7)Pi,t=∑j=1nGijVj,t−∑j=1nBij′θj,t,Qi,t=−∑j=1nBijVj,t−∑j=1nGijθj,t.Similarly, the line power can be expressed as(8)Pij,tl=GijVi,t−Vj,t−Bijθi,t−θj,t.


#### 2.1.2. Voltage Constraint

(9)Vimin≤Vi,t≤Vimax.


#### 2.1.3. Bound of Generation Power

(10)PGimin≤PGi,t≤PGimax,QGimin≤QGi,t≤QGimax,Pijl,min≤Pij,tl≤Pijl,max.


#### 2.1.4. Climbing Constraints

(11)PGi,t−PGi,t−1≤RUi,PGi,t−1−PGi,t≤RDi,where Gij+jBij is the element of power system admittance matrix between branches line i,j; Gij′+jBij′ represents the element of admittance matrix without self-admittance; gij+jbij is the admittance of line i,j; Vi,t is the voltage at bus i in period t ; Pi,t and Qi,t are the active power and reactive power injected by node i in period t, respectively; Pgi,t is the active power output of the gas turbine i in period t; Pi,tPtG is the active power output of the power to gas equipment i in period t; PGi,t and QGi,t are the active power and reactive power of the thermal power generator i in period t, respectively; Pdi,t and Qdi,t are the active power load demand and reactive load demand at bus i in period t, respectively; RUi and RDi are the climbing rate and landslide rate of thermal power units i; Pij,tl is the line power of line i,j in period t; PGimin and PGimax are the lower and upper limits of PGi; QGimin and QGimax are the lower and upper limits of QGi,t; ViminandVimax are the lower and upper limits of the voltage Vi,t at bus i; and Pijl,min and Pijl,max are the lower and upper limits of Pij,tl.


### 2.2. Constraints of the Natural Gas System

Generally, a natural gas system is composed of gas source points, supply pipes, compressors, load, and convert primary energy into electric energy through the gas turbines. If the natural gas system operates steadily, the steady-state natural gas flow through the pipeline can be assumed as a function of the upstream and downstream pressure. When the natural gas flow from nodei to node j, the gas system operational constraints include natural gas flow constraints of natural gas pipeline in (12)–(25), nodal gas balance in (26), and operating limits in (27)–(29).


#### 2.2.1. Natural Gas Flow Constraints of Natural Gas Pipeline

(12)fi,j,t−=12fi,j,t−fj,i,t,(13)fi,j,t2−Ki,j2=pi,t2−pj,t2,fi,j,t>0.Otherwise,(14)fi,j,t2−Ki,j2=pj,t2−pi,t2,fi,j,t<0,where fi,j,t is the natural gas flow in the pipeline i,j in period t; f¯i,j,t is the average natural gas flow in the pipeline i,j in period t; Ki,j is a pipe constant, which is related to pipe length, diameter, friction coefficient, and gas compression coefficient; pi,t is the gas pressure at the node i in period t. Noted that equations (13) and (14) can exist precisely in the form of a SOC relaxation; thus, they can be expressed as(15)fi,j,t2−Ki,j2≤pi,t2−pj,t2,fi,j,t>0,(16)fi,j,t2−Ki,j2≥pi,t2−pj,t2,fi,j,t<0.A standard technique to convexity such a model is to relax (15), which can be solved rapidly and accurately by commercial solvers. However, the solutions that are obtained from such relaxation may yield nontrivial violations of the constraint set (13).To make the SOC model tractable while obtaining high-quality solutions, Chen et al. [28] proposed an enhanced SOC relaxation for the convex relaxation (15). Defining two sets of variables, ai,j,t and bi,j,t, which are defined as the sum and difference of the pressures at both ends of each pipe,(17)ai,j,t=pi,t+pj,t,(18)bi,j,t=pi,t−pj,t.Then, defining two new sets of auxiliary variables,γi,j,t and λi,j,t, and the convex relaxation of (16) is given by(19)γi,j,tKi,j2≥λi,j,t,(20)γi,j,t≥fi,j,t2−,(21)γi,j,t≤fi,jmax+fi,jminf¯i,j,t−fi,jmaxfi,jmin,(22)λi,j,t≥ai,jminbi,j,t+bi,jminai,j,t−ai,jminbi,jmin,(23)λi,j,t≥ai,jmaxbi,j,t+bi,jmaxai,j,t−ai,jmaxbi,jmin,(24)λi,j,t≤ai,jminbi,j,t+bi,jmaxai,j,t−ai,jminbi,jmax,(25)λi,j,t≤ai,jmaxbi,j,t+bi,jminai,j,t−ai,jmaxbi,jmin.


#### 2.2.2. Gas Nodal Gas Flow Balance

(26)wgi,t+wi,tP2G−wli,t−fgi,t=fij,t.


#### 2.2.3. Operating Limits

(27)fi,jmin≤fi,j,t≤fi,j,tmax,(28)pimin≤pi,t≤pimax,(29)wgimin≤wgi,t≤wgimax,where wgi,t is the injection flow at the ith gas source point in period t, wli,t is the load consumption flow at the node i in period t; wi,tP2G is the gas injection flow of the PtG equipment i in period t; fgi,t is the natural gas consumption flow of the gas turbine i connected to the power grid in period t; pimin and pimax are the lower and upper limits of pi,t; wiminandwimax are the lower and upper limits of wgi,t; fi,jmin and fi,jmax are the lower and upper limits of fi,j,t; ai,jminandai,jmin are the lower and upper limits of ai,j,t; and bi,jminandbi,jmin are the lower and upper limits of bi,j,t.


### 2.3. Constraints of Coupling Equipment

In this section, gas turbines and PtG equipment are studied as the coupling link of the gas-electric integrated energy system.


#### 2.3.1. Gas Turbines

Assuming that the gas turbine is connected between the power system nodei and the node j of the natural gas system, the conversion relationship between natural gas input and electric power output is as follows:(30)fgi,t=k1Pgi,t2+k2Pgi,t+k3,where k1, k2, and k3 are the consumption flow coefficients of the gas turbine.


#### 2.3.2. PtG Equipment

(31)wi,tP2G=ηP2GPi,tP2G,where ηP2G is the conversion efficiency of PtG equipment.Note that the proposed model mentioned above is a deterministic optimization (DO) model. However, the uncertainties of wind power output have not been considered at all. It is possible that the solution is not feasible in the actual operation. Therefore, the robust stochastic optimization (RSO) method will be used to solve the proposed optimization problem with uncertainties in Section3.


## 2.1. Constraints of the Electricity System

### 2.1.1. Power Balance Equation Constraint

(2)Pi,t=PGi,t+Pgj,t+Pk,tw−Pi,tPtG−Pdi,t,=Vi,t∑j=1nVj,tGijcosθij,t+Bijsinθij,t,(3)Qi,t=QGi,t−Qdi,t,=Vi,t∑j=1nVj,tGijsinθij,t−Bijcosθij,t.Note that formulae (2) and (3) are a nonconvex quadratic optimization problem. In order to improve the efficiency of the model, the nonconvex nonlinear AC power flow equation is approximated to a decoupled linear power flow equation [27]. First of all, expand equations (2) and (3) to obtain the following results:(4)Pi,t=Vi,t∑j=1nVj,tGijcosθij,t+Bijsinθij,t=giiVi,t2+∑j=1,j≠ingijVi,tVi,t−Vj,tcosθij,t−bijVi,tVj,tsinθij,t,(5)Qi,t=Vi,t∑j=1nVj,tGijsinθij,t−Bijcosθij,t=−biiVi,t2+∑j=1,j≠inbijVi,tVi,t−Vj,tcosθij,t+gijVi,tVj,tsinθij,t.The nonlinear term in the above formula is approximated as follows:(6)Vi,t2≈Vi,t,Vi,tVj,tsinθij,t≈θi,t−θj,t,cosθij,t≈1,Vi,tVi,t−Vj,tcosθij,t≈Vi,t−Vj,t.Then, (4) and (5) can be simplified as(7)Pi,t=∑j=1nGijVj,t−∑j=1nBij′θj,t,Qi,t=−∑j=1nBijVj,t−∑j=1nGijθj,t.Similarly, the line power can be expressed as(8)Pij,tl=GijVi,t−Vj,t−Bijθi,t−θj,t.


### 2.1.2. Voltage Constraint

(9)Vimin≤Vi,t≤Vimax.


### 2.1.3. Bound of Generation Power

(10)PGimin≤PGi,t≤PGimax,QGimin≤QGi,t≤QGimax,Pijl,min≤Pij,tl≤Pijl,max.


### 2.1.4. Climbing Constraints

(11)PGi,t−PGi,t−1≤RUi,PGi,t−1−PGi,t≤RDi,where Gij+jBij is the element of power system admittance matrix between branches line i,j; Gij′+jBij′ represents the element of admittance matrix without self-admittance; gij+jbij is the admittance of line i,j; Vi,t is the voltage at bus i in period t ; Pi,t and Qi,t are the active power and reactive power injected by node i in period t, respectively; Pgi,t is the active power output of the gas turbine i in period t; Pi,tPtG is the active power output of the power to gas equipment i in period t; PGi,t and QGi,t are the active power and reactive power of the thermal power generator i in period t, respectively; Pdi,t and Qdi,t are the active power load demand and reactive load demand at bus i in period t, respectively; RUi and RDi are the climbing rate and landslide rate of thermal power units i; Pij,tl is the line power of line i,j in period t; PGimin and PGimax are the lower and upper limits of PGi; QGimin and QGimax are the lower and upper limits of QGi,t; ViminandVimax are the lower and upper limits of the voltage Vi,t at bus i; and Pijl,min and Pijl,max are the lower and upper limits of Pij,tl.


## 2.1.1. Power Balance Equation Constraint

(2)Pi,t=PGi,t+Pgj,t+Pk,tw−Pi,tPtG−Pdi,t,=Vi,t∑j=1nVj,tGijcosθij,t+Bijsinθij,t,(3)Qi,t=QGi,t−Qdi,t,=Vi,t∑j=1nVj,tGijsinθij,t−Bijcosθij,t.Note that formulae (2) and (3) are a nonconvex quadratic optimization problem. In order to improve the efficiency of the model, the nonconvex nonlinear AC power flow equation is approximated to a decoupled linear power flow equation [27]. First of all, expand equations (2) and (3) to obtain the following results:(4)Pi,t=Vi,t∑j=1nVj,tGijcosθij,t+Bijsinθij,t=giiVi,t2+∑j=1,j≠ingijVi,tVi,t−Vj,tcosθij,t−bijVi,tVj,tsinθij,t,(5)Qi,t=Vi,t∑j=1nVj,tGijsinθij,t−Bijcosθij,t=−biiVi,t2+∑j=1,j≠inbijVi,tVi,t−Vj,tcosθij,t+gijVi,tVj,tsinθij,t.The nonlinear term in the above formula is approximated as follows:(6)Vi,t2≈Vi,t,Vi,tVj,tsinθij,t≈θi,t−θj,t,cosθij,t≈1,Vi,tVi,t−Vj,tcosθij,t≈Vi,t−Vj,t.Then, (4) and (5) can be simplified as(7)Pi,t=∑j=1nGijVj,t−∑j=1nBij′θj,t,Qi,t=−∑j=1nBijVj,t−∑j=1nGijθj,t.Similarly, the line power can be expressed as(8)Pij,tl=GijVi,t−Vj,t−Bijθi,t−θj,t.


## 2.1.2. Voltage Constraint

(9)Vimin≤Vi,t≤Vimax.


## 2.1.3. Bound of Generation Power

(10)PGimin≤PGi,t≤PGimax,QGimin≤QGi,t≤QGimax,Pijl,min≤Pij,tl≤Pijl,max.


## 2.1.4. Climbing Constraints

(11)PGi,t−PGi,t−1≤RUi,PGi,t−1−PGi,t≤RDi,where Gij+jBij is the element of power system admittance matrix between branches line i,j; Gij′+jBij′ represents the element of admittance matrix without self-admittance; gij+jbij is the admittance of line i,j; Vi,t is the voltage at bus i in period t ; Pi,t and Qi,t are the active power and reactive power injected by node i in period t, respectively; Pgi,t is the active power output of the gas turbine i in period t; Pi,tPtG is the active power output of the power to gas equipment i in period t; PGi,t and QGi,t are the active power and reactive power of the thermal power generator i in period t, respectively; Pdi,t and Qdi,t are the active power load demand and reactive load demand at bus i in period t, respectively; RUi and RDi are the climbing rate and landslide rate of thermal power units i; Pij,tl is the line power of line i,j in period t; PGimin and PGimax are the lower and upper limits of PGi; QGimin and QGimax are the lower and upper limits of QGi,t; ViminandVimax are the lower and upper limits of the voltage Vi,t at bus i; and Pijl,min and Pijl,max are the lower and upper limits of Pij,tl.


## 2.2. Constraints of the Natural Gas System

Generally, a natural gas system is composed of gas source points, supply pipes, compressors, load, and convert primary energy into electric energy through the gas turbines. If the natural gas system operates steadily, the steady-state natural gas flow through the pipeline can be assumed as a function of the upstream and downstream pressure. When the natural gas flow from nodei to node j, the gas system operational constraints include natural gas flow constraints of natural gas pipeline in (12)–(25), nodal gas balance in (26), and operating limits in (27)–(29).


### 2.2.1. Natural Gas Flow Constraints of Natural Gas Pipeline

(12)fi,j,t−=12fi,j,t−fj,i,t,(13)fi,j,t2−Ki,j2=pi,t2−pj,t2,fi,j,t>0.Otherwise,(14)fi,j,t2−Ki,j2=pj,t2−pi,t2,fi,j,t<0,where fi,j,t is the natural gas flow in the pipeline i,j in period t; f¯i,j,t is the average natural gas flow in the pipeline i,j in period t; Ki,j is a pipe constant, which is related to pipe length, diameter, friction coefficient, and gas compression coefficient; pi,t is the gas pressure at the node i in period t. Noted that equations (13) and (14) can exist precisely in the form of a SOC relaxation; thus, they can be expressed as(15)fi,j,t2−Ki,j2≤pi,t2−pj,t2,fi,j,t>0,(16)fi,j,t2−Ki,j2≥pi,t2−pj,t2,fi,j,t<0.A standard technique to convexity such a model is to relax (15), which can be solved rapidly and accurately by commercial solvers. However, the solutions that are obtained from such relaxation may yield nontrivial violations of the constraint set (13).To make the SOC model tractable while obtaining high-quality solutions, Chen et al. [28] proposed an enhanced SOC relaxation for the convex relaxation (15). Defining two sets of variables, ai,j,t and bi,j,t, which are defined as the sum and difference of the pressures at both ends of each pipe,(17)ai,j,t=pi,t+pj,t,(18)bi,j,t=pi,t−pj,t.Then, defining two new sets of auxiliary variables,γi,j,t and λi,j,t, and the convex relaxation of (16) is given by(19)γi,j,tKi,j2≥λi,j,t,(20)γi,j,t≥fi,j,t2−,(21)γi,j,t≤fi,jmax+fi,jminf¯i,j,t−fi,jmaxfi,jmin,(22)λi,j,t≥ai,jminbi,j,t+bi,jminai,j,t−ai,jminbi,jmin,(23)λi,j,t≥ai,jmaxbi,j,t+bi,jmaxai,j,t−ai,jmaxbi,jmin,(24)λi,j,t≤ai,jminbi,j,t+bi,jmaxai,j,t−ai,jminbi,jmax,(25)λi,j,t≤ai,jmaxbi,j,t+bi,jminai,j,t−ai,jmaxbi,jmin.


### 2.2.2. Gas Nodal Gas Flow Balance

(26)wgi,t+wi,tP2G−wli,t−fgi,t=fij,t.


### 2.2.3. Operating Limits

(27)fi,jmin≤fi,j,t≤fi,j,tmax,(28)pimin≤pi,t≤pimax,(29)wgimin≤wgi,t≤wgimax,where wgi,t is the injection flow at the ith gas source point in period t, wli,t is the load consumption flow at the node i in period t; wi,tP2G is the gas injection flow of the PtG equipment i in period t; fgi,t is the natural gas consumption flow of the gas turbine i connected to the power grid in period t; pimin and pimax are the lower and upper limits of pi,t; wiminandwimax are the lower and upper limits of wgi,t; fi,jmin and fi,jmax are the lower and upper limits of fi,j,t; ai,jminandai,jmin are the lower and upper limits of ai,j,t; and bi,jminandbi,jmin are the lower and upper limits of bi,j,t.


## 2.2.1. Natural Gas Flow Constraints of Natural Gas Pipeline

(12)fi,j,t−=12fi,j,t−fj,i,t,(13)fi,j,t2−Ki,j2=pi,t2−pj,t2,fi,j,t>0.Otherwise,(14)fi,j,t2−Ki,j2=pj,t2−pi,t2,fi,j,t<0,where fi,j,t is the natural gas flow in the pipeline i,j in period t; f¯i,j,t is the average natural gas flow in the pipeline i,j in period t; Ki,j is a pipe constant, which is related to pipe length, diameter, friction coefficient, and gas compression coefficient; pi,t is the gas pressure at the node i in period t. Noted that equations (13) and (14) can exist precisely in the form of a SOC relaxation; thus, they can be expressed as(15)fi,j,t2−Ki,j2≤pi,t2−pj,t2,fi,j,t>0,(16)fi,j,t2−Ki,j2≥pi,t2−pj,t2,fi,j,t<0.A standard technique to convexity such a model is to relax (15), which can be solved rapidly and accurately by commercial solvers. However, the solutions that are obtained from such relaxation may yield nontrivial violations of the constraint set (13).To make the SOC model tractable while obtaining high-quality solutions, Chen et al. [28] proposed an enhanced SOC relaxation for the convex relaxation (15). Defining two sets of variables, ai,j,t and bi,j,t, which are defined as the sum and difference of the pressures at both ends of each pipe,(17)ai,j,t=pi,t+pj,t,(18)bi,j,t=pi,t−pj,t.Then, defining two new sets of auxiliary variables,γi,j,t and λi,j,t, and the convex relaxation of (16) is given by(19)γi,j,tKi,j2≥λi,j,t,(20)γi,j,t≥fi,j,t2−,(21)γi,j,t≤fi,jmax+fi,jminf¯i,j,t−fi,jmaxfi,jmin,(22)λi,j,t≥ai,jminbi,j,t+bi,jminai,j,t−ai,jminbi,jmin,(23)λi,j,t≥ai,jmaxbi,j,t+bi,jmaxai,j,t−ai,jmaxbi,jmin,(24)λi,j,t≤ai,jminbi,j,t+bi,jmaxai,j,t−ai,jminbi,jmax,(25)λi,j,t≤ai,jmaxbi,j,t+bi,jminai,j,t−ai,jmaxbi,jmin.


## 2.2.2. Gas Nodal Gas Flow Balance

(26)wgi,t+wi,tP2G−wli,t−fgi,t=fij,t.


## 2.2.3. Operating Limits

(27)fi,jmin≤fi,j,t≤fi,j,tmax,(28)pimin≤pi,t≤pimax,(29)wgimin≤wgi,t≤wgimax,where wgi,t is the injection flow at the ith gas source point in period t, wli,t is the load consumption flow at the node i in period t; wi,tP2G is the gas injection flow of the PtG equipment i in period t; fgi,t is the natural gas consumption flow of the gas turbine i connected to the power grid in period t; pimin and pimax are the lower and upper limits of pi,t; wiminandwimax are the lower and upper limits of wgi,t; fi,jmin and fi,jmax are the lower and upper limits of fi,j,t; ai,jminandai,jmin are the lower and upper limits of ai,j,t; and bi,jminandbi,jmin are the lower and upper limits of bi,j,t.


## 2.3. Constraints of Coupling Equipment

In this section, gas turbines and PtG equipment are studied as the coupling link of the gas-electric integrated energy system.


### 2.3.1. Gas Turbines

Assuming that the gas turbine is connected between the power system nodei and the node j of the natural gas system, the conversion relationship between natural gas input and electric power output is as follows:(30)fgi,t=k1Pgi,t2+k2Pgi,t+k3,where k1, k2, and k3 are the consumption flow coefficients of the gas turbine.


### 2.3.2. PtG Equipment

(31)wi,tP2G=ηP2GPi,tP2G,where ηP2G is the conversion efficiency of PtG equipment.Note that the proposed model mentioned above is a deterministic optimization (DO) model. However, the uncertainties of wind power output have not been considered at all. It is possible that the solution is not feasible in the actual operation. Therefore, the robust stochastic optimization (RSO) method will be used to solve the proposed optimization problem with uncertainties in Section3.


## 2.3.1. Gas Turbines

Assuming that the gas turbine is connected between the power system nodei and the node j of the natural gas system, the conversion relationship between natural gas input and electric power output is as follows:(30)fgi,t=k1Pgi,t2+k2Pgi,t+k3,where k1, k2, and k3 are the consumption flow coefficients of the gas turbine.


## 2.3.2. PtG Equipment

(31)wi,tP2G=ηP2GPi,tP2G,where ηP2G is the conversion efficiency of PtG equipment.Note that the proposed model mentioned above is a deterministic optimization (DO) model. However, the uncertainties of wind power output have not been considered at all. It is possible that the solution is not feasible in the actual operation. Therefore, the robust stochastic optimization (RSO) method will be used to solve the proposed optimization problem with uncertainties in Section3.


## 3. Solving Wind Power Output Uncertainty Based on the RSO Model

### 3.1. RSO Model

The RSO model [25] for solving wind power output uncertainty in this section can be extended from a DRO model. For a multistage optimization problem, the DRO model can be formulated as(32)minf0x,ξ0+maxPξ∈ϒΕPfx,ξ,s.t.hx,ξ≤0,∀ξ∈UΓ,where x is the set of decision variables, fx,ξ and hx,ξ are the functions of x and ξ; UΓ is the uncertainty set, Γ is the adjust scale of the uncertainty set; and Pξ represents the set of all probability distributions on ℝ. Different from the DRO model, a general RSO model is expressed as follows:(33)minf0x,ξ0+maxPξ∈ϒΕPfx,ξ,s.t.hx,ξ≤0,∀ξ∈UΓ,fx,ξ∈ΛK,I.The difference between RSO and the DRO model is the existence of the constraintfx,ξ∈ΛK,I, where Λ is the eventwise recourse adaptation. And the definition of the corresponding eventwise ambiguity set can be representable in the format [25].(34)Ω=Ρ∈ρ0RIz×Sz∼,s∼∼PEPz∼s∼∈εk∈ℚk∀k∈KPz∼∈Zss∼=s=1∀s∈SPs∼=s=ps∃p∈ρ,where, for events εk, a closed convex set Zs, s∈S, ρ⊆p∈ℝ++s∑s∈Sps=1. Random variables s˜ represent a group of random scenarios with uncertain probability. For different scenarios, the support set of random variables z˜ may be different, and under the condition of event implementation, the expectations of z˜ may also be different.


### 3.2. Wasserstein Distance-Based Ambiguity Set

The construction of the ambiguity set has an important influence on the expectation value of adjustment cost. And ambiguity set based on Wasserstein distance (WD) can still reflect the distance between two probability distributions when the support sets of two probability distribution functions do not overlap or overlap very little and can measure the “distance” between any two probability distributions more accurately, which is a more reasonable ambiguity set.For anyp∈1,∞, the type-p Wasserstein metric between two distributions Q and Q1 is defined as(35)WpQ,Q1=infπ∈ΠQ,Q1∫w1−w2pΠdw1,dw21/p,where WpQ,Q1 is the Wasserstein distance between Q and Q1; Π is the joint probability distribution of the sum of random variables; and Q and Q1 are the marginal distributions of w1 and w2.Considering a data-driven setting, such as the Wasserstein ambiguity sets centered on empirical distributionP1=1/s∑s∈sδus, a manageable distance measure ρ:ℝIu×ℝIu↦0,∞ is given, and the p−norm Wasserstein metric between two distributions P and P1 is defined by an optimization problem(36)dwPP,P1≜infℚ∈QP,P1Eℚρu,∼u∼+ifp∈1,∞,infℚ∈QP,P1ℚ−esssupρu,∼u∼+ifp=∞,where u˜∼P,u˜+∼P1; and QP,P1 is the set of all joint probability distributions with marginal distributions P and P1.Furthermore, the ambiguity set based on Wasserstein distance has the following forms:(37)ℜwpθ=P∈ρ0Uu∼∼PdwP,P1≤θp,where θ is the radius constant of Wasserstein. Thus, a new representation of the 1-norm Wasserstein ambiguity set can be provided in the form of an eventwise ambiguity set:(38)Ωw1θ=Ρ∈ρ0RIz×Su∼,v∼,s∼∼PEPv∼s∼∈S≤θPu∼,v∼∈Zss∼=s=1∀s∈SPs∼=s=1s∀s∈S.In equation (38), for different scenes s∈S, random variables u˜ and auxiliary random variables v˜ belong to the support set Zs=u,vu∈U,v≥ρu,u˜s.Similarly, the new representation of the∞−norm Wasserstein ambiguity set can be provided in the form of an eventwise ambiguity set(39)Ωw∞θ=Ρ∈ρ0RIz×Su∼,s∼∼PPu∼∈U,ρu∼,u∼s≤θs∼=s=1Ps∼=s=1s∀s∈S.Obviously, we can effectively determine the worst-case expectation over the eventwise ambiguity setΩw by solving a classical robust optimization problem. Assuming Slater’s condition holds the worst-case expectation:(40)supP∈ΩEPrTs∼Gms∼z∼+hms∼.It is equivalent to the optimal value of the following classical robust optimization problem:(41)infϒs.t.ϒ≥αTp+∑k∈KβkTuk∀p∈ρ,uk∑s∈εkps∈ℚk,αs+∑k∈KβkTz≥rTsGmsz+hms∀z∈Zs,s∈S,ϒ∈ℝ,α∈ℝs,βk∈ℝIz∀k∈K.The computability of traditional robust optimization problems depends on the uncertainty set. Using the algebraic modeling toolbox such as RSOME [29], the robust optimization problem can be automatically transformed into a polynomial-size linear or second-order cone optimization problem, which can be solved by commercial solvers such as CPLEX, GUROBI, and MOSEK.


## 3.1. RSO Model

The RSO model [25] for solving wind power output uncertainty in this section can be extended from a DRO model. For a multistage optimization problem, the DRO model can be formulated as(32)minf0x,ξ0+maxPξ∈ϒΕPfx,ξ,s.t.hx,ξ≤0,∀ξ∈UΓ,where x is the set of decision variables, fx,ξ and hx,ξ are the functions of x and ξ; UΓ is the uncertainty set, Γ is the adjust scale of the uncertainty set; and Pξ represents the set of all probability distributions on ℝ. Different from the DRO model, a general RSO model is expressed as follows:(33)minf0x,ξ0+maxPξ∈ϒΕPfx,ξ,s.t.hx,ξ≤0,∀ξ∈UΓ,fx,ξ∈ΛK,I.The difference between RSO and the DRO model is the existence of the constraintfx,ξ∈ΛK,I, where Λ is the eventwise recourse adaptation. And the definition of the corresponding eventwise ambiguity set can be representable in the format [25].(34)Ω=Ρ∈ρ0RIz×Sz∼,s∼∼PEPz∼s∼∈εk∈ℚk∀k∈KPz∼∈Zss∼=s=1∀s∈SPs∼=s=ps∃p∈ρ,where, for events εk, a closed convex set Zs, s∈S, ρ⊆p∈ℝ++s∑s∈Sps=1. Random variables s˜ represent a group of random scenarios with uncertain probability. For different scenarios, the support set of random variables z˜ may be different, and under the condition of event implementation, the expectations of z˜ may also be different.


## 3.2. Wasserstein Distance-Based Ambiguity Set

The construction of the ambiguity set has an important influence on the expectation value of adjustment cost. And ambiguity set based on Wasserstein distance (WD) can still reflect the distance between two probability distributions when the support sets of two probability distribution functions do not overlap or overlap very little and can measure the “distance” between any two probability distributions more accurately, which is a more reasonable ambiguity set.For anyp∈1,∞, the type-p Wasserstein metric between two distributions Q and Q1 is defined as(35)WpQ,Q1=infπ∈ΠQ,Q1∫w1−w2pΠdw1,dw21/p,where WpQ,Q1 is the Wasserstein distance between Q and Q1; Π is the joint probability distribution of the sum of random variables; and Q and Q1 are the marginal distributions of w1 and w2.Considering a data-driven setting, such as the Wasserstein ambiguity sets centered on empirical distributionP1=1/s∑s∈sδus, a manageable distance measure ρ:ℝIu×ℝIu↦0,∞ is given, and the p−norm Wasserstein metric between two distributions P and P1 is defined by an optimization problem(36)dwPP,P1≜infℚ∈QP,P1Eℚρu,∼u∼+ifp∈1,∞,infℚ∈QP,P1ℚ−esssupρu,∼u∼+ifp=∞,where u˜∼P,u˜+∼P1; and QP,P1 is the set of all joint probability distributions with marginal distributions P and P1.Furthermore, the ambiguity set based on Wasserstein distance has the following forms:(37)ℜwpθ=P∈ρ0Uu∼∼PdwP,P1≤θp,where θ is the radius constant of Wasserstein. Thus, a new representation of the 1-norm Wasserstein ambiguity set can be provided in the form of an eventwise ambiguity set:(38)Ωw1θ=Ρ∈ρ0RIz×Su∼,v∼,s∼∼PEPv∼s∼∈S≤θPu∼,v∼∈Zss∼=s=1∀s∈SPs∼=s=1s∀s∈S.In equation (38), for different scenes s∈S, random variables u˜ and auxiliary random variables v˜ belong to the support set Zs=u,vu∈U,v≥ρu,u˜s.Similarly, the new representation of the∞−norm Wasserstein ambiguity set can be provided in the form of an eventwise ambiguity set(39)Ωw∞θ=Ρ∈ρ0RIz×Su∼,s∼∼PPu∼∈U,ρu∼,u∼s≤θs∼=s=1Ps∼=s=1s∀s∈S.Obviously, we can effectively determine the worst-case expectation over the eventwise ambiguity setΩw by solving a classical robust optimization problem. Assuming Slater’s condition holds the worst-case expectation:(40)supP∈ΩEPrTs∼Gms∼z∼+hms∼.It is equivalent to the optimal value of the following classical robust optimization problem:(41)infϒs.t.ϒ≥αTp+∑k∈KβkTuk∀p∈ρ,uk∑s∈εkps∈ℚk,αs+∑k∈KβkTz≥rTsGmsz+hms∀z∈Zs,s∈S,ϒ∈ℝ,α∈ℝs,βk∈ℝIz∀k∈K.The computability of traditional robust optimization problems depends on the uncertainty set. Using the algebraic modeling toolbox such as RSOME [29], the robust optimization problem can be automatically transformed into a polynomial-size linear or second-order cone optimization problem, which can be solved by commercial solvers such as CPLEX, GUROBI, and MOSEK.


## 4. Numerical Simulation

### 4.1. Simulation System

The proposed model is tested on the electricity-gas IES, which connects wind power generators. In order to verify the accuracy and validity of the proposed model, two test systems are used for the simulation. The first test system is integrated IEEE 4 and 6-node gas [30] systems that connect one wind power generator, denoted as the EG-4 system, and the second system is integrated IEEE 118 and the Belgium 20-node gas [31] systems that connect six distributed wind power generators, denoted as the EG 118 system. The remaining detailed information of the IEEE 4-bus system and the IEEE 118-bus system can be obtained from the MATPOWER library [32]. All these tests are implemented in a PC (Intel i7-7700 Core CPU, 3.6-GHz, 8 GB RAM). The proposed model is carried on MATLAB R2018a by using RSOME package as modeling software.


### 4.2. EG-4 System

The structure of the EG-4 system is shown in Figure1. Besides, the capacity of each wind generator and gas turbines is 100 MW and 150 MW, respectively. ηP2G is 0.8; k1,k2,andk3 shown in (30) is 0, 0.05, and 0, respectively, in the simulation. The power load, natural gas load, and wind power predicted output curves of the electricity-gas IES are shown in Figure 2.Figure 1
The structure of the EG-4 system.Figure 2
The curves of forecasted wind power output, power load, and natural gas load.In order to study the impact of the penetration of wind power generation, adjusting cost coefficient of wind unit and the sample size on the operation of the electricity-gas IES, the following two cases are set for comparative analysis:Case 1.
Constructed the 1-norm Wasserstein ambiguity set to deal with the uncertainty of wind power output.Case 2.
Constructed the∞−norm Wasserstein ambiguity set to deal with the uncertainty of wind power output.


#### 4.2.1. Influence of Different Wind Power Penetration Rates on the Results of the RSO Model

The challenge of wind power uncertainty to the electricity-gas IES is usually directly related to wind power penetration rate, and using different ambiguity sets to deal with the uncertainty of wind power output will obtain different results. Thus, in order to study the impact of the penetration rate of wind power generation on the operation of the electricity-gas IES, the operation results under different penetration rates in the two cases are calculated, respectively, which are shown in Table1 and Figures 3–6.Table 1
Comparison of objective function value in two cases under different wind power penetration.
Penetration rate of wind power (%)Case1Case2Objective function value ($)Percentage decrease in objective function valueObjective function value ($)Percentage decrease in objective function value15.610.21 × 104—8.28 × 104—21.79.86 × 1043.438.04 × 1042.9027.09.59 × 1046.077.72 × 1046.76Figure 3
Comparison of actual wind output in two cases under different wind penetration rates.Figure 4
Comparison of the output of the PtG unit in two cases under different wind penetration rates.Figure 5
Comparison of thermal power units actual output under different wind penetration rates in Case1.Figure 6
Comparison of thermal power units actual output under different wind penetration rates in Case2.As presented in Table1, with the increase of wind power penetration rate, the economy of the whole electricity-gas IES is better in both cases. The relationship between penetration rate and optimal value is not linear because it is related to the safe and stable operation of the system and has many influencing factors that will affect optimal value.According to Figures3 and 4, it shows that with the increase of wind power penetration rate, wind power actual output is increasing in both cases. And wind power output of Case 2 is higher than that of Case 1 in the same penetration rate. When the wind power output increases, the output of PtG units also increases significantly. This is because when wind power is more injected to the grid, PtG units can be used to absorb a part of the remaining wind power and improve the utilization rate of wind power.From Figures5 and 6, it shows that the relationship between the output of thermal power unit PG4 and wind power penetration rate is a positive correlation, and the shape of its output curve is similar to that of power load curve. Moreover, the output of PG4 is greater than that of PG3, which indicates that PG4 plays a role of peak load unit in the whole system operation.


#### 4.2.2. Influence of Different Adjusting Cost Coefficients of Wind Unit on the Results of the RSO Model

In order to compare and analyze the impact of wind power adjustment cost coefficient on the wind power utilization process of the system, the adjustment cost coefficient of wind power is taken as 0, 0.5, and 1, respectively, for calculation. The corresponding wind power utilization rate and objective function values are shown in Table2.Table 2
Comparison of wind power utilization in two cases under different wind adjusting cost coefficientck,twind.
Adjusting cost coefficientck,twind (p.u.)Case1Case2Objective function value ($)Wind power utilizationObjective function value ($)Wind power utilization09 × 10494.93%——0.58.98 × 10499.91%7.71 × 10499.43%18.89 × 10499.95%7.24 × 10499.54%As shown in Table2, the utilization ratio of wind power is 94.93% without considering the adjustment cost coefficient of wind power in Case 1. When the coefficient ck,twind is 0.5, the utilization rate of wind power has reached more than 99% in both cases. Besides, when the adjustment cost coefficient ranged from 0 to 1, the objective value is decreasing with the increasing adjustment cost coefficient.


#### 4.2.3. Influence of Different Sample Sizes on the Results of the RSO Model

In order to verify the robustness and economy of the proposed model, the objective function values under different sample sizes are compared with the distribution of operating costs obtained by the Monte Carlo method (MCM). Among them, SO is to extract 2000 real scenery forecast scenes and put them into the model to calculate the operation cost.Figure7 shows the distribution of the objective function values of the model and the results obtained by the SO method under different sample sizes. Under the same sample size, the objective function value of Case 1 is higher than that of Case 2 and SO. And when the sample size changes from 1000 to 2000, the objective function values tend to be consistent. This shows that the proposed model is less dependent on the sample size and tends to the final result when the sample value is around 2000, which improves the computation efficiency and avoids the large error caused by the small sample size. As shown in the figure, no matter what the sample size is, the total generation costs calculated from the real scenario is always below the objective function value of the model. It shows that decision-making based on the model can deal with the scene of the actual operation of the system and has robustness. As the sample size increases, the distance between them decreases. This is because with the rise of sample size, the ambiguity uncertainty set is shrinking, and the extreme probability distribution tends to the real probability distribution. It also shows that the more the samples are obtained, the lower the conservative degree of the model solution is and the better the economy is.Figure 7
Comparison of objective function value by SO, Case1, and Case 2 with different historical data.Table3 shows the voltage level of the system in a certain period t = 8 when the sample size is 2000. It can be seen that the voltage level of the system is within the allowable range. Similarly, other sample sizes are similar and will not be shown.Table 3
The voltage level of the system atT = 8 when the sample size is 2000.
Bus voltage amplitude (p.u.)Case1Case2Bus 111Bus 21.0201.006Bus 31.0941.096Bus 41.0911.094To sum up, compared with Case1, Case 2 can get more economic optimal solution under the condition of safe operation of the electricity-gas IES. That is, the ambiguity set constructed by ∞−norm Wasserstein is better than that based on 1-norm Wasserstein.


### 4.3. EG-118 System

To further illustrate the effectiveness of the proposed model, a large system is analyzed based on the modified IEEE 118 node power system and Belgium 20-node natural gas system, which is denoted as the EG-118 system.In the EG-118 system, ten thermal units are connected to the power system nodes at 59, 61, 65, 66, 69, 80, 89, 100, 103, and 111, and the capacity of these generators are 200 MW; three gas turbines are connected to the power system buses at 10, 25, and 26, and these generators are connected to nodes 3, 4, and 7 of the natural gas networks, respectively; two PtG units are connected to the power system buses at 49 and 54, and these generators are connected to nodes 9 and 10 of the natural gas networks, respectively. Constructing the∞−norm Wasserstein ambiguity set to deal with the uncertainty of wind power output, when the wind power penetration rate is 15% and the sample size is 1000, the objective function value is $6.38 × 108, and output of the thermal units PG59, PG61, PG69, PG80, and PG111 are shown in Figure 8. The voltage level of the system in a certain period t = 8 is shown in Figure 9.Figure 8
The output of thermal units.Figure 9
The bus voltage level of the system in a certain periodT = 8.Figure8 shows that the thermal unit may play a role in the peak load unit in the whole system operation. And from Figure 9, it can be seen that the voltage level of the system is within the allowable range, which proves that the proposed model is feasible in practical operation.


## 4.1. Simulation System

The proposed model is tested on the electricity-gas IES, which connects wind power generators. In order to verify the accuracy and validity of the proposed model, two test systems are used for the simulation. The first test system is integrated IEEE 4 and 6-node gas [30] systems that connect one wind power generator, denoted as the EG-4 system, and the second system is integrated IEEE 118 and the Belgium 20-node gas [31] systems that connect six distributed wind power generators, denoted as the EG 118 system. The remaining detailed information of the IEEE 4-bus system and the IEEE 118-bus system can be obtained from the MATPOWER library [32]. All these tests are implemented in a PC (Intel i7-7700 Core CPU, 3.6-GHz, 8 GB RAM). The proposed model is carried on MATLAB R2018a by using RSOME package as modeling software.


## 4.2. EG-4 System

The structure of the EG-4 system is shown in Figure1. Besides, the capacity of each wind generator and gas turbines is 100 MW and 150 MW, respectively. ηP2G is 0.8; k1,k2,andk3 shown in (30) is 0, 0.05, and 0, respectively, in the simulation. The power load, natural gas load, and wind power predicted output curves of the electricity-gas IES are shown in Figure 2.Figure 1
The structure of the EG-4 system.Figure 2
The curves of forecasted wind power output, power load, and natural gas load.In order to study the impact of the penetration of wind power generation, adjusting cost coefficient of wind unit and the sample size on the operation of the electricity-gas IES, the following two cases are set for comparative analysis:Case 1.
Constructed the 1-norm Wasserstein ambiguity set to deal with the uncertainty of wind power output.Case 2.
Constructed the∞−norm Wasserstein ambiguity set to deal with the uncertainty of wind power output.


### 4.2.1. Influence of Different Wind Power Penetration Rates on the Results of the RSO Model

The challenge of wind power uncertainty to the electricity-gas IES is usually directly related to wind power penetration rate, and using different ambiguity sets to deal with the uncertainty of wind power output will obtain different results. Thus, in order to study the impact of the penetration rate of wind power generation on the operation of the electricity-gas IES, the operation results under different penetration rates in the two cases are calculated, respectively, which are shown in Table1 and Figures 3–6.Table 1
Comparison of objective function value in two cases under different wind power penetration.
Penetration rate of wind power (%)Case1Case2Objective function value ($)Percentage decrease in objective function valueObjective function value ($)Percentage decrease in objective function value15.610.21 × 104—8.28 × 104—21.79.86 × 1043.438.04 × 1042.9027.09.59 × 1046.077.72 × 1046.76Figure 3
Comparison of actual wind output in two cases under different wind penetration rates.Figure 4
Comparison of the output of the PtG unit in two cases under different wind penetration rates.Figure 5
Comparison of thermal power units actual output under different wind penetration rates in Case1.Figure 6
Comparison of thermal power units actual output under different wind penetration rates in Case2.As presented in Table1, with the increase of wind power penetration rate, the economy of the whole electricity-gas IES is better in both cases. The relationship between penetration rate and optimal value is not linear because it is related to the safe and stable operation of the system and has many influencing factors that will affect optimal value.According to Figures3 and 4, it shows that with the increase of wind power penetration rate, wind power actual output is increasing in both cases. And wind power output of Case 2 is higher than that of Case 1 in the same penetration rate. When the wind power output increases, the output of PtG units also increases significantly. This is because when wind power is more injected to the grid, PtG units can be used to absorb a part of the remaining wind power and improve the utilization rate of wind power.From Figures5 and 6, it shows that the relationship between the output of thermal power unit PG4 and wind power penetration rate is a positive correlation, and the shape of its output curve is similar to that of power load curve. Moreover, the output of PG4 is greater than that of PG3, which indicates that PG4 plays a role of peak load unit in the whole system operation.


### 4.2.2. Influence of Different Adjusting Cost Coefficients of Wind Unit on the Results of the RSO Model

In order to compare and analyze the impact of wind power adjustment cost coefficient on the wind power utilization process of the system, the adjustment cost coefficient of wind power is taken as 0, 0.5, and 1, respectively, for calculation. The corresponding wind power utilization rate and objective function values are shown in Table2.Table 2
Comparison of wind power utilization in two cases under different wind adjusting cost coefficientck,twind.
Adjusting cost coefficientck,twind (p.u.)Case1Case2Objective function value ($)Wind power utilizationObjective function value ($)Wind power utilization09 × 10494.93%——0.58.98 × 10499.91%7.71 × 10499.43%18.89 × 10499.95%7.24 × 10499.54%As shown in Table2, the utilization ratio of wind power is 94.93% without considering the adjustment cost coefficient of wind power in Case 1. When the coefficient ck,twind is 0.5, the utilization rate of wind power has reached more than 99% in both cases. Besides, when the adjustment cost coefficient ranged from 0 to 1, the objective value is decreasing with the increasing adjustment cost coefficient.


### 4.2.3. Influence of Different Sample Sizes on the Results of the RSO Model

In order to verify the robustness and economy of the proposed model, the objective function values under different sample sizes are compared with the distribution of operating costs obtained by the Monte Carlo method (MCM). Among them, SO is to extract 2000 real scenery forecast scenes and put them into the model to calculate the operation cost.Figure7 shows the distribution of the objective function values of the model and the results obtained by the SO method under different sample sizes. Under the same sample size, the objective function value of Case 1 is higher than that of Case 2 and SO. And when the sample size changes from 1000 to 2000, the objective function values tend to be consistent. This shows that the proposed model is less dependent on the sample size and tends to the final result when the sample value is around 2000, which improves the computation efficiency and avoids the large error caused by the small sample size. As shown in the figure, no matter what the sample size is, the total generation costs calculated from the real scenario is always below the objective function value of the model. It shows that decision-making based on the model can deal with the scene of the actual operation of the system and has robustness. As the sample size increases, the distance between them decreases. This is because with the rise of sample size, the ambiguity uncertainty set is shrinking, and the extreme probability distribution tends to the real probability distribution. It also shows that the more the samples are obtained, the lower the conservative degree of the model solution is and the better the economy is.Figure 7
Comparison of objective function value by SO, Case1, and Case 2 with different historical data.Table3 shows the voltage level of the system in a certain period t = 8 when the sample size is 2000. It can be seen that the voltage level of the system is within the allowable range. Similarly, other sample sizes are similar and will not be shown.Table 3
The voltage level of the system atT = 8 when the sample size is 2000.
Bus voltage amplitude (p.u.)Case1Case2Bus 111Bus 21.0201.006Bus 31.0941.096Bus 41.0911.094To sum up, compared with Case1, Case 2 can get more economic optimal solution under the condition of safe operation of the electricity-gas IES. That is, the ambiguity set constructed by ∞−norm Wasserstein is better than that based on 1-norm Wasserstein.


## 4.2.1. Influence of Different Wind Power Penetration Rates on the Results of the RSO Model

The challenge of wind power uncertainty to the electricity-gas IES is usually directly related to wind power penetration rate, and using different ambiguity sets to deal with the uncertainty of wind power output will obtain different results. Thus, in order to study the impact of the penetration rate of wind power generation on the operation of the electricity-gas IES, the operation results under different penetration rates in the two cases are calculated, respectively, which are shown in Table1 and Figures 3–6.Table 1
Comparison of objective function value in two cases under different wind power penetration.
Penetration rate of wind power (%)Case1Case2Objective function value ($)Percentage decrease in objective function valueObjective function value ($)Percentage decrease in objective function value15.610.21 × 104—8.28 × 104—21.79.86 × 1043.438.04 × 1042.9027.09.59 × 1046.077.72 × 1046.76Figure 3
Comparison of actual wind output in two cases under different wind penetration rates.Figure 4
Comparison of the output of the PtG unit in two cases under different wind penetration rates.Figure 5
Comparison of thermal power units actual output under different wind penetration rates in Case1.Figure 6
Comparison of thermal power units actual output under different wind penetration rates in Case2.As presented in Table1, with the increase of wind power penetration rate, the economy of the whole electricity-gas IES is better in both cases. The relationship between penetration rate and optimal value is not linear because it is related to the safe and stable operation of the system and has many influencing factors that will affect optimal value.According to Figures3 and 4, it shows that with the increase of wind power penetration rate, wind power actual output is increasing in both cases. And wind power output of Case 2 is higher than that of Case 1 in the same penetration rate. When the wind power output increases, the output of PtG units also increases significantly. This is because when wind power is more injected to the grid, PtG units can be used to absorb a part of the remaining wind power and improve the utilization rate of wind power.From Figures5 and 6, it shows that the relationship between the output of thermal power unit PG4 and wind power penetration rate is a positive correlation, and the shape of its output curve is similar to that of power load curve. Moreover, the output of PG4 is greater than that of PG3, which indicates that PG4 plays a role of peak load unit in the whole system operation.


## 4.2.2. Influence of Different Adjusting Cost Coefficients of Wind Unit on the Results of the RSO Model

In order to compare and analyze the impact of wind power adjustment cost coefficient on the wind power utilization process of the system, the adjustment cost coefficient of wind power is taken as 0, 0.5, and 1, respectively, for calculation. The corresponding wind power utilization rate and objective function values are shown in Table2.Table 2
Comparison of wind power utilization in two cases under different wind adjusting cost coefficientck,twind.
Adjusting cost coefficientck,twind (p.u.)Case1Case2Objective function value ($)Wind power utilizationObjective function value ($)Wind power utilization09 × 10494.93%——0.58.98 × 10499.91%7.71 × 10499.43%18.89 × 10499.95%7.24 × 10499.54%As shown in Table2, the utilization ratio of wind power is 94.93% without considering the adjustment cost coefficient of wind power in Case 1. When the coefficient ck,twind is 0.5, the utilization rate of wind power has reached more than 99% in both cases. Besides, when the adjustment cost coefficient ranged from 0 to 1, the objective value is decreasing with the increasing adjustment cost coefficient.


## 4.2.3. Influence of Different Sample Sizes on the Results of the RSO Model

In order to verify the robustness and economy of the proposed model, the objective function values under different sample sizes are compared with the distribution of operating costs obtained by the Monte Carlo method (MCM). Among them, SO is to extract 2000 real scenery forecast scenes and put them into the model to calculate the operation cost.Figure7 shows the distribution of the objective function values of the model and the results obtained by the SO method under different sample sizes. Under the same sample size, the objective function value of Case 1 is higher than that of Case 2 and SO. And when the sample size changes from 1000 to 2000, the objective function values tend to be consistent. This shows that the proposed model is less dependent on the sample size and tends to the final result when the sample value is around 2000, which improves the computation efficiency and avoids the large error caused by the small sample size. As shown in the figure, no matter what the sample size is, the total generation costs calculated from the real scenario is always below the objective function value of the model. It shows that decision-making based on the model can deal with the scene of the actual operation of the system and has robustness. As the sample size increases, the distance between them decreases. This is because with the rise of sample size, the ambiguity uncertainty set is shrinking, and the extreme probability distribution tends to the real probability distribution. It also shows that the more the samples are obtained, the lower the conservative degree of the model solution is and the better the economy is.Figure 7
Comparison of objective function value by SO, Case1, and Case 2 with different historical data.Table3 shows the voltage level of the system in a certain period t = 8 when the sample size is 2000. It can be seen that the voltage level of the system is within the allowable range. Similarly, other sample sizes are similar and will not be shown.Table 3
The voltage level of the system atT = 8 when the sample size is 2000.
Bus voltage amplitude (p.u.)Case1Case2Bus 111Bus 21.0201.006Bus 31.0941.096Bus 41.0911.094To sum up, compared with Case1, Case 2 can get more economic optimal solution under the condition of safe operation of the electricity-gas IES. That is, the ambiguity set constructed by ∞−norm Wasserstein is better than that based on 1-norm Wasserstein.


## 4.3. EG-118 System

To further illustrate the effectiveness of the proposed model, a large system is analyzed based on the modified IEEE 118 node power system and Belgium 20-node natural gas system, which is denoted as the EG-118 system.In the EG-118 system, ten thermal units are connected to the power system nodes at 59, 61, 65, 66, 69, 80, 89, 100, 103, and 111, and the capacity of these generators are 200 MW; three gas turbines are connected to the power system buses at 10, 25, and 26, and these generators are connected to nodes 3, 4, and 7 of the natural gas networks, respectively; two PtG units are connected to the power system buses at 49 and 54, and these generators are connected to nodes 9 and 10 of the natural gas networks, respectively. Constructing the∞−norm Wasserstein ambiguity set to deal with the uncertainty of wind power output, when the wind power penetration rate is 15% and the sample size is 1000, the objective function value is $6.38 × 108, and output of the thermal units PG59, PG61, PG69, PG80, and PG111 are shown in Figure 8. The voltage level of the system in a certain period t = 8 is shown in Figure 9.Figure 8
The output of thermal units.Figure 9
The bus voltage level of the system in a certain periodT = 8.Figure8 shows that the thermal unit may play a role in the peak load unit in the whole system operation. And from Figure 9, it can be seen that the voltage level of the system is within the allowable range, which proves that the proposed model is feasible in practical operation.


## 5. Conclusion

In this paper, an RSO model is used to solve the dynamic optimal power flow model for the electricity-gas IES considering wind power uncertainty, where the ambiguity set of wind power output is constructed based on Wasserstein distance and then affined to an eventwise ambiguity set. The main conclusions are as follows:(1)
It is a feasible method to describe the uncertainty of scenery as Wasserstein packet with Wasserstein distance as radius and the empirical distribution of wind prediction error as the center. With the increase of historical data, the ambiguity uncertainty set decreases, and the model solution tends to be a stochastic optimization solution. It shows that the proposed model has a better performance than the traditional robust optimization. The lower degree of conservatism can greatly reduce the operation cost under the premise of ensuring the safety of the system.(2)
The AC power flow model is approximated as a decoupled linear power flow model to improve the feasibility of the generation plan obtained from the model. The simulation results show that the voltage obtained by the model meets the requirements of the system.The proposed RSO model is a single-stage optimization problem, and the future work will turn to the dynamic OPF problem of the multistage RSO model in the electricity-gas IES.



---

*Source: 8879906-2020-10-12.xml*