Co-optimisation and Settlement of Power-Gas Coupled System in Day-ahead Market under Multiple Uncertainties

The interdependency of power systems and natural gas systems is being reinforced by the emerging power-to-gas facilities (PtGs), and the existing gas-fired generators. To jointly improve the efficiency and security under diverse uncertainties from renewable energy resources and load demands, it is essential to co-optimise these two energy systems for day-ahead market clearance. In this paper, a data-driven integrated electricity-gas system stochastic co-optimisation model is proposed. The model is accurately approximated by sequential mixed integer second-order cone programming, which can then be solved in parallel and decentralised manners by leveraging generalised Benders decomposition. Since the price formation and settlement issues have rarely been investigated for integrated electricity-gas systems in an uncertainty setting, a novel concept of expected locational marginal value is proposed to credit the flexibility of PtGs that helps hedging uncertainties. By comparing with a deterministic model and a distributionally robust model, the advantage of the proposed stochastic model and the efficiency of the proposed solution method are validated. Detailed results of pricing and settlement for PtGs are presented, showing that the expected locational marginal value can fairly credit the contribution of PtGs and reflect the system deficiency of capturing uncertainties.


Motivation
Power-to-gas (PtG) is quite effective in storing large quantity of excess renewable electricity compared with conventional powerto-power energy storage technologies [1]. Given the high energy density of methane and the great potential of natural gas network as storages [2], PtG has been considered a promising technique in sustainable energy systems [1,3]. Besides, natural gas-fired units (GfUs), despite being traditional facilities, contribute an increasingly large share of the electricity generation [4,5]. The development of PtGs and the growth of GfUs tightly couple the electric power system with the natural gas system [6]. The electric power system and the natural gas system are conventionally operated as individual systems without sufficient coordinations, as they belong to different energy sectors. However, the intensified coupling has resulted in an integrated electricity-gas system (IEGS), for which coordinated operation become inevitable. Moreover, the liberalization of both the electricity market and the natural gas market [5,7,8], together with the interactive safety and reliability requirements of IEGS [6,9,10], are appealing for a security-constrained co-optimisation regime and corresponding settlement methods.
The challenges of co-optimizing IEGS in day-ahead markets include: i) the uncertainties from both renewable generations and electricity/gas demands, i) the non-convexity of the natural gas flow model, and iii) the requirement of decentralised decision making. Therefore, it is necessary to develop a model that simultaneously addresses the above-mentioned issues with desired accuracy and reliability. Also, efficient solution algorithm should be developed.
Another practical challenge is the pricing issue or the settlement of these two energy sectors. Settlement of IEGS is a rather new topic, especially when the uncertainties of renewable generations and load demands are accounted for. Under an uncertainty environment, the traditional price formation mechanism in day-ahead markets must be systematically reevaluated and improved, because the original pricing regime may not be equitable and incentive enough for market participants who provide flexibilities and reserves.

Problem Modeling and Solution Algorithm:
The stochastic day-ahead scheduling problem of IEGS is investigated by [11], in which the natural gas flow problem is solved independently by Newton-Raphson substitution method to generate natural gas usage cuts. However, PtGs and the line-pack effect are ignored. Ref. [12] proposes an interval-optimisation-based model for IEGS to handle the wind power uncertainty, which is then solved directly by a mixed integer nonlinear programming (MINLP) solver. A robust unit commitment (UC) model is developed in [13] to deal with the uncertainty of transmission line outage. Again, the demand uncertainty in the gas system is not considered, and both the line-pack effect and the gas compressor model are omitted in order for problem tractability. Ref. [14] deals with the optimal gas-power flow problem without considering the on/off statues of generators. Only the wind power uncertainty on the power system side is considered therein, which is addressed by distributionally robust optimisation. Ref. [15] proposes a hybrid scenario and information gap based optimisation model for the day-ahead co-optimisation of IEGS under multiple uncertainties, and the MINLP is solved with a commercial solver. Ref. [16] proposes an uncertainty management framework for IEGS, which leveraging chance-constrained optimisation and robust optimisation. The transient gas pipeline flows are accurately modeled in [16].
To address the non-convexity of the problem and enable decentralised solutions, apart from the linearisation via Taylor series expansion [17] and the second-order cone reformulation used in [13,14], Ref. [18] proposes using mixed integer second-order cone constraints to enhance the approximation of the non-convex gas flow equation. More recently, Ref. [19] proposes an outer approximation with equality relaxation method to cope with the non-convexity issue. In [20], the shadow price is utilised to coordinately optimize IEGS in day-ahead markets. In the robust IEGS model of [21], the non-convex natural gas problem is reformulated as a mixed integer linear programming (MILP), and the non-convex sub-problem of the robust optimisation model is solved distributedly via the alternating direction method of multipliers (ADMM) with heuristics. In a subsequent work [22], the authors introduce price-sensitive demandresponses, and the uncertainty is handled by distributionally robust optimisation based on the linearised natural gas model.

Pricing and Settlement:
Regarding the pricing and settlement issues, the authors in [23] propose a method for pricing the gas capacity reserved to GfUs. However, the non-convex gas transmission constraints are approximated by some linear cuts, and constraints in stochastic scenarios are discarded. The strategic offering and equilibrium problem of coordinated IEGS markets is investigated in [24], whereas the line-pack effect and gas nodal pressure variables are omitted for problem tractability. A scenario-based model is proposed in [25] to determine the optimal operation strategy of GfUs and PtGs in energy and regulation markets. Further, the Shapley value is employed to allocate the payoff among these facilities.
The concept of cost of uncertainty is developed in [26] to characterize the impact of uncertainty on the dispatch cost, but the value of flexible resources is not evaluated. In [27], the authors make use of the derivative of a robust UC model to construct the uncertainty marginal price, which quantifies the value of reserve and the cost of uncertainty in the day-ahead market. A recent work in [28] deals with the problem of pricing transmission overload and generation violation caused by random renewable generations and demands. Therein, a distributionally robust chance-constrained optimal power flow model is developed, which renders uncertaintycontained locational marginal prices that determine how the revenue should be distributed to conventional generators. More recently, Ref. [29] proposes a chance-constrained stochastic pricing method for linear electricity markets, in which the price is formed by a scenario-independent mathematical programming reduced from the chance-constrained model.

Contribution and Paper Organization
In this paper, a day-ahead co-optimisation problem of IEGS is investigated, considering the uncertainties of both renewable generations and electricity/gas demands. Moreover, the price formation and settlement issue is studied with a focus on PtGs, and the economic efficiency of PtGs is also analysed. The proposed co-optimisation method and settlement regime are validated by thorough numerical results and comparisons with a deterministic model and a distributionally robust model.
The detailed technical contributions of this paper include: 1)A stochastic day-ahead market model is developed for the integrated electricity-gas system, which precisely accounts for the natural gas flow constraints, line-pack effect, PtGs, as well as correlated uncertainties. The stochastic model makes use of data-driven scenarios so that the natures of multiple uncertainties could be better retained.
2)The stochastic model is approximated by sequential mixed integer second-order cone programming (MISOCP), which is shown to be highly precise. Based on generalised Benders decomposition, the convex sub-problems are further decoupled and solved by the electric power system operator and the natural gas system operator decentrally. The stochastic model and the ensemble solution method are shown to have advantages over state-of-the-arts in terms of dealing with the uncertainty, the non-convexity, and the decentralised decision making issues.

3)A novel concept of expected locational marginal value (E-LMV) is
proposed for price formation in the electricity-gas market, which has advantages in crediting PtGs equally and ensuring cost recovery of such flexibility providers in a power-gas coupled market with production and demand uncertainties. Moreover, the revenue adequacy condition of the day-ahead natural gas market is analysed for the first time.
The remainder of this paper is organised as follows: Section 2 and Section 3 introduce the natural gas system model, and the electric power system model, respectively. Section 4 establishes the stochastic model for the power-gas coupled system, and introduces the novel pricing method. Section 5 presents the solution algorithms. Numerical experiments and detailed results are reported in Section 6. Section 7 concludes with discussions.

Natural Gas System Model
This section presents a dynamic/multi-period optimal flow model of the natural gas system. Typical components are modeled including gas compressors, gas storages, GfUs and PtGs. The gas traveling velocity and compressibility are accounted for [9], as gas travels much slower than electricity and it can be stored in pipelines. Further, we assume by convention that the state variables of the natural gas system are stable within each 1-hour scheduling time slot [30].

Notation for Natural Gas system
N , T Set of nodes in the natural gas system and set of scheduling time periods. G Src n , G Str n Sets of natural gas suppliers and gas storages, connected at node n. G Cmp n , G Pipe n Sets of active pipelines (with gas compressors) and passive pipelines (without gas compressors) connected with node n; node n is the outlet of pipeline (m, n), or the inlet of pipeline (n, m).

sgn(·)
Sign function that returns -1 for negative input, 0 for zero, and 1 for positive input.  Non-generation-related natural gas demand at node n at time t [Mscm/h].

Model Formulation
The GS model is formulated as: f Str s,τ ≤ Ss ∀s ∈ G Str n , t ∈ T (1d) e (m,n),t = K lp (m,n) (π m,t + π n,t ) /2 ∀(m, n) ∈ G Pipe n , t ∈ T (1k) e (m,n),t = ∆ · f Pipe (m,n),t + f Pipe The objective function accounts for the gas volume from suppliers and the net gas consumption of storages. Constraints (1a) and (1b) define flow limits and daily quantity limits of gas sources. Constraints (1c) and (1d) define flow limits and capacity limits of gas storages. Constraint (1e) restricts the gas pressure of each node to be within a safety range. For active pipelines, compression ratios are limited by constraint (1f), while gas consumptions and line packs are defined respectively by constraints (1g) and (1h). For passive pipelines, the general flow equation (1i) expresses the relationship between the pressure gradient and the gas flow, which can be evaluated via (1j); Equality (1k) indicates that the line pack is proportional to the average pressure, and the line pack should also complies with the mass conservation (1l). Constraint (1m) imposes a requirement on line-pack level in the last scheduling period. Constraint (1n) enforces gas balance at each node. Constraint (1g) adopts a simplified gas consumption function for the compressor [18,30] instead of the original one, which is highly nonlinear on the gas flow through and the compression ratio [30]. Constraints (1i) and (1k) can also be applied to active pipelines after such pipelines are separated into two segments from the location of compressors, but this is necessary only when the length of a pipeline is considerable. Two key parameters of the natural gas flow model, i.e., K gf (m,n) and K lp (m,n) , are calculated according to the equations derived in the appendix of [9]. To obtain K lp (m,n) , the friction factor of pipeline is yielded from the Nikuradse equation first, which is detailed in [31]. Parameters used to calculate K gf (m,n) and K lp (m,n) are available online [32].

Electric Power System Model
This section presents a basic security-constrained unit commitment (SCUC) model for the electric power system. The reserve requirements are omitted herein since stochastic programming is used in this paper. Nevertheless, constraints for the reserve are retained in a deterministic model, which is adopted as benchmark in case studies.

Notation for Electric Power System
G, G CfU Sets of all units and coal-fired units (CfUs). L Set of transmission lines.

E, E ref
Sets of buses and reference bus.
Element on the m-row and the n-th column of the nodal susceptance matrix [kV 2 S]. κ G g,n , κ PtG v,n 0-1 coefficient indicating whether unit g or PtG v is connected at bus n. P g , P g Minimum and maximum production levels of unit g Binary variables indicating whether the unit is on. u g,t , v g,t Binary variables indicating whether the unit is started up and shut down.

Model Formulation
The SCUC model is formulated as: The objective function accounts for the start-up and shut-down costs of CfUs, and the generation cost of CfUs. Constraints In the electric power system model, the start-up and shut-down costs, as well as the on/off variables of GfUs are omitted. This is due to two facts. First, GfUs are quick-start units that can change their intra-day on/off statues, so it is inappropriate to fix their statues day-ahead. Second, the on/off statues can be ignored in the optimisation model without affecting the engineering behavior of GfUs (because the start-up time and minimum production level of GfUs are quite short/low), while incorporating binary variables into the convex dispatch problem will complicate the stochastic counterpart of this problem a lot (e.g., the dispatch problem becomes a mixed integer programming, to which many decomposition algorithms are no longer applicable). In practice, we can simply add a constant term to the objective function to account for the daily average start-up and shut-down costs of GfUs though.
It is worth mentioning that in the implementation of the models, slack variables indicating load shedding and renewable generation curtailment are introduced to the gas/power balance equations, and the penalty costs are augmented to the objective functions accordingly.

Integrated Electricity-Gas System
It is assumed that the generators can be divided into two groups, i.e., CfU and GfU. Thus, we have G GfU = G\G CfU . Besides, we have G PtG = n∈N G PtG n for PtG facilities. The natural gas system and the power system are coupled via the following equations: where η GfU g , η PtG v are the efficiencies of GfU g and PtG v given by 0.43 and 0.58 respectively [14,33], and Hg is the heating rate of natural gas given by 1.08×10 4 MW/Mscm.
The coupling parameters are regarded as decision variables in IEGS, so it is necessary to add bounds for them, e.g., where p GfU g , p PtG v are the capacities of GfU g and PtG v, respectively.
Combining models (1), (2), coupling constraints (3) and the bounds of coupling variables (4), the integrated electric-gas system model (IEGS) can be obtained. For brevity, we denote by x the binary variables, by y the continuous variables, and by c I , c C the cost vectors associated with them. Eventually, IEGS can be written as, The only non-convex part in IEGS is the general flow equaiton (1i), which is represented by set Y GF in Problem (5).

Uncertainty Modeling
To address the variabilities and uncertainties of renewable energy resources and load demands, renewable generations as well as electricity/gas demands are viewed as random variables, and a stochastic-programming-based model is developed. Conventionally, stochastic programming relies on the probability distribution of random variables. In practice, however, the probability distribution may not exactly exist or the parameters cannot be obtained [34]. In recent years, non-parametric statistical methods have been introduced to the power and energy society [34,35], which help drawing an empirical distribution from historical data without the necessity of assuming any types of distribution for random variables.
In the proposed data-driven method, we first extract the forecast errors from historical data by subtracting the day-ahead forecast values from the real-time values, then use a scenario reduction method to select some representative error scenarios, and finally employ the reduced error scenarios to construct the scenarios by adding the errors to the day-ahead forecast value [36]. A Wassersteinmetric-based scenario reduction algorithm [37] is used for scenario reduction. The Wasserstein metric, also known as the Earth Mover's distance, is a function that defines how close two probability distributions are [38], and is more suited for measuring the distance of distributions than many other metrics such as the Euclidean distance. It is worth noting that the reduced scenario set obtained from this algorithm preserves the correlations between high-dimensional random variables [37]. It is worth to mention that many other techniques can be introduced to improve the statistical performance of scenario selections. For example, aside from probability metric methods, importance sampling, which aims at selecting scenarios that best represent the average cost impact of uncertainty on the problem [39], should be a promising alternative.
In what follows, each realization (scenario) of random nodal injections W n,t , D n,t and F Load s,t are denoted as ξ. Moreover, we denote by Ω the index set of ξ, ξω the ω-th scenario, yω the ω-th recourse variable, and σω the probability of the ω-th scenario. In two-stage stochastic programming, the second-stage recourse variable is a function of the first-stage decision and the random variable. Therefore, Y − is written as Y − (x, ξ), and the stochastic integrated electric-gas system model (S-IEGS) can be formulated as follows, in which the on/off statues of CfUs are optimised according to the reduced scenario set, and the second-stage dispatch decision regarding each scenario is determined accordingly. The price function in S-IEGS is assumed to be in line with that in IEGS. However, it is possible to formulate S-IEGS as a two-settlement process, i.e., attach the pre-dispatch quantity under the forecast scenario with price c C , and multiply the adjusted productions under each scenario with intra-day deviation penalties [40]. The reasons why stochastic programming is preferred in this paper to address the uncertainties in IEGS are threefold: 1)Existing works devoted to stochastic-programming-based cooptimisation problems of IEGS are still limited [11,15]. 2)As shown in Section 5 and 6.1, taking the advantage of stochastic programming, the solution procedure ends up iteratively solving some separable convex problems, the convergence and optimality of which are guaranteed. 3)Although cutting-edge techniques like (distributionally) robust optimisation can also deal with uncertainties, they make the MINLP problem rather complicated, so that approximation algorithms (not only for the physical model itself) [22] and heuristics [21] become inevitable.
To support the viewpoints above, distributionally robust optimisation is adopted for comparison. The distributionally robust integrated electric-gas system model (DR-IEGS) can be formulated as follows, where Ξ is the feasible region of ξ, P 0 (Ξ) denotes the set of all probability measures on a sigma algebra of Ξ, and the subset P is known as the ambiguity set in distributionally robust optimisation [22,28]. For tractability, only linear moment constraints are considered in the ambiguity set as in Ref. [22,41]. It should be noted that linear moment is not capable of modeling the correlation of uncertainties. The distributionally robust model (7) is also a data-driven approach. Historical data is used to construct the ambiguity set, among which the model aims to seek a worst-case distribution. The main difference between S-IEGS and DR-IEGS is that the optimal decision is derived based on the generated scenarios in Problem (6), whereas the optimal decision is achieved regarding the worst-case distribution in Problem (7).
Assuming that historical data is available to both S-IEGS and DR-IEGS, one can show that the stochastic model has the advantage over the distributionally robust model in terms of tractability and in-sample/out-of-sample performances. These will be demonstrated in Section 5 and Section 6.

Pricing PtGs in Day-ahead Market under Uncertainties
The main role that PtGs play in the integrated system is to consume surplus renewable generations and produce natural gas. Therefore, the contributions of PtGs are twofold: i) reducing the penalty cost (or the environmental cost) of renewable generation curtailments, and ii) supplying additional natural gas. It is necessary to quantify such contributions, especially in a competitive market. One common method is using the locational marginal prices (LMPs), which are the optimal Lagrangian multipliers of the optimisation problem that determine the costs of producing one extra unit of resource at different loccations [42]. If we associate with the gas balance equation (1n) and the power flow equation (2h) Lagrangian multipliers λ and µ respectively, then the "net" LMP (or LMP simply) of PtG v that defined in $/MW from the PtG's perspective is where λ n,t is the multiplier of (1n) for node n at time t, and µ m,t is the multiplier of (2h) for bus m at time t. Both λ n,t and µ n,t can be either positive, zero, or negative. It can be proved that when Problem (5) is solved to optimality [Since Problem (5) is a MINLP, solving it to optimality is defined herein as: fixing the binary variables as their optima, and re-solving the NLP problem to optimality (maybe local optimality) to obtain the optimal multipliers.]: i) the electric power consumed by PtG v is non-zero if and only if ψ v,t is non-positive; ii) ψ v,t is negative if and only if the capacity of PtG is inadequate. The former holds since otherwise the conversion would increase the total cost. The latter holds since otherwise the PtG production level can be improved to further reduce the total cost, which is contradictory with the fact that the current solution is optimal. The second observation suggests that PtG can only profit from congestion under the LMP-based pricing regime.
Evidently, the above-mentioned LMP only reflects the marginal value of PtG under a certain scenario (i.e., the forecast scenario), and it doesn't accounts for the flexibility service that PtG could provide after the realization of uncertainty. Due to the significant randomness in day-ahead markets, it is crucial to price the flexible resources provided by PtGs that mitigate the uncertainties [27,28]. As such, a novel concept of expected locational marginal value (E-LMV) is proposed in this paper. E-LMV can be formed with the byproduct of solving S-IEGS: where the subscript ω of ψ v,t,ω and p PtG v,t,ω indicates that they are derived from the ω-th scenario. Intuitively, E-LMV represents the expectation of payment that is entitled to PtG, regarding its potential recourse actions after uncertainties reveal. By taking the expectation value of multiple LMPs, E-LMV provides a payment scheme that is closer to the "true" (in terms of mathematical expectation) LMP, and therefore is suited for a market with considerable uncertainties.
E-LMVs can be defined similarly for the other participants in the day-ahead market. For example, E-LMV of RES at bus n at time t is given by Ultimately, the day-ahead market is settled based on E-LMVs. We have the following proposition for E-LMVs (see Appendix for the proof and further discussions), which suggests that the money collected by system operators from consumers is more than that should be paid to suppliers.

Proposition 1.
Supposing there is no gas compressors in the natural gas system, E-LMVs ensure revenue adequacy for the integrated electric-gas system.
For the distributionally robust model, we propose using the extremal distribution P * to derive E-LMV. Since a distributionally robust optimisation problem always possesses a discrete extremal distribution, E P * [Ψ v,t ] can be calculated using the extremal distribution as in Eqn. (9). Intuitively, E-LMV yielded from the distributionally robust model should be higher than that from the stochastic model; this will be verified in Section 6.

Solution Algorithm
In this section, we first introduce a method to address the nonconvexity issue for the natural gas flow model, and then present the overall solution algorithm for S-IEGS.

Convexification of Nonlinear General Flow Equation
The most challenging part of Problem (6) is the non-convexity of the general flow equation, as detailed in Eqn. (1i). Techniques for tackling this difficulty can be divided into: i) nonlinear programming (NLP) methods that solve the problem with interior point methods, etc.; ii) MILP reformulation and second-order cone programming (SOCP) approximation [14,18] that aim to approximate with high accuracy the non-convex problem using tractable mathematical programmings; and iii) intelligent algorithms like particle swarm optimisation, genetic algorithm, and neural networks [43][44][45][46]. Noting that SOCP approximation enjoys higher computational efficiency, and it is such that decomposition methods could be easily implemented, we adopt it in this paper.
The first row in Eqn. (11) defines a second-order cone: Q 1 (m,n),t = π m,t ,f Pipe (m,n),t ,π n,t π m,t ≥ f Pipe (m,n),t ;π n,t , whereπ m,t = K gf (m,n) π m,t andπ n,t = K gf (m,n) π n,t . The second row in Eqn. (11) results in a DC (difference of convex functions) programming that is difficult to solve in general. According to [47], DC programming can be approximately solved by a penalty convexconcave procedure (PCC). Specifically, the concave items are linearised at the current points, yielding a convex problem (SOCP in this paper; hence a sequential SOCP method), and then sequentially, the convex problem is solved to update the points for linearization. To ensure feasibility, a positive slack variable is needed: The intersection of Q 1 (m,n),t and Q 2 (m,n),t equivalently forms the feasible set of constraint (11) only if s + (m,n),t vanishes. For brevity, we define for each scenario the convex approximation of Y GF as Q, which is given by Q = y π m,t ,f Pipe (m,n),t ,π n,t ∈ Q 1 (m,n),t ∩ Q 2 (m,n),t ∀(m, n) ∈ G Pipe , t ∈ T .
Moreover, for ease of exposition, a normalised slack variable is defined ass + , the [(m, n), t]-th entry of which is given by s + (m,n),t / (K gf 2 (m,n) π 2 m,t ).

Generalised Benders Decomposition with PCC
Incorporating PCC into the generalised Benders decomposition procedure, an algorithm for solving S-IEGS can be obtained, as detailed in Algorithm 1. The Benders sub-problem is modified to avoid the necessity of solving a dual SOCP problem. Specifically, by introducing equality constraint (13a), it can be proved that the optimal dual variable associated with this constraint, which is available from offthe-shelf solvers, is sufficient to construct a Benders cut. Besides, in order for a valid cut, strong duality must hold for the Benders sub-problem, which in turn requires that Problem (13) and its dual have strictly feasible solutions, i.e., Q 1 (m,n),t and Q 2 (m,n),t have non-empty interior [48]. In computational practice, the feasibility condition is ensured by introducing slack variables to the power/gas balance equations (and penalty costs to the objective function accordingly), while the non-empty interior condition is guaranteed by the slack variable of Q 2 (m,n),t . The proposed algorithm has some desirable properties: 1)It is separable with regard to each scenario, and hence Problem (13) can be solved in parallel. 2)Problem (13) or its separated problems are convex, and thus can be solved in a decentralised manner by the electric system operator and the natural gas system operator, e.g., using ADMM. 3)If at each Benders iteration, Problem (13) can be distributedly solved to optimality, since MILP (12) involves merely the electric power repeat 6: Solve the current approximation, Problem (13)  Fig. 1: Flowchart of the solution algorithm.
system model, then, without confidential information of each system being revealed, Algorithm 1 converges and returns a UC solution.
For ease of reading, the framework of the whole solution algorithm is provided, which is shown in Fig. 1. The outer loop of the algorithm is the generalised Benders decomposition that iterates from the MILP master problem and the convex sub-problem. The Benders sub-problem is parallelizable, which means |Ω| scenarios could be addressed with PCC meanwhile. As mentioned above, the convex sub-problems can be decomposed into a linear programming (LP) of the power system dispatch problem and an SOCP of the gas flow problem, and then coordinated with ADMM.

Solution Method for Distributionally Robust Model
For comparison purpose, the distributionally robust model will also be solved. Yet, the convexification method and PCC algorithm cannot be easily extended to the distributionally robust model. One reason is that the convexified model is nonlinear, and thus the state-of-the-art method, linear decision rule (LDR) is inapplicable [22,41]. Another obstacle is that if we choose fully adaptive recourse instead of LDR, then the solution procedure requires dualizing the second-stage problem, making it unclear how to sequentially penalize the primal constraints.
To this end, Taylor series expansion is applied to linearize Eqn. (1i) for the distributionally robust model [22]. Although the linearised model is favorable for developing solution algorithm, it is less tight than the sequential SOCP method. Due to the abovementioned limitations, it is recognised that distributionally robust optimisation is not so attractive to the already complicated nonconvex IEGS problem.
The distributionally robust model is solved by an extremal distribution generation method proposed in [49]. The by-product of the solution method is an extremal distribution, which is then used for E-LMV calculation.

Case Studies
In this section, numerical experiments are carried out to validate the effectiveness of the proposed stochastic-programming-based model, the efficiency of the proposed solution method, and the advantage of the pricing method.
The test system is obtained by combining the IEEE 39-bus system and the Belgium 20-node gas system. The configuration of the integrated system is exactly as shown in Fig. 2 [33], and detailed data is available online [32]. Two 1200-MW wind farms are located at Bus 32 and Bus 33, resulting in a wind power penetration rate of 24.6%. In order to hedge against the volatile wind power generation and help consuming extra wind power, two 200-MW PtGs are installed near the wind farms, and the gas is injected into Node 13 and Node 14 of the gas system, respectively. The GfUs located at Bus 30, Bus 36, and Bus 37 are supplied by the gas extracted from Node 5, Node 2, and Node 13, respectively.
The day-ahead forecast and real-time data series of wind farm outputs and load demands over one year are adopted [50]. After scaling, we generate error scenarios with 85% of the data series (the dayahead forecast errors of wind power and load demands are assumed to be ±50% and ±10% respectively), and randomly remain 15% of them for out-of-sample tests. According to current practise, the penalty costs of wind curtailment and electric/gas load shedding are set higher in order to mimic the environmental cost, and reduce the loss of load, respectively. Without loss of generality, in the case studies, the price of wind curtailment is set to 10 times of the mean cost of power generation in the test system, namely 142 $/MWh; the prices of electric load shedding and gas load shedding are set to 200 times of the mean cost of power generation and the mean gas price in the test systems, namely 2840 $/MWh and 396 $/MBTU, respectively. The optimisation problems are built in GAMS 26.1.0 and solved by CPLEX 12.8. The relative convergence tolerance of CPLEX and those in Algorihtm 1 are all set as 10 -4 . All runs are executed on an Intel i5 CPU machine running at 1.80 GHz with 8 GB of RAM.

Performances of Proposed Algorithm
The efficiencies of the proposed algorithm is verified on multiple cases. The Benders loop converges with predefined accuracy (i.e., 10 -4 , and it converges to a zero gap in some cases) after 42 to 66 iterations. The PCC loop takes about 16 iterations, and the slack variables in Q 2 (m,n),t usually vanish (see Fig. 3), indicating that the solution is feasible to the primal MINLP. Despite being less computationally expensive, the linearised model used by DR-IEGS always produces non-zero residuals of the relaxed gas flow equations. So DR-IEGS seldom achieves a feasible solution to the primal MINLP, as also reported in [22].
The accuracy of Algorithm 1 is demonstrated via Table 1. For the nonlinear gas model, PCC finds a solution extremely close to the one returned by IPOPT, albeit it becomes more time-consuming due to a smaller step-size of the penalty factor (i.e., ς=1.02). For IEGS, Algorithm 1 finds a solution that is only 0.061% larger than the feasible solution returned by COUENNE, which exhaustedly runs out of time.
The total computational time of solving S-IEGS is reported in Table 2. Since the Benders sub-problem is separable, when leveraging parallel computations, the algorithm can actually terminate within 30 minutes even for the 100-scenario case (the average runtime of each scenario ranges from 952.57 seconds to 1727.08 seconds), thus meeting the time requirement of day-ahead markets. In order to test the scalability of the proposed algorithm, we replace the 39-bus system with the IEEE 118-bus system. Numerical results show that if we only impose power flow limits on critical transmission lines instead of all lines as in engineering practise, the S-IEGS problem is solvable within 2 hours accounting for the effect of parallel computation. Specifically, the relative gap of the Benders loop could be closed to about 10 -3 within 100 iterations, and PCC   Fig. 2: Diagram of integrated electricity-gas system (IEEE 39-bus system and Belgium 20-node gas system).
Converge to 0 at the last iteration   basically converges within 20 iterations. Although the number of iterations needed to solve S-IEGS is about 20 times (i.e., the average number PCC iterations) of that needed to solve a stochastic UC problem with similar scale, the overall computational effort turns out to be acceptable as the SOCPs could be solved quite efficiently. We also deploy the standard ADMM [51] to Problem (13), and find that the two-block SOCP can be solved to global optimality within 200 iterations, or solved to a 10 -4 gap within 20 iterations (see Fig. 4). The runtime of ADMM for the test system is several minutes. It is worth mentioning that in DR-IEGS, the sub-problem cannot be decomposed and precisely solved by the electric system operator and the natural gas system operator.
Therefore, the stochastic model and the proposed solution method is practicable and favorable in terms of efficiency, accuracy, scalability and the possibility of distributed computing.

Advantages of Proposed Stochastic Method
In this subsection, the improvement of UC decision brought by stochastic programming is evaluated. The benchmarks include a deterministic IEGS model (D-IEGS), which deals with uncertainties by operational reserves (the reserve rates for the gas system and the power system are 5% and 10% respectively), and the distributionally robust model described in Section 4.2.
The scenario reduction process is illustrated by Fig. 5. The left panel of Fig. 5 shows the 312 historical observations of wind power forecast error of two wind farms; the right panel of Fig. 5 shows the 20 reduced scenarios, in which a scenario with higher probability is plotted with a heavier line. By using the algorithm in [38], the Wasserstein distance between the reduced scenario set and the original data can be approximated. As shown in Table 3, the asymptotics of the reduced scenario sets is quite obvious, i.e., the distribution gets closer to the empirical one as the scenario size grows. As expected, the UC solution varies with the scenario size, and it "converges" as the number of scenarios becomes sufficiently large (see Table 3). In fact, only two "sub-optimal" UC solutions occur, which have distinct on/off statues over 10 or 1 time slots compared with the "optimal" one. We find that 20 scenarios might be representative enough for this case. It is observed that the extremal distribution yielded from DR-IEGS is quite "far" from the empirical distribution, and the UC solution also differs a lot with those of S-IEGS. For all the methods, after a UC decision is derived, in-sample and out-of-sample simulations are carried out to yield the expected costs under this UC solution. The simulation results for all the methods are presented in Fig. 6. The stochastic model outperforms the deterministic one slightly in terms of the amount of wind curtailments. As shown in Table 4, although the stochastic model incurs wind curtailments in the scheduling phase, the UC solution derived from it does reduce 2.17-MWh wind curtailments in simulations. Thus, the proposed method facilitates the utilization of wind power more effectively than the comparative decision making methods, and helps reducing the impact of greenhouse gas emission better. The cost saving achieved by optimizing the UC decision is about 0.12‰. The distributionally robust model minimizes the expectation of scheduling cost under the worst-case distribution, and thus the objective value and the wind curtailment level in the scheduling phase are both highest. The UC decision yielded is robust against the worstcase distribution, and results in less wind curtailment in real-time operation (see the last panel in Fig. 6). However, since the worst-case distribution rarely occurs, the UC solution is somewhat conservative and pessimistic. As shown in Table 4, the simulation cost for the distributionally robust model is highest, regardless of the lowest wind curtailment level. Another reason for the conservativeness is that the ambiguity set of DR-IEGS fails to model the correlation of random variables, and the extremal distribution contains many fast ramping events that are unlikely to occurs in reality. Although high-order moments can capture spatial and temporal correlations, incorporating them to DR-IEGS will give rise to some semidefinite programmings and bi-convex programmings, making the model more difficult to solve [52].
Throughout the computational experiment, load shedding doesn't occur in IEGS. This should be owed to the flexibility originating from gas storage stations and the line-pack effect.

Settlement of PtGs using E-LMV
To settle the day-ahead market, as usual, the UC solution yielded from S-IEGS is fed back to the deterministic model to obtain a predispatch solution and LMPs. In this way, the PtG production levels as optimally scheduled are presented in Fig. 7, together with LMPs defined by Eqn. (8). In the test system, LMPs of the power system range from 4.32 $/MW to 17.32 $/MW, while those of the natural gas system range from 7.15 $/MW to 7.45 $/MW (considering the efficiency factor, it is 4.14 $/MW to 4.32 $/MW). According to Fig. 7, PtGs convert power to gas only when ψ is zero, that is, the LMPs on the power system side and the gas system side all equal 4.32 $/MW. This verifies the claim in Section 4.3.
Noting that the minimum variable cost of generators is 10 $/MW [32], LMPs take the value of 4.32 $/MW only when the wind farms encounter overproduction. However, when overproduction occurs, absent PtGs, the LMPs of such buses would be non-positive. Therefore, it is easy to see that PtGs consume excess wind power, raise up the price, and end up getting less payment and often zero payoff. In this case, the payment to PtGs derived from LMPs is 0 k$, because congestion doesn't occur near Bus 32 and 33 under the forecast scenario.
If the market is settled using E-LMVs, the payments of PtGs at each time period are as shown in Fig. 8. In Fig. 8, the day-ahead forecast and the upper/lower envelop are also plotted. The envelope is obtained by taking the pointwise maximum/minimum of wind power levels in the scenario set, so it indicates the highest/lowest possible wind power level in the stochastic model. In this test system, the wind power capacity is 1200 MW, i.e., 100-MW higher than the summation of the PtG capacity and the transmission line capacity. Therefore, congestion occurs either when CfU at Bus 32 or 33 is scheduled OFF and the wind power exceeds 1100 MW, or when CfU at Bus 32 or 33 is scheduled ON and the wind power exceeds 900 MW or 950 MW (subtracting the minimum production level of the CfU). In the stochastic model, it is hard to seek a UC solution that incurs no congestion under all probabilistic scenarios. Therefore, payments to PtGs are more likely to occur. It can be seen from Fig. 8 that payments occur even when the highest possible wind power level is less than 1100 MW, because congestion exists under some scenarios given the optimal UC solution. However, if the payment is derived from the forecast value (i.e., the expected scenario), the payment is zero as above-mentioned. Therefore, E-LMV better reflects the expected value of PtGs than LMP of the expected scenario does. It can be expected that under a same system configuration, the more volatile and uncertain wind power is, the higher E-LMV will be.
The total credit to PtGs derived from S-IEGS is 4.03 k$. The value obtained from DR-IEGS is 27.82 k$, which is several times higher than that from S-IEGS. In fact, it may not be persuasive to settle the market based on the worst-case situation.
As defined in Eqn. (9), the ω-th scenario contributed to E[Ψ v,t ] only if ψ v,t,ω is negative, which requires that p PtG v,t,ω = p PtG v,t . Hence, the mechanism of the proposed settlement scheme is akin to the financial transmission right, but in a stochastic setting. According to Proposition 1, the payment to PtGs is balanced by the charge from volatile renewable generations and demands. The payment received by PtG owners can be spent on capacity expansion.

Long-term Marginal Value of PtGs
Using the same setting, we solve S-IEGS and run simulations for cases with different PtG capacities to assess the long-run contribution of PtGs.
According to Table 5, the marginal value of installing 100 extra MW of PtGs is remarkable when the initial capacity is 100 MW, which is given by the difference of expected costs, i.e., (3, 085.37 − Fig. 8: Payments to PtGs with respect to wind power (intervals). Although the PtG technology is still costly, the cost saving achieved by installing such facilities can be much higher than that via optimal scheduling only (4.00‰ v.s. 0.12‰). For IEGS, it is of vital importance to decide an economic PtG size. From this perspective, the results in Table 5 also suggest the applicability of S-IEGS model and the proposed algorithm to optimally sizing PtG capacities.

Conclusions and Discussions
In this paper, a data-driven stochastic model is developed to cooptimise IEGS in day-ahead markets and address multiple correlated uncertainties. The data-driven stochastic model has cost benefit compared with a deterministic model. Moreover, it is demonstrated that the stochastic model has advantage over a distributionally robust model in terms of algorithmic tractability, and also on cost efficiency due to the fact that the stochastic programming framework allows more precise modeling of the gas flow problem. The proposed algorithm ensures convergence and provides highquality solutions to the original MINLP problem, even under a decentralised computational setting. The computational time is reasonable regarding the clearing time of day-ahead market, as the algorithm framework allows parallel and distributed computing.
According to the analysis of LMPs at coupling buses/nodes, cost recovery is difficult for PtGs under a deterministic-LMP-based regime. The expected locational marginal value proposed in this paper provides an alternative to pricing PtG facilities in a day-ahead market with production and demand uncertainties, and it ensures that PtGs get sufficient payments to expand their capacities to better mitigate the volatile renewable generations. It is also demonstrated that the cost saving achieved by installing PtGs is higher than that via optimal scheduling.
The conclusion that λ f Pipe ≤ 0 now can be drawn based on Eqn. (23) and (24). The revenue, if exists, is caused by the limitations of flow rate and line-pack capacity, which are determined by Π,

Discussions
In fact, it is due to the enforced nodal pressures/flow rates, instead of the gas loss, that the revenue adequacy of gas market cannot be verified when compressors exist. This is similar to the electricity market. For example, if the rate of power flow on a transmission line is enforced to be higher than some levels, then costly power may flows to less-expensive locations, and the revenue adequacy of electricity market is not guaranteed.
For a nonlinear gas market with gas compressors, the revenue might still be non-negative in reality though it cannot be verified in theory. This is in line with the observation in numerical experiments, i.e., with the optimal Lagrangian multipliers of the SOCP model, one of the source nodes (Node 8) has a lower gas price than those at demand nodes, guaranteeing that the cost of the gas consumed by compressors can be compensated precisely (the revenue adequacy of the gas market is zero).