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

Authors: L. Bayón; J. M. Grau; M. M. Ruiz; P. M. Suárez
Journal: Mathematical Problems in Engineering
 (2011)

---

## 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.


---

## Body

## 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 tf 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 Section2 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,P (in MW), the rate of water discharge, z′ (in m3/s), and the effective head, h (in m), is given by(2.1)P=z′hG,
where G is the efficiency (in m4/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 P is represented by the linear equation: (2.2)P(z′(t))=Az′(t),
where A 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 of70% 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 P is defined piecewise as (2.3)P(z′):={A⋅z′ifz′≥0η⋅A⋅z′ifz′<0.
If we assume that b is the volume of water that must be discharged throughout the optimization interval [0,T], the following boundary conditions will have to be fulfilled: (2.4)z(0)=0,z(T)=b.
Besides the previous statement, we consider z′(t) to be bounded by technical constraints (2.5)qmin⁡≤z′(t)≤qmax⁡,∀t∈[0,T].
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,T]. Profit is obtained by multiplying the hydraulic production of the pumped-storage hydroplant by the clearing price, π(t), at each hour, t. Taking our objective functional F(z) in continuous time form, the problem is (2.6)max⁡zF(z)=max⁡z∫0TL(t,z(t),z′(t))dt=max⁡z∫0Tπ(t)P(z′(t))dt,
on(2.7)Ω={z∈Ĉ1[0,T]∣z(0)=0,z(T)=b;qmin⁡≤z′(t)≤qmax⁡,∀t∈[0,T]}.
A standard Lagrange-type optimal control problem can be mathematically formulated as follows:(2.8)max⁡(u,z)∫0TL(t,z(t),u(t))dt=max⁡(u,z)∫0Tπ(t)P(u)d,z′=u;z(0)=0,z(T)=b,umin⁡≤u(t)≤umax⁡.
We define the Hamiltonian in normal form: (2.9)H(t,z,u,λ):=L(t,z,u)+λu=π(t)P(u)+λu.
The resulting Hamiltonian, H, is linear in the control variable, u, 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 u, the optimality condition (maximize H  u) leads to (2.10)u*(t)={umax⁡ifΦ(x,λ)>0usingifΦ(x,λ)=0umin⁡ifΦ(x,λ)<0
and u* is undetermined if Φ(x,λ)≡Hu=0. The function Φ is called the switching function. If Φ(x*(t),λ(t))=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 umax to umin or vice versa are called switching times. If Φ(x*(t),λ(t))=0 for every t in some subinterval [t′,t′′] of [0,T], then the original problem is called a singular control problem and the corresponding trajectory for [t′,t′′], 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(2.11)H(t,u,λ):={[A⋅π(t)+λ]uifu≥0[η⋅A⋅π(t)+λ]uifu<0,
and the derivative of H with respect to u  (Hu) presents discontinuity at u=0, which is the point at which a sudden change in Hu is produced, as it is the border between the power generation zone (positive values of u) and the pumping zone (negative values of u).The classical gradient ofH at u is defined only when H is differentiable at u. 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, ∂f, can be calculated as a convex hull of (almost) all converging sequences of the gradients(2.12)∂f(x)=co⁡{lim⁡∇f(xi):xi⟶x,xi∉S,xi∉Ωf}.
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 functionHu(t,·,λ) is discontinuous in u=0, and that if u≠0, then Hu+≡Hu-, we have that the partial Clarke’s gradient is (2.13)∂uH={A⋅π(t)+λifu>0[Hu+,Hu-]ifu=0η⋅A⋅π(t)+λifu<0
with Hu+=A·π(t)+λ and Hu-=η·A·π(t)+λ. The switching function is Φ(z,λ)≡∂uH, and the optimality condition leads to (2.14)u*(t)={umax⁡=0ifA⋅π(t)+λ>0using=0if0∈∂uH=[A⋅π(t)+λ,η⋅A⋅π(t)+λ]umin⁡=0ifη⋅A⋅π(t)+λ<0.
On the other hand, the costate equation of Pontryagin’s Maximum Principle allows us to obtain (2.15)λ′=-Hz=0⟶λ=λ0(const).
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 (2.16)u*(t)={umax⁡=0ifA⋅π(t)>-λ0using=0if-λ0∈[A⋅π(t),η⋅A⋅π(t)]umin⁡=0ifη⋅A⋅π(t)<-λ0.


## 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.Figure1 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,π(t) must be interpolated to obtain a continuous function.(ii)
Second, for a givenλ, we have to determine the switching times: t1,t2,…. These instants are calculated solving
(3.1)A⋅π(t)=-λ;η⋅A⋅π(t)=-λ.(iii)
Third, the optimal valueλ0 must be determined in order for
(3.2)zλ(T)=∑i=1Nuδiu⋅qmax⁡+∑i=1Nlδil⋅qmin⁡=b
with δiu and δil being the duration of the ith bang-bang segment in the upper and lower bound, respectively, Nu and Nl the number of such segments, and zλ(T) the final volume obtained for each λ.(iv)
To calculate an approximate value ofλ0, 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(4.1)P(t)=0.000126821z′(t)
with A=0.000126821. We shall also consider the technical constraints: qmin=-283866 (m3/h), qmax=394258 (m3/h). When the efficiency is η=1.15, these constraints, respectively, correspond to Pmin=-41.4,Pmax=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 ofT=24 h. was considered. The clearing price, π(t) (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, π(t), is only known at each hour (t=1,2,…,24). In this paper the known values of π(t) were linearly interpolated with good results.Table 1
The clearing price,π(t).
tπ(t)tπ(t)tπ(t)tπ(t)176.93769.4713104.081990.00268.20875.7914100.0020106.89368.209105.901580.5021103.00460.0010106.501678.2322100.00555.0111110.001775.932386.93656.2812108.461878.232479.93The 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:b=2·106 (m3) 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η.
ηProfit (euros)Pumped water (m3)% improvementwithout pumping27145.2001.3530282.51.49123·10611.61.3030896.41.61463·10613.81.2531567.51.7438·10616.31.2032300.1.87975·10618.91.1533105.52.07863·10621.9Table 3
Optimal solution varyingb.
b (m3)Profit (euros)Profit (euros)% improvement(without pumping)(with pumping)1·10613753.120633.950.02·10627145.23230018.93·10640067.643318.38.14·10652017.653733.13.3As 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 around80%, the improvement in profit is quite substantial, being around 20%.In Figure2 we present the solution of the algorithm when the efficiency is η=1.20. The optimal value obtained for -λ0 is 0.01139601030789622, and the switching times are (4.2)t1=1.2345,t2=7.85646,t3=8.46727,t4=14.52,t5=18.9881,t6=22.7759.
Figure 3 presents the optimal hydro-power, P, for three cases: (a) without pumping, (b) with η=1.2, and (c) with η=1.15.Figure 2
Value of-λ0.Optimal hydro-power,P(t).
(a)(b)(c)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 π(t). 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 π(t) 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 Figure4). 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 (4.3)Error=|zλ(T)-b|<tol
with tol=50 (m3). The secant method has provided satisfactory results using these initial values: (4.4)λmin⁡=min⁡Aπ(t);λmax⁡=max⁡ηAπ(t).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 b. We compare with the solution without pumping.As can be seen (see Table1), 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·106 (m3) 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.


## 4.1. Influence of the Efficiency

First, we consider the restriction on the volume:b=2·106 (m3) 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η.
ηProfit (euros)Pumped water (m3)% improvementwithout pumping27145.2001.3530282.51.49123·10611.61.3030896.41.61463·10613.81.2531567.51.7438·10616.31.2032300.1.87975·10618.91.1533105.52.07863·10621.9Table 3
Optimal solution varyingb.
b (m3)Profit (euros)Profit (euros)% improvement(without pumping)(with pumping)1·10613753.120633.950.02·10627145.23230018.93·10640067.643318.38.14·10652017.653733.13.3As 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 around80%, the improvement in profit is quite substantial, being around 20%.In Figure2 we present the solution of the algorithm when the efficiency is η=1.20. The optimal value obtained for -λ0 is 0.01139601030789622, and the switching times are (4.2)t1=1.2345,t2=7.85646,t3=8.46727,t4=14.52,t5=18.9881,t6=22.7759.
Figure 3 presents the optimal hydro-power, P, for three cases: (a) without pumping, (b) with η=1.2, and (c) with η=1.15.Figure 2
Value of-λ0.Optimal hydro-power,P(t).
(a)(b)(c)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 π(t). 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 π(t) 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 Figure4). 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 (4.3)Error=|zλ(T)-b|<tol
with tol=50 (m3). The secant method has provided satisfactory results using these initial values: (4.4)λmin⁡=min⁡Aπ(t);λmax⁡=max⁡ηAπ(t).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 b. We compare with the solution without pumping.As can be seen (see Table1), 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·106 (m3) 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.



---

*Source: 561592-2011-08-09.xml*