Robust Day-Ahead Scheduling of Electricity and Natural Gas Systems via a Risk-Averse Adjustable Uncertainty Set Approach

Abstract

The requirement for energy sustainability drives the development of renewable energy technologies and gas-fired power generation. The increasing installation of gas-fired units significantly intensifies the interdependency between the electricity system and natural gas system. The joint scheduling of electricity and natural gas systems has become an attractive option for improving energy efficiency. This paper proposes a robust day-ahead scheduling model for electricity and natural gas system, which minimizes the total cost including fuel cost, spinning reserve cost and cost of operational risk while ensuring the feasibility for all scenarios within the uncertainty set. Different from the conventional robust optimization with predefined uncertainty set, a new approach with risk-averse adjustable uncertainty set is proposed in this paper to mitigate the conservatism. Furthermore, the Wasserstein–Moment metric is applied to construct ambiguity sets for computing operational risk. The proposed scheduling model is solved by the column-and-constraint generation method. The effectiveness of the proposed approach is tested on a 6-bus test system and a 118-bus system.

1. Introduction

The climate change and global warming concerns coupled with high oil prices and government support are driving increasing integration of variable renewable energies into the power system [1]. Renewable energies have many advantages including reducing environmental pollution, bringing economic benefits, improving energy efficiency, and enhancing life quality [2].

The ongoing shale gas revolution is expected to slash nature gas prices and attract investment into gas-fired power generation. Take California’s generation mix as an example, natural gas provided approximately 53 percent of total installed capacity which is the largest source of energy to meet the load [3]. The increasing installation of gas-fired units significantly intensifies the interdependency between the electricity system and natural gas system. On one hand, gas-fired units can offer flexible dispatch and quick ramping capability [4]. On the other hand, the gas outages or gas shortage damage the reliability of power system as gas-fired units go off-line [5].

The joint scheduling of electricity and natural gas infrastructures has become an attractive option for improving the energy efficiency of the whole system [6]. In recent years, there have been many works dedicated to the joint optimal scheduling of the integrated energy system. Some studies are summarized as follows:

(1) Gas Flow Modeling: In order to ensure the computational tractability, the Weymouth equation is widely adopted to describe the relationship between gas flow and gas pressure under the steady-state assumption. Liu et al. propose a security-based methodology for the solution of short-term SCUC with considering the impact of the natural gas transmission system by using natural gas transmission subproblem to check the feasibility as well as the security [7]. In reference [8], authors develop a novel gas flow model such that line pack, i.e., the stored gas in the pipeline, is taken into account. That paper contributes with an extended incremental method that linearizes the nonlinear term in the Weymouth equation. However, such method introduces a large number of additional variables/constraints, posing a high computational burden to the solver. In paper [9], the Weymouth equation is convexified via a second order cone relaxation method, but the exactness of the relaxation cannot be guaranteed. In references [10,11], natural gas transmission is simplified by linear constraints, making the complex gas transmission problem tractable.

(2) Uncertainty Modeling: The increasing integration of variable renewable energy (VRE) into the power system brings uncertainties to the operation, which may result in load shedding, grid instability and potentially, cascading outages. In the scope of this paper, the VREs’ power output is the only source of uncertainties. How to manage the VREs’ uncertainties is a hot issue in the research field of joint scheduling. Alabdulwahab et al. build a scenario-based stochastic scheduling model for the electricity and natural gas systems [12]. A decomposition method is designed to check the feasibility of gas transmission problem in all the scenarios. In reference [13], the operation of a GB integrated gas and electricity network considering the wind power uncertainty is investigated using three operational methods: deterministic, two-stage stochastic programming, and multistage stochastic programming. Based on the interval optimization, Bai et al. propose a joint scheduling model with considering an optimal demand response strategy [14]. With resorting to the robust optimization method, He et al. propose a robust co-optimization scheduling model to study the coordinated optimal operation of the integrated energy systems [15]. That paper is tackled via an alternating direction method of multipliers (ADMM) by iteratively solving a power system subproblem. In this way, the privacy of data is protected.

In the existing literature focusing on joint scheduling, stochastic programming and robust optimization are two popular tools for addressing uncertainties. Stochastic programming assumes that the distribution governing the uncertainty data is known or can be estimated in advance. In contrast, robust optimization assumes that the probabilistic information is not available, and only the upper and lower bounds of uncertainty variables are determined [16]. Therefore, robust optimization may lead to an over-conservative decision. Recently, a novel optimization tool, named ’distributionally robust optimization’, is proposed to bridge the gap between stochastic programming and robust optimization. It constructs an ambiguity set containing the possible distributions of random variables. Then, it optimizes the expected cost under the worst-case distribution within the ambiguity set. In references [17,18], authors establish the moment-based ambiguity set where the moment information of possible distributions is consistent with that of the empirical distribution. A statistical distance-based metric– Wasserstein metric is employed in paper [19,20,21] to construct the ambiguity set.

This paper proposes a robust day-ahead scheduling model for electricity and natural gas system, which minimizes the total cost including fuel cost, spinning reserve cost and cost of operational risk while ensuring the feasibility for all scenarios within the uncertainty set. The proposed model is a two-stage robust optimization model. The base-case plan (base-case power output of units and preserved spinning reserve capacity) and the parameters of adjustable uncertainty set are determined in the first stage, while carrying out redispatch in the second stage after the realizations of uncertainties are revealed. Follow the study in [10], a simplification model is utilized to model the gas transmission with considering the line pack. The proposed model adopts an adjustable uncertainty set approach proposed in [22]. However, to hedge against the ambiguity of distribution, this paper proposes a risk-averse adjustable uncertainty set approach that computes the operational risk under the worst-case distribution by a distributionally robust method. Different from our previous work [23], this paper builds the ambiguity set via a Wasserstein–Moment metric [24,25]. Comparative studies are carried out to compute operational risk among the Wasserstein–Moment metric, Wasserstein metric, and Moment metric. In addition, this paper aims to coordinate different energy infrastructures but not the power system modeled in [23]. Our study has the following major contributions:

(1) We co-optimize energy and reserve for day-ahead scheduling in electricity and natural gas system. The gas transmission problem is simplified for achieving computational tractability.

(2) A risk-averse adjustable uncertainty set approach is proposed to mitigate the conservatism of robust optimization. To overcome distribution ambiguity, we construct the ambiguity set using a Wasserstein–Moment metric. Then, operational risk is the expected value under the worst-case distribution within such ambiguity set.

This paper is organized as follows: Section 2 presents the proposed methodology that includes the risk-averse adjustable uncertainty set approach and the robust day-ahead scheduling model for electricity and natural gas system; Numerical results are provided in Section 3.

2. Methodology

The flowchart of the proposed method is shown in Figure 1. First, we collect the historical data of forecasting errors. Then, based on the Wasserstein–Moment Metric, DRO models are constructed to compute the operational risk. By tuning the bounds of adjustable uncertainty set, we can form the operational risk curve for each VRE plant. Secondly, we build the robust day-ahead scheduling model for electricity and natural gas systems using the risk-averse adjustable uncertainty set approach. Such a model can be solved via the column-and-constraint generation method.

2.1. Risk-Averse Adjustable Uncertainty Set Approach

As the forecasting system has been an effective tool to predict the power output of VRE plant j, the forecasting error is the random variable which is denoted as ${\xi }_j$. Suppose that we have a historical dataset $\Xi$ for ${\xi }_j$ with sample size N, i.e., $\Xi =\{{\widehat{\xi }}_{j,1},...{\widehat{\xi }}_{j,N}\}$. Let ${\widehat{\mu }}_j$ and ${\widehat{\sigma }}_j$ represent the sample mean and the sample standard deviation of ${\xi }_j$, respectively. We can construct an adjustable uncertainty set ${\mathcal{V}}_j({\delta }_j)$:

$${\mathcal{V}}_j({\delta }_j)=[{\widehat{\mu }}_j-{\delta }_j{\widehat{\sigma }}_j,{\widehat{\mu }}_j+{\delta }_j{\widehat{\sigma }}_j]$$

(1)

where ${\delta }_j$ is a variable used to tune the scope of ${\mathcal{V}}_j({\delta }_j)$. According to the definition of robust optimization, for all the scenarios within the ${\mathcal{V}}_j({\delta }_j)$, the redispatch feasibility should be ensured. In other words, any scenario falling outside the ${\mathcal{V}}_j({\delta }_j)$ may cause the operational loss ([22,23]). In statistical decision theory, the risk function is defined as the expected value of a given loss function as a function of the decisions in the face of uncertainty. Following such definition, when fixing the decision ${\delta }_j$ as ${\widehat{\delta }}_j$, the operational risk is defined as an expected value of the operational loss under the empirical distribution ${\widehat{\mathbb{P}}}_N$, shown as follow:

$$R_j^{*}={\mathbb{E}}_{{\widehat{\mathbb{P}}}_N}\{\max (\rho ({\xi }_j-{\widehat{\mu }}_j-{\widehat{\delta }}_j{\widehat{\sigma }}_j),\tau ({\widehat{\mu }}_j-{\widehat{\delta }}_j{\widehat{\sigma }}_j-{\xi }_j),0)\}$$

(2)

where $\rho$ and $\tau$ are the wind power curtailment price and load shedding price, respectively. Nonetheless, it may result in a suboptimal solution, in particular in the case when limited historical data are available. To solve this problem, we construct an ambiguity set $\mathcal{B}$ to contain the possible distribution with calculating the operational risk under the worst-case distribution in such ambiguity set. The corresponding operational risk ${\widehat{R}}_j$ is:

$${\widehat{R}}_j=\underset{\mathbb{P}\in \mathcal{B}}{sup}\mathbb{E}\{\max (\rho ({\xi }_j-{\widehat{\mu }}_j-{\widehat{\delta }}_j{\widehat{\sigma }}_j),\tau ({\widehat{\mu }}_j-{\widehat{\delta }}_j{\widehat{\sigma }}_j-{\xi }_j),0)\}$$

(3)

This paper introduces a Wasserstein–Moment metric ([24,25]) to construct an ambiguity set ${\mathcal{B}}_{\epsilon }({\widehat{\mathbb{P}}}_N)$ where ${\widehat{\mathbb{P}}}_N$ represents the empirical distribution formed by $\Xi$.

$${\mathcal{B}}_{\epsilon }({\widehat{\mathbb{P}}}_N)=\{\mathbb{P}\in \Phi ;d_W(\mathbb{P},{\widehat{\mathbb{P}}}_N)⩽\epsilon ,{\mathbb{E}}_{\mathbb{P}}[({\tilde{\xi }}_j-{\widehat{\mu }}_j)({\tilde{\xi }}_j-{\widehat{\mu }}_j)⩽\gamma {{\widehat{\sigma }}_j}^2]\}$$

(4)

where $d_W()$ is the function to compute the Wasserstein distance between two distributions; $\Phi$ denotes the set containing all the possible distributions of $\tilde{\omega }$; $\gamma$ is a parameter to control conservatism, which can be obtained by the cross-validation method [17]. In this paper, the parameter $\epsilon$ is chosen as follow [20]:

$$\epsilon (N)=S\sqrt{2log(1/(1-\beta ))/N}$$

(5)

where S represents the diameter of the support of random variable; the parameter $\beta$ is set as 0.9, implying that the ambiguity set is able to cover the true distribution with the probability of $90\%$. According to Theorem 1, the operational risk w.r.t. ${\mathcal{B}}_{\epsilon }({\widehat{\mathbb{P}}}_N)$ can be obtained by solving a semi-definite linear program.

Theorem 1.

Suppose that random variable $\widetilde{{\xi }_j}$ has a convex support $\Omega =\{\mathit{D}\widetilde{{\xi }_j}⩽\mathit{g}\}$ where $\mathit{D}$ and $\mathit{g}$ are column vectors. The worst expectation of operational risk (3) is equivalent to:

$$\underset{\lambda ⩾0,\eta ⩾0,{\varrho }_n^{(1)}⩾0,{\varrho }_n^{(2)}⩾0,{\varrho }_n^{(3)}⩾0,{\mathit{\kappa }}_n^{(1)}⩾\mathbf{0},{\mathit{\kappa }}_n^{(2)}⩾\mathbf{0},{\mathit{\kappa }}_n^{(3)}⩾\mathbf{0}}{min}\frac{1}{N}\sum _{n=1}^N{\theta }_n+\lambda \epsilon +\eta \gamma {\widehat{\sigma }}^2$$

(6)

$$s.t.{\mathit{W}}_n\succcurlyeq 0,{\mathit{Y}}_n\succcurlyeq 0,{\mathit{Z}}_n\succcurlyeq 0,\forall n$$

(7)

$${\mathit{W}}_n=\left[\begin{array}{cc}\eta & (-\rho +{\varrho }_n^{(1)}+{\mathit{D}}^T{\mathit{\kappa }}_n^{(1)})/2-\eta \widehat{\mu }\\ (-\rho +{\varrho }_n^{(1)}+{\mathit{\kappa }}_n^{(1)T}\mathit{D})/2-\eta \widehat{\mu }& {\theta }_n+\rho ({\widehat{\mu }}_j+{\widehat{\delta }}_j{\widehat{\sigma }}_j)-{\varrho }_n^{(1)}{\widehat{\xi }}_{j,n}+\eta {\widehat{\mu }}^2-{\mathit{\kappa }}_n^{(1)T}\mathit{g}\end{array}\right],\forall n$$

(8)

$${\mathit{Y}}_n=\left[\begin{array}{cc}\eta & (\tau +{\varrho }_n^{(2)}+{\mathit{D}}^T{\mathit{\kappa }}_n^{(2)})/2-\eta \widehat{\mu }\\ (\tau +{\varrho }_n^{(2)}+{\mathit{\kappa }}_n^{(2)T}\mathit{D})/2-\eta \widehat{\mu }& {\theta }_n+\tau (-{\widehat{\mu }}_j+{\widehat{\delta }}_j{\widehat{\sigma }}_j)-{\varrho }_n^{(2)}{\widehat{\xi }}_{j,n}+\eta {\widehat{\mu }}^2-{\mathit{\kappa }}_n^{(2)T}\mathit{g}\end{array}\right],\forall n$$

(9)

$${\mathit{Z}}_n=\left[\begin{array}{cc}\eta & ({\varrho }_n^{(3)}+{\mathit{D}}^T{\mathit{\kappa }}_n^{(3)})/2-\eta \widehat{\mu }\\ ({\varrho }_n^{(3)}+{\mathit{\kappa }}_n^{(3)T}\mathit{D})/2-\eta \widehat{\mu }& {\theta }_n-{\varrho }_n^{(3)}{\widehat{\xi }}_{j,n}+\eta {\widehat{\mu }}^2-{\mathit{\kappa }}_n^{(3)T}\mathit{g}\end{array}\right],\forall n$$

(10)

$$-\lambda ⩽{\varrho }_n^{(1)}⩽\lambda ,-\lambda ⩽{\varrho }_n^{(2)}⩽\lambda ,,-\lambda ⩽{\varrho }_n^{(3)}⩽\lambda ,\forall n$$

(11)

The proof can be derived from a lemma provided in the appendix of our prior work [25]. If we parameterize the adjustable variable $\delta$ as a linear interpolation of a sequence of fixed values, ${\widehat{\delta }}_{j,1}$, ${\widehat{\delta }}_{j,2}$,...,${\widehat{\delta }}_{j,N_a}$. Accordingly, by using Theorem 1, we obtain a sequence of possible operational risk, ${\widehat{R}}_{j,1}$, ${\widehat{R}}_{j,2}$,...${\widehat{R}}_{j,N_a}$. Despite the fact that $R_j$ is a nonlinear function of ${\delta }_j$, it can be approximated by the piecewise linear function shown as follows:

$$R_j={\widehat{R}}_{j,1}+\sum _{m=1}^{n_a-1}({\widehat{R}}_{j,m+1}-{\widehat{R}}_{j,m}){\zeta }_{j,m}$$

(12)

$${\delta }_j={\widehat{\delta }}_{j,1}+\sum _{m=1}^{n_a-1}({\widehat{\delta }}_{j,m+1}-{\widehat{\delta }}_{j,m}){\zeta }_{j,m}$$

(13)

$$0⩽{\zeta }_{j,m}⩽1,\forall m$$

(14)

$${\zeta }_{j,m+1}⩽{\nu }_{j,m},{\nu }_{j,m}⩽{\zeta }_{j,m},m=1,2,...n_a-2$$

(15)

where ${\zeta }_m$ and ${\nu }_m$ are auxiliary continuous variable and binary variable w.r.t. the mth interval of piecewise linear curve. Constraint (15) ensures that there will be at most one index m where we have $0<{\zeta }_{j,m}<1$. Any given point ${\widehat{\delta }}_{j,m}⩽{\delta }_j⩽{\delta }_{j,m+1}$ can be uniquely represented by (12), as a linear combination of two consecutive breakpoints weighted by the associated variables ${\zeta }_{j,m}$. The nonlinear function is then approximated by the model as convex combination of the function values evaluated at the breakpoints weighted by the associated variables ${\zeta }_{j,m}$. Hereinafter, for simplicity, model (12)–(15) is written to the general formulation:

$$R_j=F_j({\delta }_j)$$

(16)

Through (16), we can plot an operational risk curve which describes the $R_j$ as the function of ${\delta }_j$. Obviously, increasing ${\delta }_j$ enlarges the scope of the adjustable uncertainty set, resulting in a lower operational risk. That is, $R_j$ is a non-increasing function of ${\delta }_j$.

2.2. Robust Day-Ahead Scheduling Model for Electricity and Natural Gas System

The proposed model aims to minimize the total cost including fuel cost, spinning reserve cost and cost of operational risk while satisfying the operating constraints. Since the operational risk is explicitly added in the objective function, the decision-maker achieves the optimal decision by making a tradeoff between economics and reliability. Denote ${\mathcal{G}}_g$ to represent the set of gas-fired units. Let ${\mathcal{G}}_{ng}$ be the set of non-gas thermal units. $\mathcal{W}$ represents the set of VRE plant. ${\mathcal{G}}_w$ denotes the set of gas wells. The objective function is written as below:

$$\min \{\sum _{t=1}^T(\sum _{i\in {\mathcal{G}}_{ng}}(c_{it}+{\overline{c}}_i{\overline{r}}_{it}+{\underline{c}}_i{\underline{r}}_{it})+\sum _{i\in {\mathcal{G}}_g}({\overline{c}}_i{\overline{r}}_{it}+{\underline{c}}_i{\underline{r}}_{it})+\sum _{k\in {\mathcal{G}}_w}{C}_k{g}_{kt}^w+\sum _{j\in \mathcal{W}}{R}_{jt})\}$$

(17)

where $c_{it}$ is the fuel cost of unit i at time t; $C_k$ is the gas price for gas well k; ${\underline{r}}_i$ and ${\overline{r}}_i$ are the downward spinning reserve and upward spinning reserve, respectively; ${\underline{c}}_i$ and ${\overline{c}}_i$ are the spinning reserve price; $R_{jt}$ is the operational risk w.r.t. VRE plant j at time t, which is calculated as discussed in Section 2.

2.2.1. First-Stage Constraints

In the first stage, the decision-maker optimizes the decision related to the base-case dispatch plan. The constraints are given as follows:

(1) Fuel cost of non-gas thermal unit i at time t.

The fuel cost curve is approximated by the piecewise linear function.

$$c_{it}⩾a_{is}^c{p}_{it}+b_{is}^c{u}_{it},i\in {\mathcal{G}}_{ng},i\in {\mathcal{G}}_g,\forall t,\forall s$$

(18)

where $a_{is}^c$ and $b_{is}^c$ are the slope and intercept of the s-th segment; $p_{it}$ is the base-case power output of unit i; $u_{it}$ is the binary variable representing the state of unit i.

(2) Power balance constraints.

$$\sum _{i\in {\mathcal{G}}_g,{\mathcal{G}}_{ng}}{p}_{it}+\sum _{j\in \mathcal{W}}{p}_{jt}^w=\sum _{b\in \mathcal{B}}{p}_{bt}^d,\forall t$$

(19)

where $p_{jt}^w$ is the forecasted power output of VRE plant j; $p_{bt}^d$ is the forecasted power of load at bus b; $\mathcal{B}$ is the set of buses.

(3) Power output limits:

$$p_{it}+{\overline{r}}_{it}⩽{\overline{p}}_i{u}_{it},{\underline{p}}_i{u}_{it}⩽p_{it}-{\overline{r}}_{it},i\in {\mathcal{G}}_{ng},i\in {\mathcal{G}}_g,\forall t$$

(20)

where ${\underline{p}}_i$ and ${\overline{p}}_i$ are the lower limit and upper limit of power output of unit i, respectively.

(4) Spinning reserve limits:

$$0⩽{\overline{r}}_{it}⩽{\overline{R}}_i,0⩽{\underline{r}}_{it}⩽{\underline{R}}_i,i\in {\mathcal{G}}_{ng},i\in {\mathcal{G}}_g,\forall t$$

(21)

where ${\underline{R}}_i$ and ${\overline{R}}_i$ are the lower limit and the upper limit of spinning reserve of unit i, respectively.

(5) Transmission line capacity contraints

$${\underline{f}}_l⩽\sum _{i\in {\mathcal{G}}_{ng},{\mathcal{G}}_g}{\pi }_{li}{p}_{it}+\sum _{j\in \mathcal{W}}{\pi }_{lj}{p}_{jt}^w+\sum _{b\in \mathcal{B}}{\pi }_{lb}{p}_{bt}^d⩽{\overline{f}}_l,\forall l,\forall t$$

(22)

where ${\pi }_{li}$, ${\pi }_{lj}$, and ${\pi }_{lb}$ are the power transfer distribution factors from unit i, VRE plant j, and load at bus b to line l, respectively. ${\underline{f}}_l$ and ${\overline{f}}_l$ are capacity limits of line l.

(6) Gas consumption of gas-fired units.

For gas-fired unit i, its gas consumption $g_{it}^c$ depends on power output $p_{it}$. In order to develop the piecewise linear approximation of such functional relationship, let us divide the axis $p_{it}$ over the range of function to $n_b$-1 segments using m = 1,2,...$n_b$ breakpoints ${\widehat{x}}_i^1$, ${\widehat{x}}_i^2$,... ${\widehat{x}}_i^{n_b}$. Each mth interval is associated with a continuous variable $p_{itm}$ and a binary variable ${\nu }_{itm}$.

$$g_{it}^c=\sum _{m=1}^{n_b-1}(a_{im}^g{p}_{itm}+b_{im}^g{\nu }_{itm})$$

(23)

$${\widehat{x}}_i^m{\nu }_{itm}⩽p_{itm}⩽{\widehat{x}}_i^{m+1}{\nu }_{itm}$$

(24)

$$\sum _{m=1}^{n_b-1}{p}_{itm}=p_{it}$$

(25)

$$\sum _{m=1}^{n_b-1}{\nu }_{itm}=u_{it}$$

(26)

where $a_{im}^g$ and $b_{im}^g$ denote the slope and intercept of the m-th segment. Constraint (26) ensures that only one binary variable ${\nu }_{itm}$ can take the value 1, and thus only one segment is active. Constraint (24) and (25) enforce that $p_{itm}$ is equal to $p_{it}$ if mth segment is active, otherwise $p_{itm}=0$.

(7) Gas nodal balance constraint.

$$\sum _{i\in {\phi }_d(b_g)}{g}_{it}^c+g_{b_{g}t}^d=\sum _{w\in {\phi }_w(b_g)}{g}_{kt}^w+\sum _{l_g\in \chi (b_g)}{G}_{b_g{l}_{g}t},\forall b_g,\forall t$$

(27)

where $g_{b_{g}t}^d$ is the gas load at node $b_g$; $\chi (b_g)$ represents the set of pipelines connected with gas node $b_g$; ${\phi }_d(b_g)$ denotes the set of gas-fired units plugged at node $b_g$; ${\phi }_w(b_g)$ is the set of gas well at node $b_g$; $G_{b_g{l}_{g}t}$ is the gas flow injecting into node $b_g$ from pipeline $l_g$. We note that $G_{b_g{l}_{g}t}$ is positive when gas flow runs from pipeline $l_g$ to $b_g$, otherwise it is a negative value.

(8) Line pack temporal constraints.

Line pack refers to the quantity of natural gas stored in the pipeline. Line pack can be used to smooth the real-time gas demand variation and maintain gas pressure stability.

$$L{P}_{l_{g}t}-L{P}_{l_{g}t-1}+G_{O(l_g)l_{g}t}+G_{D(l_g)l_{g}t}=0,\forall l_g,\forall t$$

(28)

where $O(l_g)$ denotes the origin node of pipeline $l_g$; $D(l_g)$ denotes the destination node of pipeline $l_g$; $L{P}_{l_{g}t}$ is the linepack of pipeline $l_g$ at time t.

(9) Line pack/gas supply/outgoing gas flow limits constraints.

$${\underline{LP}}_{l_g}⩽L{P}_{l_{g}t}⩽{\overline{LP}}_{l_g},\forall l_g,\forall t$$

(29)

$${\underline{g}}_k^w⩽g_{kt}^w⩽{\overline{g}}_k^w,\forall k,\forall t$$

(30)

$${\underline{G}}_{O(l_g)l_g}⩽G_{O(l_g)l_{g}t}⩽{\overline{G}}_{O(l_g)l_g},{\underline{G}}_{D(l_g)l_g}⩽G_{D(l_g)l_{g}t}⩽{\overline{G}}_{D(l_g)l_g},\forall l_g,\forall t$$

(31)

where the upper bound and lower bound of a variable is written with an overline and an underline.

(10) Initial line pack constraints.

$$L{P}_{l_{g}0}=L{P}_{l_{g}T}={\widehat{LP}}_{l_g},\forall l_g$$

(32)

Constraints (32) enforce the initial line pack is equal to the final line pack to facilitate the daily operation.

(11) Constraints associated with operational risk, i.e., constraints (16).

2.2.2. Second-Stage Constraints

After VREs’ power outputs are realized, system operators carry out corrective actions e.g., redispatching and load shedding. According to the definition of robust optimization, any scenario within the adjustable uncertainty set will not violate the operating constraints.

(1) Real-time released spinning reserve.

The real-time adjustment of generators should be within its spinning reserve preserved in the first stage.

$$0⩽p_{it}^{+}⩽{\overline{r}}_{it},0⩽p_{it}^{-}⩽{\underline{r}}_{it},i\in {\mathcal{G}}_{ng},i\in {\mathcal{G}}_g,\forall t$$

(33)

where $p_{it}^{+}$ and $p_{it}^{-}$ are the upward corrective power and downward corrective power, respectively.

(2) Real-time power balance.

$$\sum _{i\in {\mathcal{G}}_g,{\mathcal{G}}_{ng}}(p_{it}+p_{it}^{+}-p_{it}^{-})+\sum _{j\in \mathcal{W}}(p_{jt}^w+{\xi }_{jt})=\sum _{b\in \mathcal{B}}{p}_{bt}^d,\forall {\xi }_{jt}\in {\mathcal{V}}_{jt}({\delta }_{jt}),\forall t$$

(34)

where ${\mathcal{V}}_{jt}$(${\delta }_{jt}$) represents the adjustable uncertainty set w.r.t. VRE plant j at time t.

(3) Real-time transmission line capacity.

$${\underline{f}}_l⩽\sum _{i\in {\mathcal{G}}_{ng},{\mathcal{G}}_g}{\pi }_{li}(p_{it}+p_{it}^{+}-p_{it}^{-})+\sum _{j\in \mathcal{W}}{\pi }_{lj}(p_{jt}^w+{\xi }_{jt})+\sum _{b\in \mathcal{B}}{\pi }_{lb}{p}_{bt}^d⩽{\overline{f}}_l,\forall {\xi }_{jt}\in {\mathcal{V}}_{jt}({\delta }_{jt}),\forall l,\forall t$$

(35)

(4) Real-time gas consumption of gas-fired unit i

For gas-fired unit i, the real-time gas consumption is assumed to be a linear function of the real-time power output at the active segment (${\nu }_{itm}=1$).

$$a_{im}^g(p_{it}+p_{it}^{+}-p_{it}^{-})+b_{im}^g{u}_{it}-(1-{\nu }_{itm})M⩽{\tilde{g}}_{kt}^w⩽a_{im}^g(p_{it}+p_{it}^{+}-p_{it}^{-})+b_{im}^g{u}_{it}+(1-{\nu }_{itm})M,i\in {\mathcal{G}}_g,\forall t,\forall m$$

(36)

where M is a sufficient large constant; ${\tilde{g}}_{kt}^w$ is the real-time gas consumption. We mention that decision $\{p_{it},u_{it},{\nu }_{itm}\}$ are obtained in the first stage. It can be seen from constraint (36) that for the active segment (${\nu }_{itm}=1$), we have

$${\tilde{g}}_{kt}^w=a_{im}^g(p_{it}+p_{it}^{+}-p_{it}^{-})+b_{im}^g{u}_{it}$$

(37)

while constraint (36) is relaxed for other inactive segments.

In the second stage, there are constraints associated with the real-time operating of the gas system, which is similar to the constraints shown in (27)–(32). For simplicity, we do not present them here.

The proposed joint scheduling model is written into the compact form:

$$\underset{\mathit{x}\in \mathit{X},\mathit{y},\mathit{R},\mathit{\delta }}{\min }\{{\mathit{c}}^T\mathit{x}+{\mathbf{1}}^T\mathit{R}\}$$

(38)

$$\mathrm{s}.\mathrm{t}.R_{jt}=F_{jt}({\delta }_{jt}),\forall j,\forall t$$

(39)

$$\mathit{B}\mathit{y}⩽\mathit{b}-\mathit{A}\mathit{x}-\mathit{D}\mathit{\xi },\forall \mathit{\xi }\in \mathcal{V}$$

(40)

$$\mathcal{V}=\{\widehat{\mathit{\mu }}-\mathit{\delta }\circ \widehat{\mathit{\sigma }}⩽\mathit{\xi }⩽\widehat{\mathit{\mu }}+\mathit{\delta }\circ \widehat{\mathit{\sigma }}\}$$

(41)

where $\mathit{x}$ is the vector of variables in the first stage; $\mathit{y}$ is the vector of variables in the second stage; $\mathit{X}$ is the feasible region formed by the constraints in the first stage; $\mathit{R}$ is the vector of the operational risks $\{R_{jt}\}$; $\mathit{\xi }$ is the vector of random variables; $\mathit{\delta }$ is the vector of the parameters $\{{\delta }_{jt}\}$; $\widehat{\mathit{u}}$ and $\widehat{\mathit{\sigma }}$ are vectors of sample means $\{{\widehat{\mu }}_j\}$ and sample standard deviations $\{{\widehat{\sigma }}_j\}$; $\mathcal{V}$ is the union of uncertainty sets $\{{\mathcal{V}}_{jt}\}$; ’∘’ denotes the Hadamard product. Matrices $\mathit{A}$, $\mathit{B}$ and $\mathit{D}$, and $\mathit{b}$ represent the parameters in the second-stage constraints. Constraints (40) correspond to the second-stage constraints.

2.3. Solution Methodology

The proposed joint scheduling model can be solved by column and constraints generation algorithm (C&CG) which has been proved to be an efficient solution methodology for two-stage robust optimization model [26]. Next, we briefly discuss the solution methodology. Interested readers can refer to our previous work [23] for detailed information.

Given the fixed first-stage decision $\{\widehat{\mathit{x}},\widehat{\mathit{R}},\widehat{\mathit{\delta }}\}$, the C&CG subproblem aims to identify the worse-case scenario by solving the following model:

$$\underset{\mathit{\xi }}{\max }\underset{\mathit{y},\mathit{s}}{\min }\{{\mathbf{1}}^T\mathit{s}\}$$

(42)

$$\mathrm{s}.\mathrm{t}.\mathit{B}\mathit{y}⩽\mathit{b}-\mathit{A}\widehat{\mathit{x}}-\mathit{D}\mathit{\xi }+\mathit{s}$$

(43)

$$\widehat{\mathit{\mu }}-\widehat{\mathit{\delta }}\circ \widehat{\mathit{\sigma }}⩽\mathit{\xi }⩽\widehat{\mathit{\mu }}+\widehat{\mathit{\delta }}\circ \widehat{\mathit{\sigma }}$$

(44)

$$\mathit{s}⩾\mathbf{0}$$

(45)

where $\mathit{s}$ is a vector of relaxation variables. With introducing the uncertainty budget $\Gamma$ and $\mathrm{H}$ to mitigate the conservatism of robust optimization, the above model is reformulated to the following model by dual theory:

$$\underset{\mathit{\kappa },{\mathit{z}}^{+},{\mathit{z}}^{-}}{\max }-{\mathit{\kappa }}^T[\mathit{b}-\mathit{A}\widehat{\mathit{x}}-\mathit{D}(\widehat{\mathit{\mu }}+\widehat{\mathit{\delta }}\circ \widehat{\mathit{\sigma }}\circ {\mathit{z}}^{+}-\widehat{\mathit{\delta }}\circ \widehat{\mathit{\sigma }}\circ {\mathit{z}}^{-})]$$

(46)

$$\mathbf{0}⩽\mathit{\kappa }⩽\mathbf{1}$$

(47)

$${\mathit{B}}^T\mathit{\kappa }=\mathbf{0}$$

(48)

$${\mathit{z}}^{+}+{\mathit{z}}^{-}⩽\mathbf{1}$$

(49)

$$\sum _{j\in \mathcal{W}}(z_{jt}^{+}+z_{jt}^{-})⩽\Gamma ,\forall t$$

(50)

$$\sum _{t=1}^T(z_{jt}^{+}+z_{jt}^{-})⩽\mathrm{H},\forall \mathrm{j}$$

(51)

where binary column vectors ${\mathit{z}}^{+}=\{z_{jt}^{+}\}$ and ${\mathit{z}}^{-}=\{z_{jt}^{-}\}$ are indicators of the bounds of $\mathit{\xi }$; $\mathit{\kappa }$ is the vector of dual variables; $\Gamma$ is the spatial budget; $\mathrm{H}$ is the temporal budget. When the subproblem is solved in the kth iteration, denoting the optimal solution of ${\mathit{z}}^{+}$ and ${\mathit{z}}^{+}$ as ${\widehat{\mathit{z}}}^{+}$ and ${\widehat{\mathit{z}}}^{-}$, respectively, a set of extra variables ${\mathit{y}}^k$, and additional constraints are generated as shown in (21) and added to the master problem.

$$\mathit{B}{\mathit{y}}^k⩽\mathit{b}-\mathit{A}\mathit{x}-\mathit{D}(\widehat{\mathit{\mu }}+\mathit{\delta }\circ \widehat{\mathit{\sigma }}\circ {\widehat{\mathit{z}}}^{+}-\mathit{\delta }\circ \widehat{\mathit{\sigma }}\circ {\widehat{\mathit{z}}}^{-})$$

(52)

Suppose that K iterations have been proceeded, the C&CG master problem is given as follows:

$$\underset{\mathit{x}\in \mathit{X},\mathit{y},\mathit{R},\mathit{\delta }}{\min }\{{\mathit{c}}^T\mathit{x}+{\mathbf{1}}^T\mathit{R}\}$$

(53)

$$\mathrm{s}.\mathrm{t}.R_{jt}=F_{jt}({\delta }_{jt}),\forall j,\forall t$$

(54)

$$\mathit{B}{\mathit{y}}^k⩽\mathit{b}-\mathit{A}\mathit{x}-\mathit{D}(\widehat{\mathit{\mu }}+\mathit{\delta }\circ \widehat{\mathit{\sigma }}\circ {\widehat{\mathit{z}}}^{+}-\mathit{\delta }\circ \widehat{\mathit{\sigma }}\circ {\widehat{\mathit{z}}}^{-}),k=1,2,\dots K$$

(55)

The above iterative process will continue until the optimal value of subproblem is 0, which means that all the scenarios within the adjustable uncertainty set are feasible.

3. Numerical Results

Numerical tests on the integrated energy system including a 6-bus power system and a 6-node natural gas system is carried out as an illustrative example. A larger system containing an IEEE 118-bus power system and a 10-node natural gas system are used to verify the computational efficiency of the proposed method. The semi-definite linear programs are solved with MOSEK, and mixed integer linear programming models are solved with GUROBI. All the tests are performed on a PC, which has 8GB memory and i7-4790 CPU.

3.1. Six-Bus Power System with a 6-Node Natural Gas System

We use the test system in [15] with a trivial modification. The topologies are shown in Figure 2. The 6-bus power system includes non-gas thermal units G2–G4, gas-fired units G1 and G5, a wind farm connecting to bus 4, seven transmission lines, and three loads. Fuel price of G2–G4 is 5 $/MBtu. The 6-node gas system comprises two gas wells, five pipelines, and four gas loads. The wind power curtailment prices and load shedding prices are assumed to vary with time periods, which are plotted in Figure 3. The total power loads and total gas loads are also illustrated in Figure 3. Data for this integrated energy system are provided in [27]. Budget $\Gamma$ and $\mathrm{H}$ are set as 1 and 24, respectively.

A wind farm is connected to the power system at bus 4 with a capacity of 240 MW. The forecasted power output is plotted in Figure 4. Following the setting in [22], the forecasting error of wind generation is assumed to follow a Gaussian distribution where the mean value is 0 and the standard deviation is:

$${\sigma }_{jt}=0.1·p_{jt}^w·(1+e^{(T-t)})$$

(56)

where $p_{jt}^w$ is the forecasted power output of VRE plant j. We make the assumption on the ’true’ distribution to sample data for constructing ambiguity sets. In this way, instead of the ’true’ distribution, the sampled data are available for the decision-making and thus the proposed approach essentially does not rely on the specified distributions.

The simulation procedure is given as follows: (1) Sample N = 5000 data from the assumed ’true’ Gaussian distribution for constructing ambiguity sets; (2) The typical points on the operational risk curve are determined by solving (6)–(11); (3) Build the optimization model and solve it; (4) By keeping the base-case plan fixed, calculate the expected total cost under the ’true’ distribution to test the out-of-sample performance of the decision.

In order to investigate the conservatism of Wasserstein–Moment-metric (WM-metric), we introduce moment-metric (M-metric) [17] and Wasserstein-metric (W-metric) [23]. All tests are based on the same sample set for fair comparison. Monte Carlo Simulation (MCS) with ${10}^6$ samples is employed to provide the benchmark for comparison. Then, based on several metrics, the operational risk curves at hour 1 and hour 24 are plotted in the Figure 5.

From Figure 5, it can be seen that WM-metric is the most aggressive one among the studied metrics. The reason is that both of statistical distance and moment information are taken into account, reducing the ambiguity of random variables. In general, the more probabilistic information is used to construct the ambiguity set, the worst-case distribution will be less conservative.

Based on the above operational risk curves, the following four cases are used to investigate the performance of WM-metric in making dispatch decision:

Case 1: Conventional two-stage robust optimization model where the parameters of predefined uncertainty sets are selected as in [28] (95% confidence intervals are adopted);

Case 2: Proposed scheduling model with WM-metric;

Case 3: Proposed scheduling model with W-metric;

Case 4: Proposed scheduling model with M-metric;

The unit commitment decision is illustrated in Figure 6. It is observed that Cases 2–4 result in the same decision while Case 1 commits more unit in hour 1, 2, 3 and 23. It indicates that the conventional two-stage robust model with unadjustable uncertainty set leads to a more conservative decision compared with the proposed model. Cases 2–4 yield different spinning reserve solutions. Figure 7 shows the total upward spinning reserve for each hour. Case 2 preserves significantly less spinning reserve than other cases due to its less conservative estimation of operational risk. The downward spinning reserve in these cases is not shown here because we come to similar conclusions.

The results of costs are shown in

Table 1. Results of costs (k$).

Case	Production Cost	Estimated Total Cost		Simulated Total Cost
		Operational Risk Cost	Total Cost	Expected Loss Cost	Total Cost
Case 1	589.7	\	589.7	0.2	589.9
Case 2	581.0	6.6	587.6	2.5	583.5
Case 3	589.3	1.5	590.8	0.3	589.6
Case 4	588.4	2.3	590.7	0.5	588.9

. We compare two types of total costs. The first type is the estimated total cost which is the sum of the production cost (fuel cost and spinning reserve cost) and operational risk cost. Actually, the first type cost is the objective value of (17). By fixing the base-case dispatch plan, we solve the re-dispatch problem with ${10}^5$ samples generated from the assumed distribution, and record the expected risk cost including load shedding cost and wind power curtailment cost. The second type total cost, named simulated total cost, is equal to the sum of the production cost and expected risk cost.

From Table 1, it can be seen that Case 2 achieves the lowest production cost as well as the total cost. The reason is that the spinning reserve scheduled in Case 2 is less than that in other cases as shown in Figure 7. Compared with other cases, the production cost is reduced in Case 2 at the price of increasing the risk cost. However, Case 2 compromises reliability due to the lowest total cost. It is suggested that using WM-metric is able to mitigate the conservatism of the risk-averse adjustable uncertainty set approach so that Case 2 makes the best tradeoff between economics and reliability among all the cases.

Next, we investigate the effect of gas shortage on the unit commitment. We assume that some contingencies restrict the capacity of the pipeline directed from node 1 to node 2 so that gas supply cannot satisfy the demand of the gas-fired unit at node 1. As a result, this gas-fired unit cannot be committed but replaced by costly non-gas thermal units. Under such assumption, the unit commitment decision is presented in Figure 8 . Cases 2–4 lead to the same feasible decision while Case 1 is not feasible, implying that the adjustable uncertainty set approach can be a promising alternative method in particular cases when the conventional robust methods cannot provide a feasible solution.

Table 2. Results of costs for gas shortage example (k$).

Case	Production Cost	Estimated Total Cost		Simulated Total Cost
		Operational Risk Cost	Total Cost	Expected Loss Cost	Total Cost
Case 1 (infeasible)	\	\	\	\	\
Case 2	891.0	13.1	904.1	8.1	899.1
Case 3	889.9	20.4	910.3	9.2	899.1
Case 4	892.9	17.9	910.8	6.7	899.6

summarizes the costs for each case. Case 2 also achieves the lowest total cost as in Table 1. From the comparison between Table 1 and Table 2, we can find that the total costs are greatly increased because of the shutdown of cost-efficient unit 1. It indicates that the gas shortage severely reduces the operation efficiency of the integrated enery system.

3.2. IEEE 118-Bus Power System with a 10-Node Natural Gas System

We use this system to test the computational performance of the proposed method. The modified IEEE-118 bus power system consists of 46 non-gas thermal units, eight gas-fired units, six wind farms, 186 branches, and 91 loads. We simplify the gas system in [15] by reducing storage facilities and compressors. The modified gas system includes ten gas nodes, three gas wells, ten pipelines, and 12 gas loads. Six wind farms are connected to the system at buses 12, 17, 49, 59, 80, and 92 with a capacity of 200 MW for each bus. Uncertainty budget $\Gamma$ is set as 3 to reduce conservatism, which means that there are at most three wind farms simultaneously experience fluctuations. The four cases studied in the previous subsection are carried out, the results are listed in

Table 3. Results of costs in a larger system (k$).

Case	Production Cost	Estimated Total Cost		Simulated Total Cost
		Operational Risk Cost	Total Cost	Expected Loss Cost	Total Cost
Case 1	10,262.4	\	10,262.4	1.1	10,263.5
Case 2	10,248.1	13.1	10,261.2	5.2	10,253.3
Case 3	10,254.5	9.9	10,264.4	3.1	10,257.6
Case 4	10,257.2	7.4	10,264.6	2.8	10,260.0

. Like the results in the 6-bus test system, the scheduling model with WM-metric is less conservative than other models and achieves a better balance between reliability and economics.

We obtain the operational risk curves in the preparation stage. These problems are highly parallelizable among different time periods and wind farms, and therefore, we can speed up the computation via solving these problems in a parallel manner. Owing to the parallel manner, time costs are listed in

Table 4. Computational efficiency.

Case	Time of Preparation (s)	Time of C & CG (s)
Case 1	\	91.2
Case 2	921	125.1
Case 3	54	114.4
Case 4	10.9	102.9

. For all the cases, the time costs spent on solving scheduling models are comparable. Nonetheless, Case 2 spends more time forming operational risk curves than other cases because it is inefficient to tackle semi-definite programming constraints.

4. Conclusions

This paper proposes a robust day-ahead scheduling model for electricity and natural gas systems, which minimizes the total cost including fuel cost, spinning reserve cost and cost of operational risk while ensuring the feasibility for all scenarios within the uncertainty set. From the numerical results, we conclude that:

(1)

The proposed scheduling model can coordinate different energy infrastructures to improve the energy efficiency of the electricity and natural gas system. The gas shortage has a significant adverse impact on the economic efficiency.

(2)

The risk-averse adjustable uncertainty set approach can be a promising alternative method in particular cases when the conventional robust methods cannot provide a feasible solution.

(3)

Compared with the other two metrics, WM-metric helps to reduce the conservatism of risk-averse adjustable uncertainty set approach.

(4)

The proposed scheduling model can be solved efficiently via the C&CG algorithm.

The proposed method can be used by the system operator to execute day-ahead joint scheduling (or market clearing) of the electricity and natural gas systems. The performance of the proposed method highly depends on the data quality because available data are used to construct ambiguity sets for computing operational risk. In this regard, the data cleaning and the outliers detection method can be utilized to improve data quality. The proposed method does not specify the type of VRE plant, but deals with the random forecasting error. Therefore, it also can be applied to the system that includes other renewable units such as solar PV and biomass units, etc.