Optimization of a Pumped-Storage Fixed-Head Hydroplant: The Bang-Singular-Bang Solution

Abstract

We consider the problem of the optimization of the functioning of a pumped-storage hydroplant. The problem can be mathematically formulated as an optimal control problem, and when the considered hydromodel is of the fixed-head type, an added complication arises: the solution is of the bangsingular-bang type. In this paper, we propose a simple and efficient optimization algorithm to find the solution.

1. Introduction

This paper has a clear aim: optimizing the functioning of a pumped-storage hydroplant, with the context of the decentralized Spanish electricity market [1].

In the literature (see [2, 3]), numerous simplifications are performed to model plants of this kind, in order to obtain an extremely simple mathematical formulation which permits the use of commercial packages of, for instance, linear programming. Readers are referred to [4], which constitutes an excellent list of over 100 of relevant research papers on the hydroplant modeling.

However, the optimization of pumped-storage hydroplants is a complex Dynamic Programming problem [5], with given time and final state. The problem can be mathematically formulated as a standard Lagrange-type optimal control problem [6]. In a previous paper [7] the authors considered a variable-head model for the pumped-storage hydroplant. The solution can be obtained using Pontryagin’s Maximum Principle [5]. In another previous paper [8] the authors considered a fixed-head model for a conventional hydroplant, without pumping capacity. The solution is of the bang-bang type, and the absence of singular arcs was proven. However, when the considered model is fixed head and pumped storage simultaneously, an added complication arises: the solution is of the bang-singular-bang type.

Different methods for determining optimal controls with a possibly singular part have already been developed. A popular approach introduced by Jacobson et al. [9] has been used by a number of researchers including Edgar and Lapidus [10, 11] and more recently by Chen and Huang [12]. This method involves solving the singular/bang-bang optimal control problem as the limit of a series of nonsingular problems. It is important to establish the limitations of these perturbation-based methods for practical problems. In fact, the convergence criterion described in [9] requires that the perturbation parameter, 𝜀, be sufficiently small; however, numerical difficulties result when 𝜀 approaches a zero limit. The reader is referred to [12–14] for further details. Maurer et al. [15] presented a numerical scheme for computing optimal bang-bang controls. They assume that every component of the optimal control is bang-bang and that there are only a known finite number of switching times. Assuming that the optimal control structure is known, [16] formulates a new finite-dimensional optimization problem involving the initial states, the switching times, and the final time 𝑡_𝑓 as optimization variables. Another idea is to guess a possible control structure and start with some additional bang-bang and singular arcs. The optimization approach will then lead to a solution where the lengths of all redundant arcs will be optimized to zero. Alternatively, the direct monotone structural evolution method [17] additionally takes advantage of the minimum principle optimality conditions to obtain the correct control structure. For the study of sufficient conditions, the reader is referred, for example, to [18, 19].

Faced with the existence of diverse general methods for bang-singular-bang problems, in this paper we have developed a specific algorithm for our problem. Our method needs no prior knowledge of the number and location of the bang-singular-bang arcs neither does it handle any parameter (like, e.g., discretization or a penalization factor) that has an influence on convergence or the precision of the result.

In this paper we propose a simple and efficient optimization algorithm to find the bang-singular-bang solution to our hydroproblem. Though we treat only a specific hydraulic problem, it should be noted that our method may be applied to other problems with the same characteristics. The paper is organized in the following way. In Section 2 we set out the corresponding mathematical problem. Section 3 summarizes the optimization algorithm based upon the previous mathematical fundaments. Finally, the results obtained in an example are presented in Section 4, and the conclusions reached in this study are discussed in Section 5.

2. Mathematical Optimization of a Pumped-Storage Hydroplant

A basic physicallybased relationship between the active power generated by a hydroplant, 𝑃 (in MW), the rate of water discharge, 𝑧^′ (in m³/s), and the effective head, ℎ (in m), is given by𝑧𝑃=^′ℎ𝐺, (2.1) where 𝐺 is the efficiency (in m⁴/h·MW) (see [20, 21]). For a large capacity reservoir, it is practical to assume that the effective head is constant over the optimization interval. Here the fixed-head hydroplant model is defined and 𝑃 is represented by the linear equation: 𝑃𝑧^′(𝑡)=𝐴𝑧^′(𝑡), (2.2) where 𝐴 represents the efficiency and diverse parameters related to the geometry of the station (see [20] for further details). Pumped-storage is a well-known type of hydroplants used for load balancing. The method stores energy in the form of water, which is pumped from a lower to an upper reservoir. At times of low electrical demand, electric power is used to pump water into the upper reservoir. During periods of high electrical demand, water is released back into the lower reservoir through a turbine, thereby generating electricity.

Taking into account the conversion losses of the pumping process and evaporation losses from the exposed water surface, a maximum of 70% or 85% of the electrical energy used to pump the water into the elevated reservoir can be regained. Thus, we must introduce the efficiency, 𝜂, in the model. Despite these losses, when scheduling is optimized, the system increases revenue by selling more electricity during periods of peak demand, when electricity prices get higher. Hence, when pumped-storage plants are considered, the function 𝑃 is defined piecewise as 𝑃𝑧^′∶=𝐴⋅𝑧^′if𝑧^′≥0𝜂⋅𝐴⋅𝑧^′if𝑧^′<0. (2.3) If we assume that 𝑏 is the volume of water that must be discharged throughout the optimization interval [0,𝑇], the following boundary conditions will have to be fulfilled: 𝑧(0)=0,𝑧(𝑇)=𝑏. (2.4) Besides the previous statement, we consider 𝑧^′(𝑡) to be bounded by technical constraints 𝑞_min≤𝑧^′(𝑡)≤𝑞_max,∀𝑡∈[0,𝑇]. (2.5) No transmission losses will be considered in our study. From the perspective of a generation company and within the framework of the new electricity market, said losses are not relevant, as power generators currently receive payment for all the energy they generate in power plant bars.

In this section, we focus on the new short-term problem faced by a generation company, in a deregulated electricity market, when preparing its offers for the day-ahead market. Our model of the spot market explicitly represents the price of electricity as a known exogenous variable.

In our problem, the objective function is given by hydraulic profit over the optimization interval, [0,𝑇]. Profit is obtained by multiplying the hydraulic production of the pumped-storage hydroplant by the clearing price, 𝜋(𝑡), at each hour, 𝑡. Taking our objective functional 𝐹(𝑧) in continuous time form, the problem is max_𝑧𝐹(𝑧)=max_𝑧^𝑇₀𝐿𝑡,𝑧(𝑡),𝑧^′(𝑡)𝑑𝑡=max_𝑧^𝑇₀𝜋(𝑡)𝑃𝑧^′(𝑡)𝑑𝑡, (2.6) onΩ=𝑧∈𝐶¹[0,𝑇]∣𝑧(0)=0,𝑧(𝑇)=𝑏;𝑞_min≤𝑧^′(𝑡)≤𝑞_max,∀𝑡∈[0,𝑇]. (2.7) A standard Lagrange-type optimal control problem can be mathematically formulated as follows:max_()𝑢,𝑧^𝑇₀𝐿(𝑡,𝑧(𝑡),𝑢(𝑡))𝑑𝑡=max_()𝑢,𝑧^𝑇₀𝜋(𝑡)𝑃(𝑢)𝑧𝑑,_′(0)(𝑇)𝑢=𝑢;𝑧=0,𝑧=𝑏,_min≤𝑢(𝑡)≤𝑢_max. (2.8) We define the Hamiltonian in normal form: 𝐻(𝑡,𝑧,𝑢,𝜆)∶=𝐿(𝑡,𝑧,𝑢)+𝜆𝑢=𝜋(𝑡)𝑃(𝑢)+𝜆𝑢. (2.9) The resulting Hamiltonian, 𝐻, is linear in the control variable, 𝑢, and results in an optimal singular/bang-bang control policy. In general, the application of Pontryagin's Maximum Principle [5] is not well suited for computing singular control problems as it fails to yield a unique value for the control. It is well known [6] that when the Hamiltonian is linear in 𝑢, the optimality condition (maximize 𝐻  𝑢) leads to 𝑢^∗(𝑡)=𝑢_max()𝑢ifΦ𝑥,𝜆>0_sing()𝑢ifΦ𝑥,𝜆=0_minifΦ(𝑥,𝜆)<0 (2.10) and 𝑢^∗ is undetermined if Φ(𝑥,𝜆)≡𝐻_𝑢=0. The function Φ is called the switching function. If Φ(𝑥^∗(𝑡),𝜆(𝑡))=0 only at isolated time points, then the optimal control switches between its upper and lower bounds, which is said to be a bang-bang-type control (i.e., the problem is not singular). The times when the OC switches from 𝑢_max to 𝑢_min or vice versa are called switching times. If Φ(𝑥^∗(𝑡),𝜆(𝑡))=0 for every 𝑡 in some subinterval [𝑡^′,𝑡^′′] of [0,𝑇], then the original problem is called a singular control problem and the corresponding trajectory for [𝑡^′,𝑡^′′], a singular arc. The case when Φ vanishes over an interval is more troublesome, because the optimality condition is vacuous.

In our problem, however, an added complication arises: the Hamiltonian is defined piecewisely as𝐻()[(𝑡)](𝑡)𝑡,𝑢,𝜆∶=𝐴⋅𝜋+𝜆𝑢if𝑢≥0𝜂⋅𝐴⋅𝜋+𝜆𝑢if𝑢<0, (2.11) and the derivative of 𝐻 with respect to 𝑢  (𝐻_𝑢) presents discontinuity at 𝑢=0, which is the point at which a sudden change in 𝐻_𝑢 is produced, as it is the border between the power generation zone (positive values of 𝑢) and the pumping zone (negative values of 𝑢).

The classical gradient of 𝐻 at 𝑢 is defined only when 𝐻 is differentiable at 𝑢. However, when nondifferentiable objective functions arise in optimization problems, the generalized (or Clarke's) gradient (see [22, 23]) must be considered. Remember that the Clarke’s gradient, 𝜕𝑓, can be calculated as a convex hull of (almost) all converging sequences of the gradients𝜕𝑓(𝑥)=colim∇𝑓𝑥_𝑖∶𝑥_𝑖⟶𝑥,𝑥_𝑖∉𝑆,𝑥_𝑖∉Ω_𝑓. (2.12) With the aim of obtaining a numerical solution, we first attempt to determine the structure of the solution that is, the sequence of the bang-bang and the singular parts.

Bearing in mind that the function 𝐻_𝑢(𝑡,⋅,𝜆) is discontinuous in 𝑢=0, and that if 𝑢≠0, then 𝐻⁺_𝑢≡𝐻⁻_𝑢, we have that the partial Clarke’s gradient is 𝜕_𝑢(𝑡)𝐻=𝐻𝐴⋅𝜋+𝜆if𝑢>0_+𝑢,𝐻⁻_𝑢(𝑡)if𝑢=0𝜂⋅𝐴⋅𝜋+𝜆if𝑢<0 (2.13) with 𝐻⁺_𝑢=𝐴⋅𝜋(𝑡)+𝜆 and 𝐻⁻_𝑢=𝜂⋅𝐴⋅𝜋(𝑡)+𝜆. The switching function is Φ(𝑧,𝜆)≡𝜕_𝑢𝐻, and the optimality condition leads to 𝑢^∗(𝑡)=𝑢_max(𝑡)𝑢=0if𝐴⋅𝜋+𝜆>0_sing=0if0∈𝜕_𝑢(𝑡)(𝑡)𝑢𝐻=𝐴⋅𝜋+𝜆,𝜂⋅𝐴⋅𝜋+𝜆_min=0if𝜂⋅𝐴⋅𝜋(𝑡)+𝜆<0. (2.14) On the other hand, the costate equation of Pontryagin’s Maximum Principle allows us to obtain 𝜆^′=−𝐻_𝑧=0⟶𝜆=𝜆₀(const). (2.15) With the previous mathematical development, we can determine the optimal solution: the bang-singular-bang segments and the boundary on which the solution is situated 𝑢^∗(𝑡)=𝑢_max=0if𝐴⋅𝜋(𝑡)>−𝜆₀𝑢_sing=0if−𝜆₀∈(𝑡)(𝑡𝑢𝐴⋅𝜋,𝜂⋅𝐴⋅𝜋)_min=0if𝜂⋅𝐴⋅𝜋(𝑡)<−𝜆₀. (2.16)

3. Algorithm

Based on the above theoretical results, in this section we describe the optimization algorithm.

The algorithm presents a series of advantages. First of all, our method needs no prior knowledge of the number and location of the bang-singular-bang arcs. Moreover (as we will see in the next section), it shows a rapid convergence to the optimal solution, and it can be run in a relatively short time due to the simplicity of the operations to be performed in this method.

Figure 1 illustrates the switching law, which is in fact the basis for the proposed method. The algorithm that leads to the optimal solution comprises the following steps.(i)First, 𝜋(𝑡) must be interpolated to obtain a continuous function.(ii)Second, for a given 𝜆, we have to determine the switching times: 𝑡₁,𝑡₂,…. These instants are calculated solving 𝐴⋅𝜋(𝑡)=−𝜆;𝜂⋅𝐴⋅𝜋(𝑡)=−𝜆. (3.1)(iii) Third, the optimal value 𝜆₀ must be determined in order for 𝑧_𝜆(𝑇)=^𝑁𝑢_𝑖=1𝛿^𝑢_𝑖⋅𝑞_max+^𝑁𝑙_𝑖=1𝛿^𝑙_𝑖⋅𝑞_min=𝑏 (3.2) with 𝛿^𝑢_𝑖 and 𝛿^𝑙_𝑖 being the duration of the 𝑖th bang-bang segment in the upper and lower bound, respectively, 𝑁_𝑢 and 𝑁_𝑙 the number of such segments, and 𝑧_𝜆(𝑇) the final volume obtained for each 𝜆.(iv)To calculate an approximate value of 𝜆₀, we propose a classic iterative method (like, e.g., bisection or the secant method).

Figure 1 

Illustration of the switching law.

In the next section the proposed method is applied to an example.

4. Example

A program was written using the Mathematica package to apply the results obtained in this paper to an example of a hydraulic system made up of one fixed-head hydroplant. The hydraulic model is 𝑃(𝑡)=0.000126821𝑧^′(𝑡) (4.1) with 𝐴=0.000126821. We shall also consider the technical constraints: 𝑞_min=−283866 (m³/h), 𝑞_max=394258 (m³/h). When the efficiency is 𝜂=1.15, these constraints, respectively, correspond to 𝑃_min=−41.4,𝑃_max=50 (MW).

In this paper, we focus on the problem that a generation company faces when preparing its offers for the day-ahead market. Thus, the classic optimization interval of 𝑇=24 h. was considered. The clearing price, 𝜋(𝑡) (euros/h·MW), corresponding to one day was taken from the Spanish electricity market [1] (see Table 1). Note that in real electricity markets, the clearing price, 𝜋(𝑡), is only known at each hour (𝑡=1,2,…,24). In this paper the known values of 𝜋(𝑡) were linearly interpolated with good results.

Table 1 

The clearing price, 𝜋(𝑡).

The solution may be constructed in a simple way by taking into account the above algorithm. We shall perform two tests.

4.1. Influence of the Efficiency

First, we consider the restriction on the volume: 𝑏=2⋅10⁶ (m³) as a fixed datum and then we proceed to analyze the influence of the efficiency, 𝜂. In Table 2 we present the optimal profit when the efficiency falls within the 1.35 to 1.15 range. We show the amount of pumped water and the improvement in profit with respect to the case without pumping.

Table 2 

Optimal solution varying 𝜂.

Table 3 

Optimal solution varying 𝑏.

As can be appreciated, significant improvements in profit of between 12–14% are obtained for even very poor efficiencies (in the 65–75% range). For the normal efficiency of pumped plants of around 80%, the improvement in profit is quite substantial, being around 20%.

In Figure 2 we present the solution of the algorithm when the efficiency is 𝜂=1.20. The optimal value obtained for −𝜆₀ is 0.01139601030789622, and the switching times are 𝑡₁=1.2345,𝑡₂=7.85646,𝑡₃𝑡=8.46727,₄=14.52,𝑡₅=18.9881,𝑡₆=22.7759. (4.2) Figure 3 presents the optimal hydro-power, 𝑃, for three cases: (a) without pumping, (b) with 𝜂=1.2, and (c) with 𝜂=1.15.

Figure 2 

Value of −𝜆0.

(a)

(b)

(c)

Figure 3 

Optimal hydro-power, 𝑃(𝑡).

It can clearly be seen that when 𝜂=1.2 (standard efficiency of 80%), the plant pumps during the most off-peak period in the clearing price 𝜋(𝑡). However, when 𝜂=1.15, and hence the efficiency increases to 85%, the plant also pumps in the small off-peak period in the clearing price 𝜋(𝑡) that usually occurs around 17.00 h.

The singular arcs correspond to what, in a previous paper [24], we called a smooth transition. This is a qualitative aspect of the solution of the problem. The discontinuity of the derivative of the Lagrangian does not translate as discontinuity in the derivative of the solution. In fact, it is verified that the derivative of the extremal where the minimum is reached presents an interval of constancy.

The algorithm runs very quickly (see Figure 4). In case (b), with 𝜂=1.2, 7 iterations were needed, with the CPU time required by the program being 1.1 sec on a personal computer (Pentium IV/2 GHz). The secant method was used to calculate the approximate value of 𝜆 for which Error=𝑧_𝜆(𝑇)−𝑏<tol (4.3) with tol=50 (m³). The secant method has provided satisfactory results using these initial values: 𝜆^min=min𝐴𝜋(𝑡);𝜆^max=max𝜂𝐴𝜋(𝑡). (4.4)

Figure 4 

Convergence of the algorithm.

4.2. Influence of the Volume Available

Second, we consider the efficiency 𝜂=1.2 as a fixed datum, and then we proceed to analyze the influence of the restriction on the volume 𝑏. We compare with the solution without pumping.

As can be seen (see Table 1), the more the water available in the reservoir, the lower the profit obtained from pumping it. In fact, if the average flow of the river was such that volumes of water of 4⋅10⁶ (m³) were available, it would not be of interest to consider building a pump plant.

This type of study is very interesting when designing the most suitable type of plant for each hydrographic basin.

5. Conclusions

In this paper we have presented a simple and efficient optimization algorithm to optimize the functioning of a pumped-storage hydroplant. When the hydroplant is of the fixed-head type, the optimal solution presents a very particular form: the bang-singular-bang solution. Faced with the existence of diverse general methods for bang-singular-bang problems, in this paper we have developed a specific algorithm for our problem. Our algorithm allows the optimal solution to be obtained easily. The obtained results provide real-time information to determine which configuration is preferable in each specific real situation of the electricity market. Our theoretical results have been illustrated with numerical examples.

As to future perspectives, in spite of having presented a hydraulic example in this paper, it should be noted that our method may be applied to other problems with the same characteristics, with a constant value of the adjoint. Another interesting future line of study could be that of adapting the technique present in [19] to our hydraulic problem. In this paper the authors present an application to a free-flying robot, where it has been shown that the singular arc of the control can be interpreted as a bang-bang control of a so-called augmented smooth control problem.

Acknowledgment

This work was supported by the Spanish Government (MICINN, Project: MTM2010-15737).