Eﬀicient exact A* algorithm for the single unit hydro unit commitment problem

—The Hydro Unit Commitment problem (HUC) specific to hydroelectric units is part of the electricity production planning problem, called Unit Commitment Problem (UCP). More specifically, the studied case is that of the HUC with a single unit, denoted 1-HUC. The unit is located between two reservoirs. The horizon is discretized in time periods. The unit operates at a finite number of points defined as pairs of the generated power and the corresponding water flow. Several constraints are considered. Each reservoir has an initial volume, as well as window resource constraints, defined by a minimum and maximum volume per time period. At each time period, there is an additional positive, negative or zero intake of water in the reservoirs. The case of a price-taker revenue maximization problem is considered. An efficient exact A* variant, so called HA*, is proposed to solve the 1-HUC accounting for window constraints, with a reduced search space and a dedicated optimistic heuristic. This variant is compared to a classical Resource Constrained Shortest Path Problem (RCSPP) algorithm and a Mixed Integer Linear Programming formulation solved with CPLEX. Results show that the proposed algorithm outperforms both concurrent alternatives in terms of computational time in average on a set of realistic instances, meaning that HA* exhibits a more stable behavior with more instances solved.


I. INTRODUCTION
An electricity producer aims at meeting the demand at any time.This is because electricity can hardly be stored, meaning that any excess is lost.Scheduling the short-term production in order to meet the demand defines the Unit Commitment Problem (UCP), which is solved for the day ahead.At Electricité de France (EDF), large-scale UCP instances are solved by Lagrangian decomposition [1], yielding subproblems of the same nature (thermal, hydraulic, solar,...).As each sub-problem has different constraints, specific approaches are developed for each of them.Among these sub-problems, the Hydro Unit Commitment (HUC) has received a lot of attention, due to the large size of its instances.Indeed, instances of the HUC involve valleys, which can be constituted of up to twenty units, linked with reservoirs.
At EDF, the HUC is modeled as a Mixed Integer Linear Program (MILP) and solved with CPLEX [2].This MILP considers operating points, which are pairs (water flow, corresponding generated power).In practice, the HUC is not solved to optimality within the time limit set, as such representation induces an exponential number of solutions and large computational times.More recent work [3] have pointed out the interest of solving the HUC using a Lagrangian relaxation algorithm, where sub-problems are single unit HUC (1-HUC).The main benefit of this relaxation is that the 1-HUC can be solved with dynamic programming, while the master problem handles the coupling constraints.It is shown that such relaxation can lead to overall better results than solving the HUC as an MILP, which emphasizes the relevance of an efficient algorithm to solve the 1-HUC.
In this paper, we consider the 1-HUC with a single unit located between two reservoirs.A diagram, taken from [4] and shown in Figure 1, is sketching the 1-HUC.The principle of hydroelectric production is the following: the water from the upstream reservoir flows into the downstream reservoir through the unit, thus driving the turbines of the unit, which in turn power the generator to produce electricity.The unit operates on M turbining points, N pumping points and an idle operating point.With I = {−N, . . ., 0, . . ., M }, each operating point i ∈ I is defined as a pair formed by a water flow D i and a generated power P i .Both D i and P i are positive (resp.negative) for turbining (resp.pumping) operating points, i.e., with i > 0 (resp.i < 0), and are 0 for the idle operating point i = 0.The operating points are defined in a cumulative fashion meaning that if a unit is at turbining (resp.pumping) operating point i, then order constraints apply, involving all points 1 ≤ j < i (resp.−1 ≥ j > i) to also be operated.At each time period, the unit cannot turbine and pump simultaneously.The time horizon is discretized into T time periods.At each time period t, the unit turbines (resp.pumps) a water flow and produces (consumes) an amount of energy that is considered to be constant for the duration of the time period.Resource window constraints state that the volume of each reservoir n ∈ {1, 2} lies between a lower bound V n t and an upper bound V n t that are time-dependent.At each time period, the reservoir n receives an additional intake of water A n t .The additional intake can be positive to represent rain, melting snow etc, or negative to represent the use of water for local agriculture etc.
The revenues take into account the unit value Φ n of the water in each reservoir n at the end of the time horizon, and the value of the energy produced or consumed at a time-dependent unit value Λ t .The problem is to maximize the total revenue, while satisfying the reservoir capacities at each time period.
In practice, water management policies require that reservoirs should meet target volumes (V n t = V n t − ϵ) at the end

Downstream reservoir Unit
Fig. 1: Diagram of the 1-HUC of the time horizon.Due to these target volumes and the finite set of operating points, the 1-HUC may not have any feasible solution.However, it is possible to adjust [5] or relax [6] the target volumes in order to obtain feasible instances.Hence, we consider in this paper only 1-HUC instances admitting feasible solutions.
In this paper the aim is to propose a dynamic programming algorithm dedicated to the 1-HUC, modeled as a Resource Constrained Shortest Path Problem (RCSPP) with resource window constraints.The algorithm we propose is an exact variant of the A* algorithm [7] used to compute the shortest path in a graph.To obtain an efficient algorithm, the bounds of the windows are first tightened by propagating them from any time period to another, and a dedicated heuristic is developed.As we propose an exact algorithm, we compare it with two other exact methods, namely an MILP solved by default CPLEX as currently done at EDF and a classical RCSPP algorithm.The corresponding evaluation of performance is done on various realistic instances.The numerical results show that our algorithm yields smaller computational time variations compared to both the MILP formulation solved by CPLEX and the RCSPP algorithm.
The remainder of the paper is organized as follows.In Section II, a literature review of dynamic programming algorithms for related problems is reported.In Section III, a formulation of the 1-HUC is proposed.In Section IV, graphical representations of the 1-HUC are provided.In Section V, the bound tightening procedure and the exact A* variant are described.In Section VI, numerical results are presented.In Section VII, concluding remarks and perspectives for further research are drawn.

II. STATE OF THE ART
In this section, we first present a literature review of dynamic programming algorithms developed for the UCP and for the HUC.As the HUC can be represented as an RCSPP, we also include in this review dynamic programming algorithms for the RCSPP.

A. Dynamic programming for the Unit Commitment Problem
A dynamic programming algorithm for a single Unit Commitment (1-UC) with ramp and min up/down constraints is presented in [8].The algorithm is based on a graph with a source vertex and several groups of T vertices.For each even (resp.odd) group, vertex t indicates that the unit is turned off (resp.on) at time period t.The arcs connect the vertices of a group to the next groups, from a time period t to a time period t ′ > t.Finding a path in this graph allows one to find an onoff schedule for the unit.The difference between the 1-UC and the 1-HUC is that there is no resource in the 1-UC, while in the 1-HUC, the presence of reservoirs with minimum and maximum volumes requires to account for the water resource.

B. Dynamic programming for the Hydro Unit Commitment
In this part we focus on dynamic programming algorithms for the HUC, most of them being cited in the survey [9].
In [1] the author presents a two phase approach to solve the HUC, solving an LP for the first phase and using a dynamic programming algorithm for the second phase.More precisely, the second phase consists of solving the 1-HUC for each unit of the valley with dynamic programming, aiming to get the closest solution to the LP solution while taking into account constraints omitted in the LP.For this phase, the considered underlying graph is as follows.The reservoir volume is discretized, yielding hundreds of possible volume values for a reservoir.Then, a vertex is defined for each volume value and each time period, thus leading to hundreds of vertices per time period.A Bellman-Ford algorithm [10] is used to find a path in this graph.Such a discretization discards a lot of realistic states.In this paper, we consider all possible states with respect to the operating points, which can be exponential for instances with a large number of time periods.With an exponential number of states, the Bellman-Ford algorithm becomes far less efficient.
In [6] a method for solving a non-linear 1-HUC with a target volume is described.To solve this problem with dynamic programming, a state diagram is constructed.In a similar fashion as in [1], evenly discretized volumes are considered, yielding a limited number of states per time periods.In order to have feasible solutions, the target volume is relaxed to match this discretization.The state diagram is constructed by generating the possibilities to reach the target volume from the initial volume, satisfying the upper and lower bounds on the volume at each time period.Starting from the state at the end of the time horizon, the dynamic programming algorithm maximizes the value of the generated power.As we consider 1-HUC instances with and without target volumes, a backward algorithm may not be practical with large volumes and no target volume.
In [3], a decomposition method for solving the HUC with shortest paths is described.The considered HUC is a valley where each unit has a finite number of operating points.The topology of the valley is not restricted to a chain, as each unit (resp.reservoir) can have a set of upstream and downstream reservoirs (resp.units).There are additional ramping constraints, namely the flow variation is limited from one time period to another.This HUC also takes into account a target volume for each reservoir, at the last time period.Note that the latter target volume is a minimal bound, meaning there is no equality constraint.The solution approach decomposes this HUC into multiple 1-HUC.Some 1-HUC without resource constraints are solved by a shortest path algorithm, while others with resource constraints are solved by a labeling algorithm defined in [11].The latter algorithm is adapted from a classical RCSPP algorithm [12] to take into account a minimum bound for the resource.It is mentioned that this labeling algorithm loses its dominance properties between two labels if one of them does not verify the minimum bound on the resource.Such case is more frequent when the target volume is considered with equality constraints, making this algorithm less efficient.
There are also other problems solved by dynamic programming, related to the 1-HUC.In [13] a dynamic programming approach is described to solve an HUC on instances of the Itaipù unit (Brazil, Paraguay).This problem differs from ours as the only constraint is to satisfy the minimum and maximum number of turbines running at each time period.No volume is considered, therefore there are neither bounds on the volume, nor a target volume.In [14] the Hydro Unit Load Dispatch problem (HULD) is presented.This problem differs from the HUC, as the water flow is known, and solving the HULD is to provide the most economic distribution of the water through the different turbines, while verifying the flow capacity of the turbines at each time period.

C. Shortest path with resource constraints
As the HUC can be seen as a shortest path problem with a water resource bounded both from below and above, we are interested in the solution methods for the RCSPP (Resource Constrained Shortest Path Problem).
There are works on the RCSPP to solve the thermal problem on EDF instances [15].In that paper it is indicated that the resource has an upper bound but no lower bound.However, as specified in [3], the difficulty of the HUC comes from the lower bound on the volume, which prevents the use of dominance rules.
In the survey [16], a state of the art review of different shortest path variants is described.More specifically, it is indicated that there is little work on the RCSPP with equality constraints, or window constraints (such as in the HUC).Three papers are cited, namely [17] describing a heuristic, [18] presenting an integer formulation and [19] proposing a dynamic programming algorithm.As we look for an exact algorithm, we will focus on the three-phase algorithm described in [19].The presented algorithm solves the RCSPP with window constraints on acyclic graphs.The main idea, further detailed in [20], is to extend the graph, such that if multiple paths lead to the same vertex from the source, a new vertex is created for each of these paths.Once the graph has been extended in this way, the problem is solved with a pseudo-polynomial time algorithm.Such a graph extension seems difficult to apply to the 1-HUC.Indeed, graphical representations of the 1-HUC (see Section IV) solely have arcs from time period t to time period t + 1.Hence, the extension can lead up to M T vertices at time period T .Even a pseudo-polynomial algorithm on the extended graph could be impractical in the case of the 1-HUC.

III. INTEGER LINEAR PROGRAMMING
With x t,i the binary variable indicating whether the unit is at least at operating point i ∈ I at time period t ≤ T , we obtain the following formulation: max T t=1 i∈I In this formulation, the objective function maximizes the total revenue.Constraints (a1) to (b2) ensure that the minimum/maximum bounds on the volume are verified for both the upstream and downstream reservoirs at each time period.Constraints (c) and (d) correspond to the order of the operating points.Constraints (e) prevent the unit from turbining and pumping simultaneously.Lastly, constraints (f) indicate that all variables x t,i are binary.

A. Shifting all operating points
In the following, it will be convenient to only have operating points with non-negative power and flow.In [3] a modification on the flows and the volume bounds is considered to only have such operating points.First, each turbining operating point i, Operating point 0 remains unchanged.As such, all the operating points have non-negative cumulated flows.For a given t ≤ T , the bounds on the volume V and V 1 t (resp.V 2 t and V 2 t ) are shifted in order to keep the same feasible solutions: Consider an instance of the 1-HUC with M = 3 turbining operating points (D 1 = 3, P 1 = 3), (D 2 = 4, P 2 = 3), (D 3 = 2, P 3 = 2) and N = 2 pumping operating points and downstream reservoir bounds are V 2 = [60, 60, 60] and The renumbering is such that the 3 turbining operating points, of index 1 to 3 are now of index 1 + N to 3 + N , i.e., From pumping operating point −1 a turbining operating point N − 1 + 1 = 2 is created, and similarly from the pumping operating point −2 operating point ).The upstream reservoir upper bounds are modified as follows In the following, we only refer to the 1-HUC with operating points of non-negative flow and power, i.e., with I = {0, . . ., M }.In this case, only constraints (a1) to (b2), (c) and (f) are necessary to model the constraints of the 1-HUC, and no renumbering is required.

B. Rewriting the formulation
The formulation defined previously is a classical MILP model for the 1-HUC.We rewrite the formulation in order for the window constraints to appear more clearly.Constraints (a1), (a2), (b1) and (b2) can be rewritten as follows, considering There are redundancies between (a1') and (b2'), and between (a2') and (b1').Let us introduce bounds β t ′ and α t ′ in the following way with t ′ ≤ T : By using β t ′ and α t ′ , and rewriting the objective function, we obtain the formulation F 1HU C defined as follows: With this formulation, the objective function is to maximize the value of each active operating point, plus a constant.Also, one can see (a) and (b) as resource window constraints, or (a) as nested cover constraints and (b) as nested knapsack constraints.This formulation can also be improved by tightening the window constraints, as shown at the end of Section V-A.

IV. GRAPHICAL REPRESENTATION
The 1-HUC can be represented graphically in two ways, namely in a fashion similar either to the knapsack problem, or to the RCSPP.In the following we describe both representations and their dominance rules.We first introduce the cumulated flows between two time periods.
Definition 1 (Cumulated flow D t ′ ,t ).The cumulated flow D t ′ ,t is the sum of the flows from time periods t ′ to t:

A. Representation as a knapsack problem
Let G KP = (V KP , A KP ) defined as follows.Each vertex u ∈ V KP is defined as a pair u = (t, d) with t the time period, and d the cumulated flow D 1,t with variables x associated to a path that reaches u.Without loss of generality, u = (t, d) is considered only if d ∈ [β t ; α t ].The source vertex s is defined as s = (0, 0).For each vertex u = (t, d) and v = (t + 1, d + i j=0 D j ) with t < T and i ∈ {0, . . ., M } there is an arc The downside of such graph is its exponential number of vertices.However, by tightening the bounds as shown in Section V-A, it is possible to drastically reduce the number of vertices in G KP .Furthermore, we can use classical dominance rules for the longest path: Definition 2 (Dominance rule 1).Let p and q be two paths from s to a vertex u.By induction the path with the lowest value is dominated, as it cannot lead to an optimal solution.Definition 3 (Dominance rule 2).Let p be a path from s to u going through a vertex v, and q be a path from s to v. Let p s,v be the subpath of p, from s to v and p v,u the subpath of p from v to u.If the value of p s,v is larger than that of q, then q is dominated.If the value of q is larger than that of p s,v , then p is dominated by the path concatenating q and p v,u .

B. Representation as an RCSPP
Let G R = (V R , A R ) be defined as follows.Each vertex u ∈ V R is defined as a pair u = (t, i), with t the time period and i the operating point.The source s is defined as s = (0, 0).From each vertex u = (t, i), with t < T , and v = (t + 1, i ′ ) with i ′ ∈ {0, . . ., M } there is an arc The downside of G R , is that there are paths in this graph which do not verify the resource constraints, hence one needs to verify the resource constraints for each path.Note that the resource in this case is D 1,t for variables x associated to the path.On the positive side, the number of vertices is polynomial, and it is possible to use a classical dominance rule for the RCSPP as defined in [3]: Definition 4 (Dominance rule 3).If there are two partial paths from s to the same vertex u, by induction if a path has a lower value and uses more resource than another one, then this path is dominated provided both partial paths use sufficient resource to verify all lower bounds on the resource.
Note that the condition on the resource usage to ensure that the lower bounding constraints (a) are satisfied seriously weakens the dominance rules when these constraints are active.There are two other graphical representations for the 1-HUC described in the literature.In [20], the original graph defined for the RCSPP is similar to the one depicted in Fig. 2b.However, as mentioned in Section II-C, this graph is extended such that if multiple paths exist from vertex s to a vertex u, then u is duplicated such that only a single path leads to each vertex.The extended graph would be larger than the one depicted in Fig. 2a, as even with an exponential number of vertices, there are still multiple paths between s and most of the vertices.In [6], as described in Section II-B, the volume is evenly discretized.The resulting graph (3,10) (a) Graphical representation as a knapsack problem (b) Graphical representation as an RCSPP Fig. 2: Graphical representations of the 1-HUC has similar vertices as in Fig. 2a, but at each time period, vertices only exist for volumes within the set of discretized volumes.In the experimental results of [6], it is stated that the discretization is ranging from 0.3% to 0.5% of the difference between the minimum and maximum volume of the reservoir, which represents about 300 vertices at each time period.Such graph heavily differs from the graph we consider, as there can be a very large number of vertices at each time period.

V. EXACT A* VARIANT FOR THE 1-HUC
In this section, we describe the new algorithm proposed to solve the 1-HUC.The aim of this algorithm is to find the longest path in graph G KP as described in Section IV-A.A difficulty of this graph is the exponential number of vertices, which is why we resort to a variant of the A* algorithm [7].The A* algorithm is particularly efficient when the number of vertices is large as it involves an optimistic heuristic to guide the search, and to discard sub-optimal partial solutions.We denote the proposed exact variant of the A* algorithm for the 1-HUC by HA*.This algorithm involves a dedicated optimistic heuristic for the 1-HUC.To further improve the performance of this algorithm, we also present a pre-processing bound tightening procedure in order to reduce the number of vertices in G KP .
In the following, we first present the bound tightening procedure, then the optimistic heuristic, i.e., an upper bound on the optimal value and finally the HA* algorithm.

A. Tightening the bounds
The bounds α t and β t from inequalities (a) and (b) can directly be used to prune infeasible vertices.Thus, tightening these bounds may lead to a smaller number of vertices, thus reducing the search space.
In order to tighten the bounds, we use the cumulated flow D t ′ ,t as previously defined.At each time period, the flow is between 0 and M i=1 D i .For any pair of time periods (t ′ , t), t ′ < t, the lower bound D t ′ ,t = 0 and the upper bound D t ′ ,t = (t − t ′ + 1) M i=1 D i will never be violated by a feasible solution.Note that α t and β t are not bounds on the flow at time period t, but rather bounds on D 1,t .If the gap between β t and α t is large, then one may develop a large number of vertices, sometimes leading to intractable instances.For example, in the case of the HUC, it is very common to have upper bounds α t very large compared to the flows, and negative lower bounds β t .We can therefore introduce bounds αt and βt in the following way: By using bounds αt and βt , we can drastically reduce the number of vertices.However, it is still possible to further reduce it.Suppose that at time period t the bounds are such that βt > βt+1 .As the water flows are all non-negative, any vertex u = (t + 1, d) with d < βt + D t+1,t+1 cannot be part of a feasible solution.Similarly, if βt−1 < βt , any vertex u = (t − 1, d) with d < βt − D t,t cannot be part of a feasible solution.Extending this logic to the upper bounds on d, we can tighten the bounds of any time period from the bounds of any other time period, following the rules below.Let a pair of time periods (t ′ , t) with t ′ < t.Then, D 1,t must stay in Let us define αt and βt as follows: Tight bounds α * t and β * t are calculated as follows: Note that computing all bounds β * t is of complexity T 2 .Indeed, for a given t, computing β t as well as βt both require one comparison, computing βt needs T comparisons and computing β * t requires one comparison.Hence, computing all β * t is of complexity T 2 .We obtain a similar complexity for upper bounds α * t .

B. Optimistic heuristic
In the case of the 1-HUC, an optimistic heuristic overestimates the value of the objective function because we are solving a maximization problem.Let p be a path from time period 1 to t representing already taken decisions.The heuristic aims at computing an optimistic cost from time period t+1 to T .The idea of the proposed optimistic heuristic is to compute an improved linear relaxation on time periods t + 1 to T .To do so, we define quadruplets (t, i, val, f low), with t a time period, i an operating point, val the value of any arc (t−1, d), (t, d+ i j=1 D j ) in G * KP , and f low the value i j=1 D j .The aim is to progressively increase the values of variables x t,i depending on their profitability, being val/f low.
Algorithm 1 describes how to compute a linear relaxation on time periods t + 1 to T from a partial path in graph G * KP , as detailed in the following four steps.
Step 1: Initialize a fractional solution with x t ′ ,i = 1 if the partial solution represented by path p requires operating point i ≤ M at time period t ′ ≤ t, x t ′ ,i = 0 otherwise.This step is represented by lines 1 to 8 of Algorithm 1 Step 2: Initialize a list with all the quadruplets at time period t ′ ∈ [t + 1; T ] by decreasing profitability val/f low.This step is represented by lines 9 to 14 of Algorithm 1 Feasibility step 3: This step is repeated as long as lower bounding constraints (a) are not verified (the while loop at line 15 of Algorithm 1).The algorithm looks for the smallest time period t ′ such that (a) is not satisfied, and for x t ′′ ,i with t ′′ ∈ [t+1; t ′ ] maximizing profitability val/f low.The variable x t ′′ ,i is increased depending on its profitability: • If the profitability is positive, fractionally increase x t ′′ ,i as much as possible provided all upper bounds α are satisfied.• Otherwise, fractionally increase x t ′′ ,i as little as possible provided all upper bounds and the lower bound β t ′ are satisfied.
All variables x t ′′ ,i ′ < x t ′′ ,i with i ′ < i must be set to the same value as x t ′′ ,i in order to satisfy the order constraints.Also, for any i ′ > i, quadruplets (t ′′ ,i ′ ,val ′ ,f low ′ ) are updated as (t ′′ ,i ′ ,val ′ − val,f low ′ − f low).This is because as x t ′′ ,i > x t ′′ ,i ′ , one can increase x t ′′ ,i ′ without increasing x t ′′ ,i while still verifying order constraints.All variables that cannot be further increased due to the upper bounds are removed from the list.Optimality step 4: This step is repeated as long as there is a variable of positive profitability in the list of variables (again the while loop at line 15 of Algorithm 1).Select the first variable of the list, and fractionally increase its value as much as possible provided all upper bounds are satisfied.Remove from the list all variables that cannot be further increased due to the upper bound.
Property 1.The fractional solution returned by Algorithm 1 verifies order constraints.
Proof.Consider two quadruplets (t,i,val,f low) and (t,i ′ ,val ′ ,f low ′ ), with t a time period considered in the heuristic, and i ′ > i.At the start of the algorithm, x t,i = x t,i ′ = 0.If x t,i ′ is increased, then x t,i is increased by the same amount, hence order constraints are verified.If x t,i is increased, it means that val/f low > val ′ /f low ′ .Hence, val/f low > (val ′ − val)/(f low ′ − f low).Consequently, the algorithm will increase x t,i ′ only if x t,i = 1, hence order constraints are verified.
Theorem 1. Algorithm 1 defines an optimistic heuristic.Proof.Let s be an integer solution for the 1-HUC, for which variables x t,i verify constraints (a) (b) (c) and (f ).Let ŝ be a fractional solution for the 1-HUC obtained with Algorithm 1, for which xt,i verify all constraints (a), (b) (c), and constraint (f ) only for time periods 1 to t ≤ T .Consider ŝ and s to be identical for time periods 1 to t.
Let X (resp.Y) be the variables such that Clearly, for each variable x t ′ ,i ∈ X of positive profitability, a fractional value for x t ′ ,i increases the value of the objective function compared to x t ′ ,i = 0. Similarly, for each variable in Y of negative profitability, a fractional value for x t ′ ,i increases the value of the objective function compared to x t ′ ,i = 1.
Let X − ⊆ X and Y − ⊆ Y be the variables with negative profitability.Suppose |X − | > 0. By construction of ŝ, the variables of X − have value greater than 0 only in order to yield a feasible solution, with respect to a lower bound β t ′ .In such case the total flow of ŝ from time period 1 to t ′ is exactly β t ′ by construction of ŝ.As solution s also verifies all lower bounds, we deduce |Y − | > 0. Otherwise there is either a contradiction with the construction of ŝ, or s does not verify all lower bounds.Hence, the total flow of s from time period 1 to t ′ is at least β t ′ .By definition, s and ŝ are identical from time period 1 to t, consequently the only difference is on time periods t + 1 to t ′ .By construction of ŝ the variables of X − are the most profitable and have fractional value in ŝ.Consequently, their weighted value in the objective function must be higher than those of Y − in the integer solution s.
A similar proof can be made for variables in Y + ⊆ Y and X + ⊆ X the variables with positive profitability.
The value of ŝ is then greater or equal to the value of s.
It is possible to tighten the fractional solution returned by Algorithm 1 while keeping it optimistic.Clearly, an integer solution is necessarily of a total flow which is a combination of the flows from the operating points.Therefore the flow of an integer solution is necessarily a multiple of the greatest common divisor (GCD) of the operating points' flows.When the heuristic increases the value of a variable, we can increase or reduce this value so that the total flow of the returned solution remains a multiple of the GCD relative to the operating points' flows.Note that since the flows are identical from one time period to another, we can quickly compute the GCD by considering only the flows of a single time period.

C. HA* algorithm
For a 1-HUC with an objective function to maximize, the principle of HA* is the following.Consider a pool of partial solutions evaluated with the heuristic.At each iteration, the partial solution with the highest heuristic value is considered and removed from the pool.From the partial solution considered, we complement it by adding neighbors relative to its last vertex.Once a solution is found, its value is used as a bound to remove some more partial solutions from the pool.Indeed, if the solution's value is higher than a partial solution's heuristic val the value of an arc towards (t end for 32: end while 33: return the value of x value, then the partial solution can be removed from the pool.
Once the pool of partial solutions is empty, the algorithm stops and the best solution found is the optimal solution.We underline the need of a tight optimistic heuristic.If the heuristic is not optimistic, or optimistic but too loose, there are fewer cases where one can prune partial solutions while guaranteeing optimality.Hence, more vertices are developed which can exponentially increases the computational time.
For readability purposes, we introduce three structures: Definition 5 (Path structure).The path structure has three attributes: vertices the list of vertices of the path; val the value of the path with respect to the objective function; heur the optimistic heuristic value.
Definition 6 (Vertex structure).A vertex structure has two attributes: t the time period and d the cumulated flow D 1,t , as defined for the vertices of G * KP in Section IV.
Definition 7 (Arc structure).The arc structure has one attribute: val the value as defined for the arcs of G * KP in Section IV.
Algorithm 2 presents the pseudo-code of HA*, using the three previously described structures as well as OptimisticHeuristic.The considered graph is G * KP , which is G KP as illustrated in Fig. 2a  q.vertices = p.vertices ∪ u, q.val = p.val + a.val, q.heur = OptimisticHeuristic(q) if |q.vertices| = T + 1 then bestV al ← min(bestV al, q.val) remove q ′ ∈ L p with q ′ .val+ q ′ .heur≤ bestV al else dom ← F ALSE for q ′ ∈ L p do if q ′ dominates q then dom ← T RU E end if if q dominates q ′ then remove q ′ from L p end if end for if dom = F ALSE and q.val + q.heur > bestV al then add q in L p by keeping L p sorted by decreasing q.val + q.heur end if end if end for end while return the solution of value bestV al

VI. EXPERIMENTAL RESULTS
Results are computed on a single thread of an Intel(R) Core(TM) i7-9850H CPU @ 2.60GHz processor, with 2 CPUs of 8 cores, with Linux as operating system.All algorithms are developed with C++.Version 12.8 of CPLEX with default setting is used to solve the MILP formulation F 1HU C .

A. Instances
From a large set of realistic instances derived from a real EDF unit, a first set A of 13 instances is obtained.These 13 instances are retained as preliminary results have shown that formulation F 1HU C is not trivially solved.This emphasizes the need of an efficient alternative in these cases.Table II depicts for each instance the main characteristics, namely the number of time periods, the number of operating points and the presence of a constraining minimum (resp.maximum) bound on the volume at the last time period.The instances cover three cases, namely when there is only an upper bound, only a lower bound, or a target volume with a constraining upper and lower bound.In the case of an equality constraint, the instance becomes infeasible, which is out of the scope of the instances considered in this paper.Hence, the target volumes are not equality constraints, but rather tight window constraint.For instances with a target volume, more resource window constraints are obtained by propagating the bounds β * and α * from the bounds at the last time period to the previous ones.
The water flows D and powers P are in the order of 10 3 to 10 4 , with volumes in the order of 10 7 .For target volumes, the difference between the upper and lower bound is in the order of 10 3 , which is small enough with respect to the flows to yield very few vertices at time period T in G * KP .A second set B of 13 instances, similar to the first set, is also constructed.The only differences are bounds β T and α T that are shifted as follows.Let an instance with bounds β A T and α A T be in set A. A random value k ∈ [−9999; −1000] ∪ [1000; 9999] is chosen.Bounds of an instance for set B are β B T = β A T + k and α B T = α A T + k.Note that this shift is very small for these instances.Indeed, the water flows can be in the order of 10 4 , there are nearly 100 time periods and the cumulated flow D 1,T is in the order of 10 6 .The shift k is at most 1% of D 1,T .As shown in the following results, slightly modifying an instance can drastically impact the computational time.

B. Experimental results
In order to benchmark HA*, all instances are solved with HA* as well as with two alternative methods.The first alternative is a classical RCSPP algorithm [12] adapted to the 1-HUC [11].The second alternative is to use CPLEX to solve F 1HU C described in Section III.A third alternative is solving F + 1HU C , which is formulation F 1HU C described in Section III with tighter bounds described in Section V-A.However, solving F + 1HU C leads to very similar results as F 1HU C , hence results for F + 1HU C are not included in the following.All algorithms use a single thread, with a time limit of one hour.

Example 3 .
Consider the instance of Example 2. Fig.2brepresents the graph G R associated to this instance.

Fig. 3 :
Fig. 3: Total number of instances A and B solved per computational time

TABLE I :
Reducing the search space using bounds on the flows (a) Table with bounds αt and βt Table with bounds αt and βt Table with bounds α * t and β * The operating points are such that at each time period, the maximum flow is 2. By applying tighter bounds we can see that we drastically reduce the possibilities, thus the number of vertices potentially developed by dynamic programming.In TableIathe invalid values for the total flow at each time period, with respect to bounds β t , are marked with a cross.Table Ib is similar to Table Ia with tighter bounds αt and βt , the crosses being in bold to emphasize the tightening of the bounds.Table Ic follows the same representation with the tightest bounds β * t and α * t .As previously mentioned, formulation F 1HU C can be improved, considering bounds β * t and α * t instead of β t and α t in constraints (a) and (b).We denote by F + 1HU C the formulation F 1HU C with the bound tightening.In the following, we consider the graph G KP with the tightest bounds β * t and α * t , denoted by G * KP .
Algorithm 1 Algorithm OptimisticHeuristicRequire: A path p from time period 1 to t, a graph G * KP , GCD the GCD of the flow 1: Initialize a fractional solution x with all variables to 0 2: for t ′ ∈ [1; T ] and i ∈ [1; M ] do 3:if t ′ ≤ t AND p requires operating point i at time period t ′ then ′ , i) in G * ′ , i, val, f low) in L,sorted by decreasing val/f low 14: end for 15: while ∃t ′ ∈ [t + 1, T ] such that x does not verify β t ′ OR exists quadruplet in L with val > 0 do (t ′′ , i, val, f low) ← first in L with t ′′ ≤ t ′ 17:if val ≤ 0 then 18:set xt ′′ ,i to minimum such that β t ′ verified andx t ′′ ,i •f low mod GCD = 0; if β t ′ cannot be verified xt ′′ ,i ← 1set xt ′′ ,i to maximum such that all upper bounds are verified andx t ′′ ,i • f low%GCD = 0 xt ′′ ,i ′ = max(x t ′′ ,i ′ , xt ′′ ,i ) for (t ′′ , i ′ , val ′ , f low ′ ) ∈ L with i ′ ∈]i; M ] do 26: val ′ ← val ′ − val 27: f low ′ ← f low ′ − f low for (t ′′ , i, val, f low) ∈ L such thatxt ′′ ,i cannot be increased do 30: remove (t ′′ , i, val, f low) from L 31: with bounds β * t and α * t .The dominance rules used in Algorithm 2 are the dominance rules 1 and 2. KP Initialize a path p as follows: p.vertices = {(0, 0)}, p.val = 0.0, p.heur = OptimisticHeuristic(p) Initialize a list of paths L p = [p] Initialize the value of the best solution bestV al = −∞ while L p not empty do p ← first path in L p remove p from L p v ← last vertex of p.vertices for arc a from v to u with u.d ∈ [β t ; α t ] do

TABLE II :
Main characteristics of instances set A