1 Introduction

1.1 Motivation

During the last 20 years, natural gas (NG) has taken a leading role in the electric energy supply due to the increasing installation of natural gas fired power plants (NGFPPs) for electricity and steam generation. This increase has been driven by technical, economic and environmental benefits. Furthermore, in power systems with large shares of intermittent renewable generation, NGFPPs become one of the swing resources used to provide flexibility [1].

According to the International Energy Agency (IEA), power generation remains the most important driver of rising global gas demand, accounting for almost 40% of total increase in gas demand between 2014 and 2040 [2]. Therefore, as NGFPPs are one of the major natural gas consumers, there is a close interdependency between the natural gas and electric power systems.

From the viewpoint of available energy resources optimization within a region, the combined optimization of electric power and natural gas systems enables to achieve better results than in the case where these energy systems are considered separately [3, 4], and thus yields higher economic benefits to the region. Therefore, there is a need to develop new methodologies to treat the energy systems in an integrated manner, which leads to improve the decision making of systems operators and market players.

The integrated long and medium-term operational planning (LMTOP) of energy systems assesses the optimal management of energy resources and the electricity and natural gas prices forecasting. The optimal solution of the LMTOP depends strongly on the prices and availabilities of fuels (natural gas, coal, oil fuels, nuclear), and also on the availability of renewable energy resources (hydropower, solar, wind, biomass). Among others, fuel prices, energy loads, water flows, renewable energy sources are uncertainties associated with the planning procedure and needs to be properly considered in the LMTOP problem.

It should be noted that the integrated LMTOP is strongly affected the management of water for power production in large controllable hydro reservoirs. In the same way, natural gas storage facilities (e.g. depleted oil or natural gas fields, salt caverns, aquifers, liquefied natural gas reservoirs, line pack of pipelines) can be pointed out as infrastructures that have influences on the medium-term optimal decisions [5].

1.2 Literature review

Several approaches that address the integrated modeling and analysis of energy systems in comprehensive and general way have been presented. Also there are commercial models that consider electric power and natural gas systems coordinately [stochastic dual dynamic programming (SDDP)] [6, 7]. These approaches take into account multiple energy carriers; particularly electricity and NG systems interaction and the combined operation have been investigated. Nevertheless, these approaches do not consider properly the joint operation of the storage systems (water reservoirs, NG storage and line pack of pipelines). Rubio et al. [8] summarized the main approaches and models, which deal with the integrated operational planning of multiple energy carrier systems [9–13]. Ojeda-Esteybar et al. [14] presented an integrated scheduling of energy storages considering a simplified energy carrier network.

Recently, Salimi et al. [15] extended the energy hub concept of Geidl et al. [11], including more elements to improve the interaction between energy carriers. This approach embeds a simplified operational planning model to assess the elements expansion, considering physical constraints on NG and electricity networks. The evaluation is made on a 24 h period of time for typical days. Therefore, large energy reservoirs cannot be optimized.

Unsihuay-Vila et al. [16] proposed a long-term supply/inter-connection expansion planning model of integrated electricity and NG systems. The proposed model is formulated as a mixed-integer linear optimization problem which minimizes the investment and operation costs to determine the optimal location, technologies and installation times of new facilities for electric power and NG systems over a long range planning horizon. However, NG storage facilities are disregarded.

Wu et al. [17] presented a stochastic security-constrained unit commitment model for the optimization of coordinated midterm water and NG supplies. The stochastic model considers random outages of system components, load forecast errors and water inflow uncertainty. Also, Li et al. [18] proposed an integrated model for assessing the impact of interdependency of electricity and natural gas networks on power system security. The integrated model incorporates the NG network constraints into the optimal solution of security-constrained unit commitment.

Cole et al. [19] described and assessed the coordinated running of two long term scheduling models: the rice world gas trade model (RWGTM) of the NG sector and the regional energy deployment system (ReEDS) model of the U.S. electricity sector. The two models successfully converge to a solution under reference scenario conditions. This paper demonstrates that the integrated models produced better regional-level results as when running in a stand-alone form.

Nevertheless, the cited works do not tackle the problem arising from the modeling of energy reservoirs, their coordinated management, the modification in their operational planning and the economic impact over both energy systems involved.

1.3 Contribution and paper organization

This paper makes the following contributions.

  1. 1)

    This work proposes a medium term operational planning model that integrates the electric power and natural gas systems with their networks and energy storage facilities, such as hydroelectric power plants with water reservoirs and NG storage facilities.

  2. 2)

    The close interdependency between the scheduling of the stored energy resources in both energy systems is assessed and influence between each other.

  3. 3)

    Nodal electricity and NG prices are obtained from the shadow prices (marginal cost) of the nodal balance of each system. The impact on prices due to the integrated optimization of the storable energy resources is assessed.

  4. 4)

    A multi-stage decomposition technique is applied to real large-scale energy systems showing its performance to find robust solutions.

The remainder of this paper is organized as follows. The interactions between natural gas and electricity systems and the energy storage modelling are described in Sect. 2.1. Section 3 presents the modeling assumptions and mathematical formulation. Section 4 exemplifies the performance of the proposed methodology through the resolution of two study cases. Finally, Sect. 5 concludes the present work.

2 Combined management of energy resources

2.1 Interactions between natural gas and electric power systems

During the last three decades, the installation of NGFPPs has led to increasing interdependencies between electricity and natural gas sectors. From 1993 to 2013, the worldwide electricity production using NG as fuel has doubled, from around 10% to nearly 24% [2]. By 2014, the share of electricity production using NG was 51% in Argentina [20], 27% in USA [21], and 29% in the UK [22].

The interdependence between electric power and NG systems can be explained from the system operation viewpoint [23]. The NG availability for NGFPPs is constrained by the maximum capacity of gas production/importation, the limited transmission capacity of gas network, and the priority scheme for the NG supply in case of shortages, on which residential and commercial customers typically take precedence over large consumers and NGFPPs. On the other hand, the dispatch of NGFPPs determines the total amount of NG consumption and the gas flows through the pipelines. Contingencies in NG infrastructure may lead to a loss of multiple NGFPPs, jeopardizing the security of the electric power system.

These interactions can also be explained from an economic viewpoint [8]. The market arrangement with their regulatory frameworks over electric power and NG systems affects the level and the dynamics of their interdependence. Generation companies with NGFPPs take part simultaneously in both markets (producers in electricity market and consumers in gas market). Thus they are in a better position for price arbitrage between both markets. Liberalized and flexible market structures facilitate this practice, which is required to reach an electricity and gas partial economic equilibrium. According to electricity and NG market prices, and the marginal heat rate of their plants, these companies can decide to use gas and sell electricity in the electric market, or resell previously contracted gas on the gas market and purchase electricity to meet their commitments [4, 23].

There are two indexes that denote the interaction degree between electric power and NG systems. The first one is NG consumption for power generation as a share of the total NG demand; and the second one is the share of electrical energy produced by NGFPPs [23]. Both shares depend not only on the installed capacity of NGFPPs, but also on the availability of other energy resources (hydroelectricity, nuclear, renewables, etc.), the relative fuel prices (NG, coal, and oil derived products) and the flexibility of the NGFPPs for switching or mixing different fuels.

Under the light of all these describe conditions, there is a certain, strong and rising interdependency between NG and electricity sectors.

2.2 Energy storage modelling: temporal coupling

The main feature of an integrated energy system with storage units is to use the energy stored in the reservoirs associated to hydroelectric plants and/or NG storage facilities to supply the demand, thereby avoiding the use of expensive fuels in thermal power plants. However, the limited availability of the stored energy leads to a complex mathematical problem, since the energy storage facilities creates a coupling in the operational decisions between those made in the present (and their future consequences) and those decisions that should be made in future.

This temporal coupling requires solving the optimization problem for whole time horizon as a unique large scale optimization problem, increasing the number of variables and constraints, and therefore the computation time. This problem can be solved through the use of decompositions techniques which divide the overall problem into smaller sub-problems (Benders method [24] or Dantzig-Wolfe method [25]).

3 Mathematical formulation

The mathematical formulation presented in this section addresses the so-called deterministic problem, i.e., all input parameters such as, water flows, demands, and fuel prices, are single and known values. Among other methodologies, this problem formulation can be implemented within the Monte Carlo simulation method to deal with actual and inherent uncertainties associated with these parameters [26]. Therefore, a finite number of independent trials or simulations of the deterministic problem are solved in order to achieve a satisfactory level of confidence in the resulting probability distribution functions (PDF). Then, to evaluate the proposed model and to analyze the integrated optimization of the storage facilities, a deterministic problem is solved with the following assumptions.

  1. 1)

    The objective function and the nonlinear constraints are linearized in order to solve using linear programming optimization tools.

  2. 2)

    The time horizon (1–3 years) is divided into monthly or weekly stages. Each stage is divided into subperiods called blocks with different time duration (3–5 blocks).

  3. 3)

    The electricity and NG demands are represented by load duration curves approximated with a step function. Even though they could have similar profiles and hold a significant interdependence, it is assumed that they are independent.

  4. 4)

    The NGFPPs power productions are linearized using a constant net heat rate for each NGFPP. For all other thermal generation units, nonlinear production costs are linearized using a piecewise linear production curve, and linked to different fuel prices and their corresponding net heat rates (switch fuel feature).

  5. 5)

    The electric power flows are modeled through a DC flow model without losses, which takes into account the available transmission capacity imposed by the network and other operational constraints.

  6. 6)

    The NG flows are modeled through nodal NG balances and the nonlinear flow-pressure relation over the pipelines is linearized through a piecewise approximation function.

These assumptions allow achieving an appropriate trade-off between precision and computation time.

3.1 Objective function

From a centralized operation perspective for electricity and NG systems, the objective function (1) for an integrated planning of both energy carriers can be formulated as the present value minimization of electric power system operating costs and NG production costs within the planning time horizon, that is:

$${\text{OF:}}\,\hbox{min} \left[ {\sum\limits_{t} {\gamma^{t} \sum\limits_{k} {b^{k} } \left( {\sum\limits_{j} {C_{j} \cdot P_{j}^{kt} } + \sum\limits_{n} {CUE \cdot SE_{n}^{kt} } + \,\left. {\left. {\left. {\sum\limits_{g} {C_{g} \cdot W_{g}^{kt} } + \sum\limits_{n} {CUG \cdot SG_{n}^{kt} } } \right)_{k} } \right)_{t} } \right]} \right.} } \right.$$
(1)

where γ t is a discount factor of stage t obtained from an annual discount rate and b k is the number of hours of block k. The first two terms inside the parenthesis represent the electric power generation costs (C j ) of thermal plants j for block k of stage t (P kt j ) and costs of unsupplied electricity (CUE) for electricity shortage at bus n for block k of stage t (SE kt n ) respectively. The second two terms, the NG production costs (C g ) of the supplier production g for block k of stage t (W kt g ) and costs of unsupplied gas (CUG) for gas shortage at bus n for block k of stage t (SG kt n ) respectively. It should be noted that the NGFPPs (variable G kt s ) are not included in the objective function because their costs are counted in the NG system, through the production costs of NG system. This is because NGFPPs play as loads in the gas system. If NGFPPs were included in the objective function with their associated costs, the fuel costs of these units would be double counted.

3.2 Operational constraints

The objective function (1) is subject to a set of electric power system restrictions for all k blocks at each t stages:

$$- F_{m}^{\hbox{max}} \le \sum\limits_{m} \left[ M_{mn} \,\left(\sum\limits_{{i \in {\mathcal{I}}_{n} }} {( {H_{i}^{kt} } )} + \sum\limits_{{j \in {\mathcal{J}}_{n} }} {( {P_{j}^{kt} } )} + \sum\limits_{{s \in {\mathcal{S}}_{n} }} {( {G_{s}^{kt} } )} + SE_{n}^{kt} - LE_{n}^{kt}\right)\right] \le F_{m}^{\hbox{max} } \quad \forall m$$
(2)
$$\sum\limits_{{i \in {\mathcal{I}}}}^{{}} {\left( {H_{i}^{kt} } \right) + } \sum\limits_{{j \in {\mathcal{J}}}}^{{}} {\left( {P_{j}^{kt} } \right) + } \sum\limits_{{s \in {\mathcal{S}}}}^{{}} {\left( {G_{s}^{kt} } \right) + } \sum\limits_{{n \in {\mathcal{N}}}}^{{}} {\left( {SE_{n}^{kt} } \right)} = \sum\limits_{{n \in {\mathcal{N}}}}^{{}} {LE_{n}^{kt} }$$
(3)
$$H_{i}^{kt} = \rho_{i} \, \cdot Q_{i}^{kt} \,\,\,\,\,\,\,\,\,\,\,\,\forall i\,$$
(4)
$$0 \le Q_{i}^{kt} \le Q_{i}^{\hbox{max} } \,\,\,\,\,\,\,\forall i\,$$
(5)
$$\left\{ {\begin{array}{*{20}l} {0 \le P_{j}^{kt} \le P_{j}^{\hbox{max} } } \hfill & {\forall j} \hfill \\ {0 \le G_{s}^{kt} \le G_{s}^{\hbox{max} } } \hfill & {\forall s} \hfill \\ \end{array} } \right.$$
(6)

Constraints (2) and (3) correspond to the dc power flow model [27] for electrical network representation, where H kt i is the production of hydro plants i for block k of stage t; P kt j is the generation of thermal plants; G kt s is the NGFPPs production; F max m is the maximum transmission limit for line m; LE kt n is the electricity demand at bus n, for block k of stage t; M mn is the power transfer distribution factor in the line m to an increase in the injection at bus n; \({\mathcal{I}}_{n}\) is the set of hydro plants i connected to bus n; \({\mathcal{J}}_{n}\) is the set of thermal plants j connected to bus n; \({\mathcal{S}}_{n}\) is the set of NGFPPs s connected to bus n. The Lagrange multiplier (shadow price) of constraint (3) represents the electricity marginal cost of the slack bus.

Constraint (4) characterizes the hydropower production, where Q kt i represents the water flow rate of the hydropower plant i for block k of stage t and ρ i the production ratio of the hydropower plant. This parameter is function of the hydraulic head and the turbine-generator performance. To simplify the problem, ρ i is modeled as a constant average production ratio. Constraint (5) represents the maximum water flow limits of hydropower plants (Q max i ) and (6) are the maximum limits of thermal power plants (P max j ) and NGFPPs (G max s ).

The natural gas network is modelled through a nodal gas balance and the pipeline flows. The steady-state isothermal gas flow model [28, 29] is shown in (7):

$${\text{sign}}\left( {QG_{m} } \right) \cdot QG_{m}^{2} = K_{m} \cdot {\kern 1pt} \left( {S_{n}^{2} - S_{z}^{2} } \right)$$
(7)
$${\text{where}}\,\,\,\,\,\,\,\,\,\,\,\,{\text{sign}}\left( {QG_{m} } \right) = \left\{ \begin{aligned} + 1\quad {\text{if}}\, S_{n} \ge S_{z} \hfill \\ - 1\quad {\text{if}}\, S_{n} < S_{z} \hfill \\ \end{aligned} \right.$$

where QG m is the gas flow through a pipeline m; K m is the gas properties and pipeline m characteristics; S n is the upstream pressure; and S z is the downstream pressure.

These constraints are clearly nonlinear; nevertheless, as in [30, 31] the pipeline flows can be approximated as piecewise linear constraints, replacing the square pressures with auxiliary variables (π n , π z ) and modeling the nonlinear flow with a piecewise linear mixed integer problem formulation, as it shown in Fig. 1.

Fig. 1
figure 1

Piecewise linear pipeline flow approximations

Full size image
  1. 1)

    For passive pipelines:

    $$\sum\limits_{y = 1}^{Y} {\left( {Z_{my} \cdot QG_{my} + a_{my} \cdot \beta_{my} } \right)} = K_{m} \cdot \left( {\pi_{n} - \pi_{z} } \right)$$
    (8)
  2. 2)

    For active pipelines (compressor stations):

    $$\sum\limits_{y = 1}^{Y} {\left( {Z_{my} \cdot QG_{my} + a_{my} \cdot \beta_{my} } \right)} \le K_{m} \cdot \left( {\pi_{n} - \pi_{z} } \right)$$
    (9)

    where QG my is the gas flow; Z my is the slope; β my is the y-intercept of the y th piece for pipeline m. The binary variable a my is used to enable one segment of the piecewise linear function.

The bounds for the piecewise linear segments are:

$$a_{my} \cdot F_{my}^{\hbox{min} } \le QG_{my} \le a_{my} \cdot F_{my}^{\hbox{max} }$$
(10)

where F min my , F max my are the flow limits for the y th piece for pipeline m. The total pipeline flow is defined as the sum of the piecewise linear segments:

$$QG_{m} = \sum\limits_{y = 1}^{Y} {QG_{my} }$$
(11)

And because there can be only one active segment, we add an additional constraint:

$$\sum\limits_{y = 1}^{Y} {a_{my} } = 1$$
(12)

Also, the nodal balance in the NG network for all k blocks at each t stages is shown in (13):

$$\left( {\sum\limits_{{g \in \mathcal{G}_{n} }}{W_{g}^{kt} } + SG_{n}^{kt} - \sum\limits_{{s \in S_{n} }} {HR_{s}\cdot G_{s}^{kt} } + } \sum\limits_{{m \in \mathcal{M}_{n}^{+ } }} {QG_{m}^{kt} } - {\sum\limits_{{m \in\mathcal{M}_{n}^{ - } }} {QG_{m}^{kt} } - \sum\limits_{{p \in\mathcal{P}_{n} }}^{{}} {QIO_{p}^{t} } } \right) =LG_{n}^{kt} \quad \forall n$$
(13)
$$W_{g}^{ \hbox{min} } \le W_{g}^{kt} \le W_{g}^{ \hbox{max} } \,\,\,\,\,\,\,\,\,\,\,\forall g$$
(14)
$$- QG_{m}^{ \hbox{max} } \le QG_{m}^{kt} \le QG_{m}^{ \hbox{max} } \,\,\,\,\,\,\,\,\,\,\,\forall m$$
(15)
$$- QO_{p}^{ \hbox{max} } \le QIO_{p}^{t} \le QI_{p}^{ \hbox{max} } \,\,\,\,\,\,\,\,\,\,\forall p$$
(16)

where W kt g is the gas injection of the supplier g for block k of stage t; SG kt n is the gas shortage at bus n; LG kt n is gas demands at bus n; QG kt m is the pipeline flows; QIO t p is the NG storage p inflows/outflows. The gas consumption of NGFPP s is represented by the product of its production G kt s and the net heat rate HR s . In (13), HR s is modeled as a constant value. \(\mathcal{G}_{n}\) is the set of gas suppliers connected to bus n, \(\mathcal{M}_{n}^{ + }\) is the set of pipelines which flows are incoming to bus n, \(\mathcal{M}_{n}^{ - }\) is the set of pipelines which flows are outgoing to bus n and \(\mathcal{P}_{n}\) is the set of NG reservoirs connected to bus n.

The Lagrange multiplier (shadow price) of this constraint is the NG marginal cost of node n and it is an endogenous result of the optimization process. This cost represents the production cost of the NGFPPs connected to this node.

Constraint (14) represents the injection limits of NG suppliers (W min g , W max g ), and (15) the flow limits (QG max m ) of the pipeline m. Constraint (16) represents the storage inlet and outlet flow limits (QI max p , QO max p ) of the NG reservoir. These storage flow limits can vary with the NG volume stored, but in order to simplify the problem, we consider them as constants.

3.3 Temporal coupling constraints

The operational conditions outlined above are even more restricted by the constraints that link the state variables at different stages: the water balance equation of reservoirs (17) and the NG balance equation of storage facilities (18) for each stage t:

$$V_{e}^{t + 1} - V_{e}^{t} = \, - \mathop{\sum_{k}}\limits_{i \in \mathcal{I}_{e}} {b^{k} \cdot CF \cdot Q_{i}^{kt} \,} - SO_{e}^{t} + A_{e}^{t} + \left( \mathop{\sum_{k}}\limits_{i \in \mathcal{I}_{u}} {b^{k} \cdot CF \cdot Q_{i}^{kt} } + \sum\limits_{{u \in \mathcal{E}_{u} }} {SO_{u}^{t} } \right)\quad \forall e$$
(17)
$$VG_{p}^{t + 1} - VG_{p}^{t} = QIO_{p}^{t} \, \cdot \sum\limits_{k} {b^{k} } \quad \forall p$$
(18)

where V t+1 e  − V t e is the evolution of the stored volume of the water reservoir e within the stage t; Q kt i is the outflow of the hydroelectric plant I; CF is a unit conversion factor; SO e is a slack variable that represents the spilled outflow of the reservoir e and A e is the water inflow to the reservoir e. The last two terms in parenthesis represent the outflows of the upstream hydroelectric plants that discharge into the reservoir e, and the spillage of the upstream reservoirs that discharge into the reservoir e. \(\mathcal{I}_{e}\) is the set of hydroelectric power plants i associated with the water reservoir e; \(\mathcal{I}_{u}\) is the set of upstream hydroelectric plants i that discharge into reservoir e and \(\mathcal{E}_{u}\) is the set of upstream reservoirs that spill into reservoir e.

NG reservoirs can be distinguished as large-capacity reservoirs (as shown in Fig. 2) and small-capacity reservoirs [5]. The first ones are used as long term storage (seasonal cycle capability), so the inflows/outflows are reasonably constant in the short term. The second ones, (LNG storages, salt caverns, line pack of a pipeline) allow the operation of the equipment in the daily or weekly gas balance (weekly/monthly cycle capability).

Fig. 2
figure 2

Large-capacity NG reservoir modelling

Full size image

Thus, constraint (18) models large-capacity reservoirs with a single variable injection/withdrawal rate (QIO t p ) per stage t, where VG t p represents the stored volume of the natural gas storage p at the beginning of stage t. Constraint (19) models small-capacity reservoirs with a variable injection/withdrawal rate (QIO kt p ) for each block k of the stage t:

$$VG_{p}^{t + 1} - VG_{p}^{t} = \sum\limits_{k} {QIO_{p}^{kt} \cdot b_{k} } \,\,\,\,\,\,\,\,\,\,\,\,\,\,\forall p$$
(19)

Finally, constraint (20) represents the storage limits of the water reservoir and (21) the storage limits of the NG reservoir.

$$V_{e}^{ \hbox{min} } \le V_{e}^{t} \le V_{e}^{ \hbox{max} } \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\forall e$$
(20)
$$VG_{p}^{ \hbox{min} } \le VG_{p}^{t} \le VG_{p}^{ \hbox{max} } \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\forall p$$
(21)

where V min e , V max e are the minimum and maximum volumes of the water reservoir e. The base gas VG min p of the natural gas reservoir is the gas volume intended as permanent inventory to maintain adequate pressure and deliverability rate throughout the withdrawal period [5]. The gap between the maximum capacity of the reservoir (VG min p ) and the base gas is called working gas capacity of the storage facility.

3.4 Multi-stage Benders’ decomposition applied to large-scale optimization problems

When commercial MILP solvers are applied to large-scale non-convex optimization problems global optimal solution can hardly be achieved depending on the initial variables values. Moreover, the lack of robustness of MILP solvers can be checked when different initial variable values can result in different solutions. To cope with these problems, Benders decomposition can be implemented as a solution method for solving large-scale linear optimization problems [32, 33]. Multi-stage optimization problems have the following typical structure:

$$\hbox{min} \,\,\,\,\,\,\sum\limits_{t = 1}^{T} {\varvec{C}_{t} \cdot \varvec{X}_{t} }$$
(22)
$$\text{s.t.} \,\,\,\,\,\,\,\,\,\,\,\varvec{A}_{t} \cdot \varvec{X}_{t} + \varvec{E}_{t} \cdot \varvec{X}_{t - 1} \le \varvec{B}_{t}$$
(23)

These types of problems may be decomposed in T sub-problems, and solved using an iterative process. The dual dynamic programming (DDP) [33] algorithm consists of a forward and a backward pass. During the forward pass, an upper bound is obtained so that the convergence of the process can be tested. During the backward pass the Benders cuts are constructed in order to enhance the approximation of the future cost functions t + 1. The Benders cut, given as:

$$\alpha_{t + 1} + \left( {\varvec{\varPi}_{t + 1}^{{\left( {it} \right)}} \cdot \varvec{E}_{t + 1} } \right) \cdot \varvec{X}_{t} \ge \delta_{t}^{{\left( {it} \right)}}$$
(24)

where α is the future cost function; Π t is the dual variable vector of constraint (23) when the state vector is fixed at its optimal solution value at the previous stage, i.e., \(X_{t - 1} = X_{t - 1}^{*}\), it is the algorithm’s iterations counter, and δ t is a scalar term obtained as the product of B t and Π t for each iteration. The DDP algorithm and its stochastic variant, has undergone extensive research in recent years. A detailed explanation of both algorithms can be found in [32–35].

The sequential structure is portrayed in Fig. 3, in which the forward phase is solved using MILP for each stage, and the backward phase, used to generate the Benders cuts, is solved relaxing the constraints with integer variables (8)–(12) and using LP.

Fig. 3
figure 3

Operational layout of DDP

Full size image

4 Test systems

With the aim of evaluating the mathematical formulation and the proposed methodology, we present two study examples: a small-scale test system and a large-scale real system.

4.1 Case 1: small-scale system application

In order to assess the economic and operational impact of the combined optimization of energy storages, a small-scale test system has been presented. This test system allows a detailed result analysis of the integrated planning between the energy systems and their storage facilities. This test system is similar to the example used in [14].

The proposed test system consists of an electrical system and a NG system with three buses (B1, B2, and B3). The time horizon adopted for the LMTOP simulation takes 2 years divided in 24 monthly stages, with demands at each stage divided in 4 blocks of 40, 300, 270 and 120 h respectively. These systems are interconnected and interdependent in their operation.

The electrical system, as shown in Fig. 4, have a NGFPP on each bus (G 1, G 2, and G 3), one hydropower plant (H 3) connected to bus 3, and two thermal plants (P 1 and P 2), connected to bus 1. The production costs are associated with the cost of the fuels used and the specific fuel consumption rate of each unit. Additionally, it is modeled a fictitious unit (SE 1, SE 2 and SE 3) associated to the shortage cost, to represent the power generation shortage on each bus. The demands LE 1, LE 2, and LE 3 are placed at each node of the electrical system. The hydropower plant has a reservoir V 3. The stored volume at the beginning of the study period must equal the final volume, so only water intakes during the whole time horizon are optimized. The electrical network consists of three lines (F 12, F 13, and F 32) with their corresponding capacity limits.

Fig. 4
figure 4

Electric power system

Full size image

The NG system (as shown in Fig. 5) includes two production fields (W 21 and W 22) at bus B2, with production costs which consider the NG extraction cost. Artificial NG producers (SG 1, SG 2, and SG 3) connected to each node are used to simulate gas shortages. The gas demands LG 1, LG 2, and LG 3 are placed at each node. These demands don’t include the gas for electricity production. The example also simulates a NG reservoir VG 1 that can steadily withdraw from or deliver to the NG network on each monthly period. The NG network has two pipelines (QG 41 and QG 53) and two compressor stations (QG 24 and QG 23). Active and passive pipelines are modeled with a 2 step piecewise linear approximation.

Fig. 5
figure 5

NG system

Full size image

The connection between electrical and NG system occurs through the NGFPPs G 1, G 2, G 3 located at each node.

4.1.1 Study cases solved

The optimal LMTOP for the test system is solved deterministically (one scenario) taking into account the objective function (1) and constraints (2)–(21). Table 1 shows the optimal use of storable resources and their impact on the energy systems involved. For a clearer understanding of each storage resource contribution to operating cost reduction, the following study cases are proposed.

Table 1 Storage facilities of solved cases
Full size table

For Case A and B, the hydroelectric plant is modeled as a run-of-river power plant, that is, it has no storage capacity for periods longer than one stage, but it can optimize water intake within that stage. Note that in these cases, water flows do not produce spillage, so the amount of available water (available energy) to store is the same. The above explanation is essential, because if in any of the proposed cases the water is spilled, the total hydro energy would not be equal, therefore the cases cannot be compared. Case C and D, the hydroelectric plant is modeled with the reservoir V 3.

For Case B and D, the NG system includes the NG storage VG 1. Case A and C, NG storage is not considered.

The flow restriction over electric lines and pipelines could affect the scheduling of energy storages, depending on where the reservoirs are located. In order to compare only the effect of energy storage (water or NG) over the total operating costs, electric lines and pipelines capacity limits have been relaxed.

The mathematical problem is composed of 3672 variables and 1828 constraints. The problem has been programmed using MATLAB software and solved with AMPL tool. The case studies data can be found in Appendix A.

4.1.2 Problem results

Table 2 shows the economic results and energy shortages for both energy systems. The results illustrate that the modeling of all energy storage systems (water reservoir, NG storage) contribute in the long term to avoid energy shortages, thus, contributing to economic savings.

Table 2 Optimization results
Full size table

Furthermore, the integrated and coordinated use of energy storage facilities (Case D) contribute to the optimal use of scarce resources, quantitatively measured by the total operating cost of the electric power system, NG production and energy shortages, which are the lesser of all the cases.

The influence of energy storages in electricity and NG systems can be noted by comparing the results obtained with the proposed cases (with and without energy storage).

Figure 6a, b shows the storage evolution of the hydroelectric reservoir V 3 (expressed in hm3 i.e. 106 m3), streamflow A 3, hydro generation (HG) and electric demand (ED), of Case D (with hydro reservoir) and Case B (without hydro reservoir).

Fig. 6
figure 6

Evolution of V 3, A 3, HG and ED

Full size image

In Case D (Fig. 6a), the reservoir keeps and saves water when power requirement is low (stages 1–4), and the hydropower plant generates at maximum power in periods when demand requirements and thermal costs are higher (stages 6, 7, 11 and 12).

However, in Case B (Fig. 6b), the hydropower plant can only generate the available water within the stage, so its generation profile equals the streamflow profile. This leads to a suboptimal use of the water resource, because the hydro plant does not have the primary resource during times of high demand, so it is necessary to dispatch high cost generation or even to have energy shortages

Figure 7a, b shows the evolution of NG storage VG 1, (expressed in dam3 i.e. 103 m3) production of NG wells W 21 and W 22, and total NG demand (including demand of NGFPP’s) of Case D (with NG storage) and Case C (without NG storage).

Fig. 7
figure 7

Evolution of NG storage VG 1, production of NG wells W 21 and W 22

Full size image

As in the electrical system, in Case D (Fig. 7a) the gas storage VG 1 saves gas in lower demand periods, and then injects into the system when demand is higher (including the consumption of NGFPPs).

In contrast, in Case C (Fig. 7b) NG wells W 21 and W 22 supply the gas required by the demand and the remaining injection capacity is delivered to NGFPP for power generation. This implies that the NG wells are responsible for flow regulation in the NG system.

Comparing both figures, the NG injection peak achieved in Case D is greater than in Case C, due to the additional injection capacity of the NG storage in the peak demand period (stages 6, 7).

The economic influence of the reservoirs modeled on the electrical and NG systems can be seen through the comparison of nodal marginal costs of electricity and NG on each evaluated Case.

According to economic theory applied to electricity markets, the marginal cost of electricity is the variable cost of the last generation unit dispatched, called the marginal unit. In the same way, the marginal cost of NG is the variable cost of production of the NG marginal injection (gas wells or NG storage).

Figure 8 shows the weighted average marginal costs of electricity (weighed for the time duration of the block) corresponding to the bus B1 of the four cases modeled. Also, Fig. 9 shows the weighted average marginal costs of natural gas corresponding to the bus B1 of the four cases.

Fig. 8
figure 8

Weighted average marginal costs of electricity

Full size image
Fig. 9
figure 9

Weighted average marginal costs of NG

Full size image

The lack of storage facilities (water reservoir, NG storage) in the energy system produce large dispersion in the marginal costs of electricity and NG, increasing substantially at peak demand periods, as it is shown in Case A and to a lesser extent in cases B and C (Figs. 8, 9).

The optimal management of water resources produces a flattening of the nodal marginal costs of electricity (Fig. 8), generating less dispersion with a corresponding economic benefit to the sector. This statement can be seen on the comparison of the marginal costs of electricity in Cases C and A (with and without water reservoir).

It may also be noted the influence of the NG storages over the electric power system. When comparing marginal costs of electricity in Cases A and B (Fig. 8), it can be observed that the existence of NG storages (Case B) produces a decrease in the peaks of the electricity marginal cost (compared to Case A).

Similarly, in Case B the modeling of NG storage flattens the marginal costs of NG (compared to Case A), reducing the scattering and achieving an economic benefit to the sector (Fig. 9). This is because, in stages of higher NG demand, the NG storage runs as a gas producer, competing with the NG wells.

The economic influence of the water reservoir over the NG system can be noted by comparing the marginal costs of NG in cases A and C (Fig. 9). It can be seen that the existence of the water reservoir (case C) produces a flattening of the marginal costs of NG.

The integrated planning of energy storages (Case D) leads to the smallest dispersion of marginal costs of electricity and NG (Figs. 8, 9). It should be pointed that in some stages of Cases A, B and C, marginal costs are lower than those obtained in Case D. This does not mean a reduction of total costs, because in periods of peak demand, marginal costs in Cases A, B and C increase significantly over Case D.

4.2 Case 2: large-scale system application

In order to assess the proposed methodology over large-scale systems, we have applied it over the Argentinean energy system, including hydrothermal power and NG systems with their associated networks.

4.2.1 Argentina’s energy system

The Argentinean energy system is mainly composed of a hydrothermal electric power system and a NG system. Both systems have extensive networks that allow the interconnection from the main production centers to consumption areas. It should be noted that, due to the considerable share of NGFPPs in the electric generation set, there is a strong interrelationship between the systems involved.

The Argentine Power Interconnection System (APIS or SADI for its acronym in spanish) is composed by a diverse set of power plants, interconnected by a geographically extensive and weakly meshed network of 500 and 220 kV, to meet demand (131 TWh for year 2014), largely located in Buenos Aires (50%). The APIS have international links with Uruguay, Brazil and Paraguay. In order to produce the electricity needed to meet demand, several types of primary resources are used: hydroelectric (31%), NG (48%), fuel oil (8%), diesel oil (7%), uranium (4%), coal (1%), renewables (1%) (wind, solar, hydro, biomass) [36].

The Argentinean Natural Gas System (ANGS) has three main productive NG basins: Northwest, Neuquén and Austral. Also, there are two pipelines coming from Bolivia that allow them to import up to 17 MMCM/days (Mega standard cubic meters per day (106 m3/days), commonly used in NG sector). The NG supply has been increased from 2008 with the injection of regasified LNG (liquefied NG). For this purpose, two maritime ports (southern and northern Buenos Aires) were conditioned to allow the connection of LNG regasification vessels to inject gas into the NG network. For this, floating storage and regasification units (FSRU) with 14 MMCM/days of maximum injection capacity were rented. For the LNG delivery from LNG carriers to the FSRU without interrupting the NG supply into the network, a ship-to-ship LNG transfer configuration was implemented. The NG transportation system consists of a set of pipelines that start in the NG basins and converge in Buenos Aires (city gate). The NG transportation system has a total of 52 compressor stations. In all stations, NG-fired turbines burn a small portion of NG from the pipeline to generate the energy needed to run the compressors.

4.2.2 Electric power system modeling

The APIS is organized as a set of areas with similar power supply and demand characteristics (Fig. 10), linked with a simplified transmission system that simulates the 500 kV network. Figure 10 also shows the installed capacity of hydroelectric power (HP), thermoelectric power (TP) and nuclear power (NP), and the transmission capacity of power lines. We manage 25 U of NGFPPs, 47 U of NGFPPs multi-fuel (could use NG or liquid fuel), 10 U of thermoelectric power plants that use liquid fuels, 2 U of nuclear power plants, 9 U of hydroelectric plants with water reservoirs and 5 U of small hydroelectric plants. The monthly electric demand is modeled as a load duration curve. In order to simulate the power flow through the transmission network, a DC power flow without losses was applied.

Fig. 10
figure 10

APIS modeling scheme

Full size image

4.2.3 Natural gas system modeling

The ANGS, as is in APIS, is organized as a set of areas with similar gas supply and demand characteristics (Fig. 11), linked with a simplified transport system that simulates the natural gas pipeline network. In Fig. 11 also displays the production capacity of natural gas basins (expressed in MMCM/days), FSRUs, and transmission capacity of pipelines. The monthly natural gas demand is modeled as a load duration curve. In order to simulate the natural gas flow through the pipeline network, active and passive pipelines are modeled with a 2 step piecewise linear approximation.

Fig. 11
figure 11

ANGS modeling scheme

Full size image

4.2.4 Problem results

The operational planning of energy systems is solved deterministically (one scenario) for a study time horizon of 1 year, divided into 12 monthly stages. Each stage is divided into 4 blocks of different time.

The electricity and NG demands are represented by load duration curves approximated with a step function of 40, 300, 270, 120 h. It is assumed that there is independence between the two types of demand, which are similar in appearance, even more so the relationship between them is significant. Information of the electric power system was obtained from the Argentinean ISO CAMMESA [37] and the natural gas system from the Argentinean Regulator ENARGAS [38, 39].

In order to solve the problem, we use the proposed model. The optimization problem includes 15072 variables and 9924 constraints. With the decoupling procedure, each subproblem has 1256 variables and 827 constraints. The model has been implemented using a computer with i7 processor with 16 GB RAM, and programmed in MATLAB software and with MILP and LP AMPL software. The problem was solved in 21 iterations in 9.2 min. Table 3 shows main results of the problem.

Table 3 Optimization results (Argentinean case)
Full size table

5 Conclusion

The model and optimization approach presented in this work allows the assessment of the interactions between the energy storage facilities and their economic impact over their optimal systems scheduling. The simulation results of the study examples have shown the benefits of an integrated operational planning of electric power and NG systems and the close interdependency between the energy resources stored in both energy systems.

As a result of the optimization procedure, nodal electricity and NG prices can be obtained from the shadow prices (marginal cost) of the nodal balance of each system. Thus, the economic impact on prices due to the integrated optimization of the storable energy resources can be assessed.

The proposed methodology can be applied to real large-scale energy systems, like the Argentinean energy system, because the problem is conveniently decomposed according to the problem stages (typically months or weeks).

The proposed model can easily be extended to systems that include other renewable energy resources used for electric power generation. The model can also be extended to properly consider the randomness of the parameters uncertainties (e.g. water flows, energy demands, fuel costs, renewable productions).