1 Introduction

The importance of electric gas units (mainly combined cycle gas turbines, CCGTs) has been increased during the recent years due to their attractive features as flexible resources in view of the increasing renewable energy sources (RES) penetration. This is attributed to their increased efficiency, small environmental footprint and high flexibility.

The popularity of CCGTs has given rise to common energy infrastructure considerations by regulators and system analysts, identifying the strong interdependence between the electricity and gas system in technical, economic and operational terms.

Despite their common nature as energy transmission systems, the operation of the natural gas system is extremely complex, employing a large-scale, highly nonconvex and nonlinear problem structure (comprising a group of nonlinear algebraic equations), which can be modeled as a nonconvex mixed-integer nonlinear program (MINLP) [1]. Compared to this, the electricity problem, which is a unit commitment problem with optimization of energy and reserves, can be modeled as a mixed-integrater program (MIP) [2].

A review of the developed nonlinear models for combined consideration of the electricity and gas systems is given in [3] and [4]. In all models, there is a nonlinear coupling constraint linking the electricity system with the gas system.

The majority of the research in the gas system focuses on simplifications of the nonlinear equations that govern the physics of gas flow as well as the general representation of the technical components (e.g. compressor model) of the gas transmission system. A current useful and efficient way to alleviate the burden of the nonconvex MINLP gas problem is the application of mixed integer linear programming techniques for piecewise linearization. In this way, the nonlinear formulation is approached by approximated linearized functions that render the problem as a MIP, which enables the computation of global optima within fast execution times. The most comprehensive article that studies the advantages and drawbacks of various MIP formulations for the gas optimization problem is [5], where different piecewise linearization approaches are tested for different test cases in both dynamic and stationary conditions in order to derive conclusions about the more accurate and faster approximations. Based on the results, the authors argue about the effectiveness of the incremental method and further apply it in [6] and [7] in order to linearize the transient flow equations for a combined electricity-gas model that is applied on the Belgian high-calorific gas network, where the effects of initial linepack in system operation are examined. A simpler linearization approach can be seen in [8] and [9] for a stochastic unit commitment and a combined optimal power flow model, respectively. Simple electricity-gas conversion factors are applied to all studies and the compressor machine modeling is also reduced to a compression factor. The incremental method is again applied in [10] for stationary gas flow within a large decomposition framework.

A literature overview of the linearized combined models discussed above is given in Table 1, according to the linearization approach applied, the modeling of the coupling constraint and the compressor, as well as the configuration of the gas flow. The main contributions introduced in this paper with respect to the current literature are the following.

  1. 1)

    The extended incremental method is used for the linearization of the gas physics [11]. This linearization method builds on the simple incremental method, exhibiting enhanced features and advantages as compared to standard piecewise linear approximation methods.

  2. 2)

    The nonlinear and nonconvex compressor operating range is linearized by an outer approximation approach. This approach often yields a rather detailed approximation (as opposed to simple linear approximations typically used in the literature), while keeping the computational overhead relatively low compared to the piecewise linearization of multivariate nonlinearities with the extended incremental method.

  3. 3)

    The coupling constraints (originally quadratic functions) are linearized with the extended incremental method as well, whereas in the literature a simple linear conversion ratio is applied.

Table 1 Literature overview for linearized models
Full size table

Section 2 describes the nonconvex MINLP, whereas Sect. 3 presents the corresponding linearized model (MIP). In Sect. 4, several computations are presented and discussed afterwards. Lastly, in Sect. 5 the computational results are summurized.

2 Problem statement

The nonconvex MINLP is composed of an electricity (MIP) and a gas model (nonconvex MINLP) which are described in Sects. 2.1 and 2.2, respectively. The coupling of these models via nonlinear functions is explained in Sect. 2.3.

2.1 Electricity problem

In this paper, the day-ahead scheduling performed by the transmission system operator involves the solution of the unit commitment problem with simultaneous optimization of energy and reserves. The problem objective function concerns the minimization of the total as-bid cost for energy and for the provision of reserves plus the unit start-up/shut-down costs, subject to the following system operational constraints and unit technical constraints, for all dispatch periods \(t\in T\) (typically, the dispatch period is 1 h) [2].

  1. 1)

    Unit start-up/desynchronization constraints.

  2. 2)

    Unit minimum up/down time constraints.

  3. 3)

    Logical status of unit commitment constraints.

  4. 4)

    Unit power output normal and AGC limits.

  5. 5)

    Unit reserve capability constraints.

  6. 6)

    Unit ramp-up/down constraints.

  7. 7)

    The power balance equations.

  8. 8)

    System reserve requirements.

Due to the presence of binary variables for unit commitment and automatic generation control (so-called AGC) status, the electricity problem constitutes a MIP model. The objective value is defined as the electricity system cost \(C_{\mathrm{e}}\) for all dispatch periods \(t\in T\).

2.2 Gas problem

The optimization of gas transport networks has been recently described in [1] and in the survey article [12]. The main objective is to find cost optimal operations, which stick to the required pressure and flow bounds. As gas mainly flows from higher to lower pressures, it is required to install compressor machines in order to increase the gas pressure again. These compressions can be controlled by the transmission system operator, which adds discrete aspects to the problem among others. The physics of gas transport are driven by differential equations, which are simplified herein for the description of a steady-state problem. In this subsection, the gas model is analytically described. It should be noted that in all descriptions of Sect. 2.2 the dispatch period index is deliberately ignored.

The gas network is modeled via a directed graph \(G=(V,A)\) with node set \(V\) and arc set \(A\). The set of gas nodes \(u,v\in V\) is partitioned into the set of entry nodes \(V_+\) (where gas is supplied), the set of exit nodes \(V_-\) (where gas is discharged) and the set of inner nodes \(V_0\) (without any supply or discharge). Entries can be LNG (liquefied natural gas) entries \(V_{\text {lng}}\) or standard pipeline entries \(V_{\text {ent}}\). Exits can be electric gas loads (power plants) \(V_{\text {e}}\) or non-electric gas loads \(V_{\text {ne}}\). The arcs are divided into pipes \(A_{\text {pi}}\), control valves \(A_{\text {cv}}\) and compressor machines \(A_{\text {cm}}\).

For every arc \(a= (u, v)\in A\), variables \(q_{a}^{}\) are defined, expressing the gas (mass) flow in arc direction when \(q_{a}^{} > 0\) and vice versa, and having certain lower and upper bounds.

$$\begin{aligned} q_{a}^{} \in [q_{a}^{-}, q_{a}^{+}] \qquad \forall a\in A\end{aligned}$$
(1)

Every node \(u\in V\) has ingoing and outgoing arcs, defined as \(\delta ^{\text {in}}(u)\) and \(\delta ^{\text {out}}(u)\), respectively. For every node \(u\in V\), a pressure variable \(p_{u}^{}\) and a flow variable \(q_{u}^{}\) is defined with certain upper and lower bounds.

$$\begin{aligned} {\left\{ \begin{array}{ll} p_{u}^{} \in [p_{u}^{-}, p_{u}^{+}] &{} \quad \forall u\in V\\ q_{u}^{} \in [q_{u}^{-}, q_{u}^{+}] &{} \quad \forall u\in V\end{array}\right. } \end{aligned}$$
(2)

The flow variables \(q_{u}^{}\) represent the gas inputs (nonpositive) and outputs (nonnegative) at node \(u\). Flow bounds are implicitly computed from pressure bounds and are useful only for preproccessing purposes. The nodal flow variables are fixed to zero for inner nodes.

$$\begin{aligned} q_{u}^{-} = q_{u}^{+} = 0 \qquad \forall u\in V_0\end{aligned}$$
(3)

Moreover, there is a constraint for the outputs at non-electric gas loads. The forecasted gas demand \(d_{u}^{}\) for every \(u\in V_{\text {ne}}\) should be fulfilled as far as possible. If this is not possible, a so-called shedding variable \(s_{u}^{}\) measures this deviation according to the following:

$$\begin{aligned} q_{u}^{} + s_{u}^{} = d_{u}^{} \qquad \forall u\in V_{\text {ne}}\end{aligned}$$
(4)

Using the flow variables, the mass flow conservation equation is formulated as follows:

$$\begin{aligned} \sum _{a\in \delta ^{\text {in}}(u)} q_{a}^{} - \sum _{a\in \delta ^{\text {out}}(u)} q_{a}^{} = q_{u}^{}\qquad \forall u\in V\end{aligned}$$
(5)

The pipes \(A_{\text {pi}}\) are used to transport gas through the network and are specified by their length \(L_{a}\), diameter \(D_{a}\), and integral roughness \(k_a\). Here the differential equations of gas dynamics come into play. These gas dynamics within a single pipe are described by the Euler equations [13]. However, finishing well-known simplification steps, an analytic formula of pressure and flow in a pipe can be given via the Weymouth equation [1], Chapter 6.

$$\begin{aligned} p_{v}^{2} - p_{u}^{2} = \frac{L_{a} \lambda _{a} R_\text {s}z_a T}{A_{a}^2D_{a}}|q_{a}^{}|q_{a}^{}\qquad \forall a\in A_{\text {pi}}\end{aligned}$$
(6)

The specific gas constant \(R_\text {s}\) depends on the gas composition, while the friction factor \(\lambda _{a}\) is calculated from the diameter and the integral roughness using the formula of Nikuradse [14], and \(z_a\) is the average compressibility factor. Furthermore, the cross-sectional area of pipe \(a\in A_{\text {pi}}\) can be calculated with the help of the diameter as follows:

$$\begin{aligned} A_{a} = D_{a}^2 \frac{\pi }{2}\qquad \forall a\in A_{\text {pi}}\end{aligned}$$
(7)

\(z_a\) of pipe \(a\in A_{\text {pi}}\) is calculated with an average pressure on a pipe regarding the bounds, thus

$$\begin{aligned} p_{m, a}^{} = \frac{ \max \left\{ p_{u}^{-}, p_{v}^{-} \right\} + \min \left\{ p_{u}^{+}, p_{v}^{+} \right\} }{2}\qquad \forall a\in A_{\text {pi}}\end{aligned}$$
(8)

The formula of the American Gas Association is used to calculate the compressibility factor [15]:

$$\begin{aligned} z(p, T) = 1 + \alpha p\end{aligned}$$
(9)

where \(\alpha\) constitutes a parameter given as a function of the pseudocritical pressure \(p_c\) and the pseudocritical temperature \(T_c\):

$$\begin{aligned} \alpha = 0.257 \frac{1}{p_c} - 0.533 \frac{T_c}{p_cT} \end{aligned}$$
(10)

The isothermal case assumes constant temperature \(T\), which is used for that model. Now, (8)–(10) are combined for the calculation of \(z_a\). \(z_a = z(p_{m, a}^{}, T), \forall a\in A_{\text {pi}}.\)

Moreover, so-called short pipes are defined as auxiliary network elements. They are used to model multiple entries at a single entry node or to model multiple exits at a single exit node and can be interpreted as pipes \(a= (u, v)\in A_{\text {pi}}\) with length \(L_{a} = 0\) and therefore from (6), \(p_{u}^{} = p_{v}^{}\) can be concluded.

Control valves \(A_{\text {cv}}\) are used to decrease gas pressure. These elements are modeled with some discrete aspects, since they can be operated in different modes (active, bypass and closed). A closed control valve decouples the in- and outflow pressures and prohibits any gas flow through \(a= (u, v)\in A_{\text {cv}}\). The bypass mode ensures the flow through the control valve, but does not influence pressure, while the active mode ensures the desired pressure decrease within a controllable range \(\Delta _{a}^{} \in [\Delta _{a}^{-},\Delta _{a}^{+}]\):

$$\left\{ {\begin{array}{*{20}l} {a\;{\text{is}}\;{\text{active}}} \hfill & { \Rightarrow \;p_{v} = p_{u} - \Delta _{a} } \hfill \\ {a\;{\text{is}}\;{\text{in}}\;{\text{bypass}}\;{\text{mode}}} \hfill & { \Rightarrow \;p_{v} = p_{u} } \hfill \\ {a\;{\text{is}}\;{\text{closed}}} \hfill & { \Rightarrow \;q_{a} = 0} \hfill \\ \end{array} } \right.$$
(11)

Compressor machines \(A_{\text {cm}}\) are necessary to increase pressure, in order to transport gas over long distances. As control valves, a compressor machine can be either active, closed or in bypass mode. The working ranges of turbo compressors are typically described by characteristic diagrams, as shown for example in Fig. 1, where the admissible operating range is depicted in green.

Fig. 1
figure 1

Characteristic diagram of turbo compressor

Full size image

The horizontal axis in the diagram is labeled with the volumetric flow rate at inlet conditions.

$$\begin{aligned} Q_{a}^{}=\rho _{u}^{}q_{a}^{} \end{aligned}$$
(12)

where \(\rho _{u}^{}\) denotes the density at node \(u\). Density is related to pressure via the equation of state:

$$\begin{aligned} p_{u}^{} = \rho _{u}^{} R_\text {s}T_{u}^{} z_u \end{aligned}$$
(13)

In the isothermal model presented in this paper, the gas temperature at the inlet \(T_{u}^{}\) is constant. The vertical axis is labeled with the specific change in adiabatic enthalpy

$$\begin{aligned} H^\mathrm {ad}_{a} = R_\text {s}T_{u}^{} z_u \frac{\kappa }{\kappa - 1} \left[ \left( \frac{p_{v}^{}}{p_{u}^{}}\right) ^{\frac{\kappa - 1}{\kappa }} - 1 \right] \end{aligned}$$
(14)

The isentropic exponent \(\kappa =1.38\) is also a constant in the model.

In Fig. 1, the dashed lines represent isolines for the adiabatic efficiency of the compressor and the thin solid lines are isolines for compressor speed. The upper and lower bounding curve of the feasible region is given by the isolines for minimum and maximum speed. To the left, the working range is bounded by the surgeline. The right boundary curve is called chokeline. Both surgeline and chokeline are depicted with thick solid lines. All these curves result from (bi-)quadratic least squares fits of the measured data points, which are depicted as crosses in Fig. 1. The four bounding can be written as nonlinear inequalities of the form:

$$\begin{aligned} f^i_{a}(Q_{a}^{},H^\mathrm {ad}_{a}) \le 0 \qquad \forall a= (u, v)\in A_{\text {cm}},\quad i=1, 2, 3, 4 \end{aligned}$$
(15)

The power required for compression is given by

$$\begin{aligned} P_{a}^{} = \frac{q_{a}^{}H^\mathrm {ad}_{a}}{\eta _{a}^{}}\qquad \forall a= (u, v)\in A_{\text {cm}}\end{aligned}$$
(16)

where, \(\eta _{a}^{}\) denotes the adiabatic efficiency of the compressor, which is considered constant in this paper. Further details on the modeling of compressors can be found in [1], Sect. 2.

Since the objective is to minimize costs, the combination of the compressor costs and the shedding costs are calculated. It should be noted that the gas supply cost of the power plants is included in their variable cost, which is implicitly incorporated within their energy offer prices in the electricity model; thus, they are not included again in the objective function of the gas optimization problem. The compressor costs depend on the power of compressor machine \(P_{a}^{}\), while the shedding costs depend on \(s_{u}^{}\). The corresponding parameters \(c_{a}\) for \(P_{a}^{}\) and \(\gamma\) for \(s_{u}^{}\) represent the cost coefficient for compressor machine \(a\in A_{\text {cm}}\) and the cost term (penalty price) for shedded gas flow, respectively, and are used to compute the objective value. The overall optimization problem of the gas system is formulated as follows:

$$\begin{aligned} \min \quad&\sum _{a\in A_{\text {cm}}} c_{a} P_{a}^{} + \sum _{u\in V_{\text {ne}}} \gamma s_{u}^{} \end{aligned}$$
(17)
$${\text{s.t.}}\quad{\text{pressure}}\;{\text{and}}\;{\text{flow}}\;{\text{bounds}}\;({\text{1}}){-}({\text{3}})$$
(18)
$${\text{shedding}}\;{\text{at}}\;{\text{non-electric}}\;{\text{gas}}\;{\text{loads}}\;{\text{(4)}}$$
(19)
$${\text{mass}}\;{\text{flow}}\;{\text{conservation}}\;{\text{(5)}}$$
(20)
$${\text{pressure}}\;{\text{decrease}}\;{\text{in}}\;{\text{pipe}}\;{\text{(6)}}$$
(21)
$${\text{control}}\;{\text{valve}}\;{\text{model}}\;{\text{(11)}}$$
(22)
$${\text{compressor}}\;{\text{model}}\;({\text{12}}){{-}}({\text{16}})$$
(23)

The overall compressor costs \(C_{\mathrm{c}}\) and the overall costs due to shedded gas flows \(C_{\mathrm{s}}\) for all dispatch periods \(t\in T\) can now be formulated with \(C_{\mathrm{c}}= \sum _{t\in T} \sum _{a\in A_{\text {cm}}} c_{a}P_{a}^{t}\) and \(C_{\mathrm{s}}= \sum _{t\in T} \sum _{u\in V_{\text {ne}}} \gamma s_{u}^{t}.\)

Thus (17)–(23) can be extended to all dispatch periods:

$$\begin{aligned} \min \quad&C_{\mathrm{c}}+ C_{\mathrm{s}}\end{aligned}$$
(24)
$${\text{s.t.}}\quad{\text{gas}}\;{\text{constraints}}\;({\text{18}}){-}({\text{23}})\;{\text{for}}\;{\text{all}}\;{t}\in {{T}}$$
(25)

Problem (24)–(25) yields a nonconvex MINLP, since the control valves and the compressor machines include discrete variables, and (6) and (12)–(16) include nonlinear functions as equations and nonconvex inequalities. Note that Problem (24)–(25) is block diagonal, where each block is made up of the variables and constraints of a single discrete time step.

2.3 Coupling of gas and electricity problem

The two energy systems are coupled through the gas consumption of the electric gas loads. The coupling constraint can be mathematically formulated with

$$\begin{aligned} u_{u}^{t} \frac{a_{u} + b_{u}e_{u}^{t} + c_{u}(e_{u}^{t})^2}{H} = q_{u}^{t} \qquad \forall u\in V_{\text {e}}, \forall t\in T\end{aligned}$$
(26)

The binary variable \(u_{u}^{t}\) and continuous variable \(e_{u}^{t}\) express the commitment status of electric gas load \(u\in V_{\text {e}}\) and the electricity production of unit \(u\in V_{\text {e}}\), at dispatch period \(t\in T\), respectively. The fuel consumption coefficients of electric gas loads \(u\in V_{\text {e}}\), which are \(a_{u}\), \(b_{u}\), and \(c_{u}\), as well as the heating value of natural gas \(H\), are needed to calculate the conversion. For simplicity, it is assumed that the composition of the supplied gas does not vary among the system entries.

Intertemporal constraints between dispatch periods are only present in the electricity unit commitment problem in Sect. 2.1, while the gas model (24)–(25) can be separated for every dispatch period. The overall optimization problem is expressed as a minimization problem where the respective costs of the two systems are combined to express the overall objective function. It is a nonconvex MINLP:

$$\begin{aligned} \min \quad&C_{\mathrm{e}}+ C_{\mathrm{c}}+ C_{\mathrm{s}}\end{aligned}$$
(27)
$${\text{s.t.}}\quad{\text{electricity}}\;{\text{constraints}}\;{\text{in}}\;{\text{chapter}}\;{\text{2.1}}$$
(28)
$${\text{gas}}\;{\text{constraints}}\;{\text{(25)}}$$
(29)
$${\text{coupling}}\;{\text{constraints}}\;{\text{(26)}}$$
(30)

3 Linearized formulation

Two linearization techniques are used to tackle the nonlinearities of (27)–(30). The nonlinear functions in (6), (16), and (26) are linearized with the so-called extended incremental method [11]. The nonlinear and nonconvex set of (15) is linearized with outer approximation constraints [16].

Both techniques ensure the relaxation property, which means that each solution to the nonlinear model (27)–(30) is also feasible for the linearized model. This is important as pure approximations tend to be infeasible even if the underlying nonlinear model is feasible [17].

3.1 Extended incremental method

First of all, it should be mentioned that the incremental method as basis of the extended incremental method is just one of various methods to linearize nonlinear functions. MIP-relaxations may be formulated with alternative methods, for instance with a convex combination model that introduces only a number of extra binary variables and constraints that is logarithmic in the number of breakpoints [18]. However, the incremental method is used in this paper, since it performs best for gas transport problems [5, 17, 19].

The extended incremental method works by first introducing a new finitely bounded variable for each nonlinear term, and by then computing a piecewise linear approximation of it. The approximation is constructed such that an a priori given upper bound on the approximation is satisfied, while introducing as less breakpoints as possible. Using the known values of the resulting approximations error, a MIP-relaxation model for the nonlinearities is derived.

To be more precise, for a constraint of the form \(\varvec{c}^\mathrm {T}\varvec{x} + f(y) = 0\), with \(f:\mathbf R \rightarrow \mathbf R\) nonlinear, and finitely bounded \(y \in [y_{\min },y_{\max }]\), a new variable z is introduced, reformulating the constraint as \(\varvec{c}^\mathrm {T}\varvec{x} + z = 0\), and a piecewise linear approximation \(\varPhi (y)=z\) of f(y) over \([y_{\min },y_{\max }]\) with \(\max _{y_{\min } \le y \le y_{\max }} |\varPhi (y) - f(y)| \le \varepsilon\) is computed.

Assuming that \(\varPhi (y)\) is defined over n intervals, subdivided by breakpoints \(y_{\min }=b_0< b_1< \dots <b_{n-1}=y_{\max }\), the overall MIP-relaxation model of the constraint \(\varvec{c}^\mathrm {T}\varvec{x} + f(y) = 0\) is reformulated as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} y = b_0 + \sum \limits _{i=1}^{n}(b_i - b_{i-1})\delta _i \\ z = f(b_0) + \sum \limits _{i=1}^{n}\left( f(b_i) - f(b_{i-1})\right) \delta _i + e \\ \delta _{i+1} \le w_i \le \delta _i \qquad i = 1, 2, \ldots , n-1 \\ \delta _{i} \in [0,1] \qquad \qquad i = 1, 2, \ldots , n\\ w_i \in \left\{ 0,1\right\} \qquad \quad i = 1, 2, \ldots , n-1\\ e \in [-\varepsilon , \varepsilon ]\\ z \in \mathbf R \end{array}\right. } \end{aligned}$$
(31)

Each of the \(\delta\)-variables is associated with a discretization interval. With the so-called filling condition \(\delta _{i+1} \le w_i \le \delta _i\), we make sure that \(\delta _i > 0\) requires \(\delta _j = 1\) for all \(j<i\) and that \(\delta _i < 1\) requires \(\delta _j = 0\) for all \(j>i\). This way, there can be at most one index i of an interval with \(0< \delta _i < 1\). In this case, y lies within the \(i^\mathrm{th}\) interval. Assuming \(e=0\), the point (y, z) would lie on the \(i^\mathrm{th}\) line segement of the piecewise linear function. However, with \(e \in [-\varepsilon , \varepsilon ]\), the point (y, z) lies within a box of height \(2\varepsilon\) that encloses the graph of the piecewise linear function. Since \(\varepsilon\) has been chosen such that also the graph of the approximated nonlinear function is contained in the boxes, model (31) yields a relaxation.

The extended incremental method is used for the nonlinear terms in (6), (16), and (26). In order to rewrite (16) as a sum of univariate nonlinear expressions, we plug in Eq. (14) and apply the standard reformulation \(xy = -((x - y)^2\,{-}\,x^2 - y^2)/2\) for bilinear expression two times, before the extended incremental method is applied. It should be also noted that that after linearization, (26) is still a product of a binary variable \(u_{u}^{t}\) and a newly defined continuous variable according to (31). This resulting nonlinear term can be linearized without an error by the well-known bigM method.

3.2 Outer approximation

Outer approximation is a well-known linearization method for convex nonlinear programs [16]. In this paper, it is used for the linearization of (15), although a nonconvex operating range is involved. The idea is to compute a convex envelope of the operating range and to add tangential hyperplanes to the envelope during optimization. More details on this technique can be found in [1].

4 Case study

The developed model is tested on the Greek power system together with the Greek natural gas transmission system. The description of these networks is given in Sect. 4.1, while the computational results are presented in Sect. 4.2. First, the linearized formulation of Sect. 3 is tested in Sect. 4.2.1. Afterwards, the results are compared with a single electricity model in Sect. 4.2.2, where the gas network is not taken into account.

4.1 Real-world test system

The Greek power system is operated by ADMIE (Independent Power Transmission Operator S.A.) [20]. The generation mix and the marginal cost range for each unit category is presented on Table 2. In the test case, each thermal unit is considered to submit one price/ quantity offer at its minimum variable operating cost. The opportunity cost of hydro units (cost of replacement of a thermal MWh) is considered for the generation of the respective energy offers. This has been considered slightly higher than the variable cost of the most expensive thermal units. Indeed, in countries with a low portfolio of hydro resources (e.g., in Greece), such opportunity cost is significantly high, since hydro units essentially operate as peakers (peak-shaving units), replacing expensive (high-cost) thermal units. Finally, the system load is considered inelastic.

Table 2 Greek power system generating units’ data with unit type, number of unit types, installed capacity, and marginal cost range
Full size table

The Greek natural gas transport network is operated by the Greek natural gas transmission operator DESFA [21]. Figure 2 shows the network topology and Table 3 gives some statistical information about the network.

Table 3 Data of Greek gas network with 134 nodes and 133 arcs
Full size table

The system constitutes a radial network consisting of three entry points, namely two pipeline entries in the northern part of the country in Sidirokastro and Kipi, and an LNG entry in the southern part, in the island of Revythousa close to the metropolitan area of Athens. The electric gas loads are mainly located in the south, while the system is supplemented by a compressor station in Nea Messimvria (close to Thessaloniki) to maintain gas flows from the northern entries. Moreover, a control valve is located at the Lavrio branch.

Fig. 2
figure 2

Greek gas network

Full size image

4.2 Computational results

All computations have been performed on a computer with an Intel i7-4770 CPU with 4 cores running at 3.4 GHz each, and 16 GB of main memory.

The electricity model was developed with GAMS 24.4.6 [22], while the linearized gas model was built with the C++ software Framework LaMaTTo++ [23]. GAMS was used to combine these models and to linearize the coupling. As MIP solver, Gurobi 6.0.4 was used [24]. A MIP gap tolerance of 0.01‰ was selected. Moreover, a time limit of 1 h was set as the electricity unit commitment problem has to be solved in 1 h.

As instances, the days with highest non-electric gas loads \(\sum _{t\in T}\sum _{u\in V_{\text {ne}}}d_{u}^{t}\) out of all days of every quarter of the years 2013–2015 were chosen yielding in total twelve instances.

Regarding the approximation accuracy, two varying parameters have been allowed. The first parameter \(\varepsilon _1\) is the maximum error of the approximation of the nonlinear pressure loss function (6) in bar. The second parameter \(\varepsilon _2\) is the maximum error allowed for \(e_{u}^{2}\) in \(\text {MW}^2\). Four different pairs of values for \(\varepsilon _1\) and \(\varepsilon _2\) are compared in the computations that are performed as shown in Table 4.

Table 4 Approximation accuracy levels
Full size table

Moreover, preprocessing according to chapter 6 in [1] is used for the gas problem. These methods tighten the bounds of flow variables and reduce the number of linearized functions. If a squared variable \(x^2\) occurs in the problem without the use of the variable x itself, then \(x^2\) can be used without a linearization and x can be calculated a posteriori.

Regarding the entries, two different scenarios were examined.

Case I For the first case, the historical gas entry distribution of Kipi, Sidirokastro and Revythousa from 2013 to 2015 was averaged. Therefore, Kipi and Revythousa provide upto 20% of the gas each and the remaining upto 60% are supplied at Sidirokastro. In the present model, the optimizer is allowed to deviate from this distribution by at most ten percent at each entry.

Case II Another case with closed LNG entry is implemented, due to its periodic maintenance or outage. The entry distribution of Kipi and Sidirokastro can be chosen freely by the optimizer in this case. However, it is required that at least 20% of consumed gas has to be supplied from Kipi.

4.2.1 Linearized formulation

The average solution time and problem size for each accuracy level from Table 4 is given in Table 5, while the day with the highest non-electric gas load, the 18/02/2015, is presented in Table 6 in detail.

Table 5 Average number of continuous variables, discrete variables, constraints, and average solution time. The numbers in parentheses denote the number of instances that hit the time limit
Full size table
Table 6 Electricity costs \(C_{\mathrm{e}}\), compressor costs \(C_{\mathrm{c}}\), shedding costs \(C_{\mathrm{s}}\), and solution time for instance 18/02/2015
Full size table

Table 5 shows the growing amount of variables and constraints for tighter linearizations leading to higher solution times. While the number of variables and constraints only slightly deviates from the average within each row, the solution time is in general higher for higher non-electric gas loads. Moreover, Case II turns out to be harder than Case I, which can be seen from the average solution time and the higher number of instances that hit the time limit.

Table 6 shows an instance, where solution times are higher than the averages from Table 5. The overall objective function (sum of \(C_{\mathrm{e}}\), \(C_{\mathrm{c}}\), and \(C_{\mathrm{s}}\)) is always higher for tighter approximations due to reduced feasible sets. This is mainly driven by the compressor costs \(C_{\mathrm{c}}\) for Case I and by the compressor costs \(C_{\mathrm{c}}\) and especially the shedding costs \(C_{\mathrm{s}}\) for Case II. This can be verified by the average power of the compressor in every dispatch period for all instances, as shown in Fig. 3 for Case I and II. The average power is in general higher for tighter approximations and moreover higher for higher non-electric gas loads, thus especially during the evening. Higher compressor power leads to higher compressor costs due to (17)–(23).

Fig. 3
figure 3

Average power of compressor for every dispatch period and accuracy level

Full size image

As Case II prohibits incoming flow from the south, all exits, e.g. the electric gas loads in the south, have to be served by the northern entries. This leads to higher compressor costs for Case II compared to Case I, which follows from the necessary compression for the gas from the north. This can be seen in Fig. 3 again, where the average compressor power is much higher for Case II than Case I. Even if the compressor in Case II is able to fulfill the requirements of the electricity problem in Sect. 2.1, it is not able to guarantee the demands of the non-electric gas loads any more, thus \(s_{u}^{} > 0\) in (4). Figure 4 visualizes the average of all these sheddings according to Case II. As all sheddings occur in the south, therefore a zoom is performed into this region in Fig. 4. Yellow circles represent non-electric gas loads without shedding, while magenta circles represent non-electric gas loads with shedding. The higher the averaged shedding, the bigger the radius of these circles. As the LNG entry is closed for Case II, the non-electric gas loads in the south will suffer according to the model.

Fig. 4
figure 4

Southern part of the network from Fig. 2 without the control valve

Full size image

4.2.2 Comparison with single electricity model

An additional way to depict the effect of the combined optimization of the two systems is by comparing it with the single electricity model, which is the solution of the electricity problem in Sect. 2.1 without any consideration of the coupling to the gas network. The solution of the single electricity model constitutes the optimal values of the energy production of all power plants, including CCGTs as electric gas loads. The differences of these values compared to the values of the combined optimization problem are investigated. Figure 5 shows the average difference of all values per hour for all instances and accuracy levels.

Fig. 5
figure 5

Average difference of produced power by each CCGT between single electricity model and MIP for every dispatch period

Full size image

It can be seen that the average production of the electric gas loads is reduced for many dispatch periods for both cases. The sharpest decline takes place during the evening peak hours and is even more evident for Case II.

The reduced energy production can be attributed to the fact that the gas network cannot allocate the desired amount of gas to the electric gas loads according to the optimal solution of the single electricity model. The difference is more evident in Case II when the supply capacity of the gas network is further reduced.

The lost production of the gas-fired power plants is compensated by an increase from the other, more expensive electricity production sources. This leads to higher electricity costs \(C_{\mathrm{e}}\), which can be seen in Table 7 for both cases.

Table 7 Average increase of \(C_{\mathrm{e}}\) for both cases
Full size table

5 Conclusion

In this paper, the coupled optimization of electricity and natural gas systems was investigated for the real-world test system of Greece. The extended incremental method, as well as the outer approximation method, were applied for the linearization of the nonlinear equations that govern the physics of gas transport and for the coupling constraint of the two energy systems. Four different discretization accuracy levels were selected and the experiments were conducted on real system data for two distinct test cases. The computed results showed the effectiveness of the proposed combined optimization method by depicting.

  1. 1)

    The increase of objective function costs and solution times with respect to tighter approximations for all instances examined.

  2. 2)

    The appearance of shedded non-electric gas loads during certain instances, caused by the need to provide the demand of the electric gas loads, in order to maintain the electricity system balance.

  3. 3)

    The effect of the gas system’s constraints on the electricity model, which was expressed in the increase of the electricity unit commitment cost.

This work is intended to act as a basis for further investigation into computations and comparisons of the coupled optimization problem, with additional consideration of stochasticity for the electricity system and the implementation of instationary gas flows for the gas system.