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

)e 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. )e traditional stochastic optimization and robust optimization algorithms have some limitations and deficiencies in dealing with the uncertainty of wind power output. )erefore, 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. )en, 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.


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][2][3][4][5][6][7][8]. Compared with other primary energies, natural gas and wind power have the characteristics of clean and environmental protection. us, 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. e 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. e 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. is 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. e 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. e 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. e 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. erefore, this paper uses RSO algorithm to calculate the influence of the 1norm 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]. e 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. e 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.
is 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. e main contributions are as follows: (1) e 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. e simulation results verify that the Wasserstein ambiguity set based on ∞−norm has a more economic optimal solution. (3) e 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. e rest of this paper is organized as follows. Section 2 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. e conclusion is presented in Section 5.

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. e 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. e 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 where T is the total number of time periods in the dispatching cycle; N p , N gas , and N w are the sets of the thermal power units, the gas supply point, and the all distributed wind, respectively; P Gi,t is the active power of the thermal power generator of bus i in period t; C gj is the price of gas supply point i; w gj,t is the injection flow at the i th gas source point in period t; a i , b i , and c i are the consumption characteristic parameters of the thermal power unit at bus i; c wind k,t is the adjusting cost coefficient of wind unit k in period t; p s is the probability of each discrete scenario s; and P w0 k,t and P w k,t are the predicted output and actual output of wind power in period t of the wind turbine k in the scenario s.

Power Balance Equation Constraint
.
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: e nonlinear term in the above formula is approximated as follows: en, (4) and (5) can be simplified as Similarly, the line power can be expressed as

Climbing Constraints
where G ij + jB ij is the element of power system admittance matrix between branches line (i, j); G ij ′ + jB ij ′ represents the element of admittance matrix without self-admittance; g ij + jb ij is the admittance of line (i, j); V i,t is the voltage at bus i in period t ; P i,t and Q i,t are the active power and reactive power injected by node i in period t, respectively; P gi,t is the active power output of the gas turbine i in period t; P PtG i,t is the active power output of the power to gas equipment i in period t; P Gi,t and Q Gi,t are the active power and reactive power of the thermal power generator i in period t, respectively; P di,t and Q di,t are the active power load demand and reactive load demand at bus i in period t, respectively;

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 node i 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).

Natural Gas Flow Constraints of Natural
Gas Pipeline. Otherwise, where f i,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; K i,j is a pipe constant, which is related to pipe length, diameter, friction coefficient, and gas compression coefficient; p i,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 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 highquality solutions, Chen et al. [28] proposed an enhanced SOC relaxation for the convex relaxation (15). Defining two sets of variables, a i,j,t and b i,j,t , which are defined as the sum and difference of the pressures at both ends of each pipe, en, defining two new sets of auxiliary variables, c i,j,t and λ i,j,t , and the convex relaxation of (16) is given by

Operating Limits
where w gi,t is the injection flow at the i th gas source point in period t, w li,t is the load consumption flow at the node i in period t; w P2G i,t is the gas injection flow of the PtG equipment i in period t; f gi,t is the natural gas consumption flow of the gas turbine i connected to the power grid in period t; p min

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.

Gas Turbines.
Assuming that the gas turbine is connected between the power system node i and the node j of the natural gas system, the conversion relationship between natural gas input and electric power output is as follows: where k1, k2, and k3 are the consumption flow coefficients of the gas turbine.

PtG Equipment.
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. erefore, the robust stochastic optimization (RSO) method will be used to solve the proposed optimization problem with uncertainties in Section 3.

Solving Wind Power Output Uncertainty
Based on the RSO Model

RSO
Model. e 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 where x is the set of decision variables, f(x, ξ) and h(x, ξ) 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 R. Different from the DRO model, a general RSO model is expressed as follows: (33) e difference between RSO and the DRO model is the existence of the constraint f(x, ξ) ∈ Λ(K, I), where Λ is the eventwise recourse adaptation. And the definition of the corresponding eventwise ambiguity set can be representable in the format [25].
where, for events ε k , a closed convex set Z s , s ∈ [S], ρ ⊆ p ∈ R s ++ | s∈[S] p s � 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.

Wasserstein Distance-Based Ambiguity Set.
e 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 any p ∈ [1, ∞), the type-p Wasserstein metric between two distributions Q and Q 1 is defined as where W p (Q, Q 1 ) is the Wasserstein distance between Q and Q 1 ; Π is the joint probability distribution of the sum of random variables; and Q and Q 1 are the marginal distributions of w 1 and w 2 . Considering a data-driven setting, such as the Wasserstein ambiguity sets centered on empirical distribution P 1 � (1/s) s∈[s] δ u s , a manageable distance measure ρ: R I u × R I u ↦[0, ∞) is given, and the p−norm Wasserstein metric between two distributions P and P 1 is defined by an optimization problem where u ∼ P, u + ∼ P 1 ; and Q(P, P 1 ) is the set of all joint probability distributions with marginal distributions P and P 1 . Furthermore, the ambiguity set based on Wasserstein distance has the following forms: where θ is the radius constant of Wasserstein. us, a new representation of the 1-norm Wasserstein ambiguity set can be provided in the form of an eventwise ambiguity set: In equation (38), for different scenes s ∈ [S], random variables u and auxiliary random variables v belong to the support set Z s � (u, v) | u ∈ U, v ≥ ρ(u, us) .
Similarly, the new representation of the ∞−norm Wasserstein ambiguity set can be provided in the form of an eventwise ambiguity set 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: It is equivalent to the optimal value of the following classical robust optimization problem: e 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.

Simulation System.
e 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. e 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. e remaining detailed information of the IEEE 4bus 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). e proposed model is carried on MATLAB R2018a by using RSOME package as modeling software.

EG-4 System.
e structure of the EG-4 system is shown in Figure 1. Besides, the capacity of each wind generator and gas turbines is 100 MW and 150 MW, respectively. η P2G is 0.8; k1, k2, and k3 shown in (30) is 0, 0.05, and 0, respectively, in the simulation. e power load, natural gas load, and wind power predicted output curves of the electricitygas IES are shown in Figure 2.
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.

Influence of Different Wind Power Penetration Rates on the Results of the RSO Model.
e 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. us, 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 Table 1 and Figures 3-6.
As presented in Table 1, with the increase of wind power penetration rate, the economy of the whole electricity-gas IES is better in both cases. e 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 Figures 3 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. is 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 Figures 5 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.

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. e corresponding wind power utilization rate and objective function values are shown in Table 2.
As shown in Table 2, the utilization ratio of wind power is 94.93% without considering the adjustment cost coefficient of wind power in Case 1. When the coefficient c wind k,t 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.        8 Complexity Figure 7 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. is 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. is 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. Table 3 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.

Influence of Different
To sum up, compared with Case 1, Case 2 can get more economic optimal solution under the condition of safe operation of the electricity-gas IES. at is, the ambiguity set constructed by ∞−norm Wasserstein is better than that based on 1-norm Wasserstein.

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    Complexity output, when the wind power penetration rate is 15% and the sample size is 1000, the objective function value is $6.38 × 10 8 , and output of the thermal units PG59, PG61, PG69, PG80, and PG111 are shown in Figure 8. e voltage level of the system in a certain period t � 8 is shown in Figure 9. Figure 8 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.

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. e 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. e lower degree of conservatism can greatly reduce the operation cost under the premise of ensuring the safety of the system. (2) e 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. e simulation results show that the voltage obtained by the model meets the requirements of the system. e 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.

Data Availability
e data used to support the findings of this study are included within this article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.  10 Complexity