Pierre Carpentier
email: pierre.carpentier@ensta-paristech.fr

Jean-Philippe Chancelier

Vincent Leclère
email: vincent.leclere@enpc.fr

François Pacaud
email: f.pacaud@efficacity.com

J.-Ph Chancelier

Stochastic decomposition applied to large-scale hydro valleys management

Keywords: Discrete time stochastic optimal control, Decomposition methods, Dynamic programming, Energy management 2010 MSC: 90-08, 90C06, 90C15, 90C39, 49M27

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Large-scale systems and energy applications

Consider a controlled dynamical system over a discrete and finite time horizon. This system may be influenced by exogenous noises that affect its behavior. Assume that, at every instant t, the decision maker designs a control based on all the observations of noises available up to time t. We are thus looking for strategies (or policies), that is, feedback functions that map every instant and every possible history of the system to a decision to be made.

We can find typical applications in the field of energy management. Consider a power producer that owns a certain number of power units. Each unit has its own local characteristics such as physical constraints that restrain the set of feasible decisions, and induces a production cost or a revenue. The power producer control the power units so that an overall goal is met. A classical example is the so-called unit commitment problem (see [START_REF] Takriti | A stochastic model for the unit commitment problem[END_REF]) where the producer has to satisfy a global power demand at every instant. The power demand, as well as other parameters such as unit breakdowns, are random. The producer is looking for strategies that minimize the overall expected production cost, over a given time horizon. Another application, which is considered in this paper, is the management of a large-scale hydro valley: here the power producer manages a cascade of dams, and maximizes the revenue obtained by selling the energy produced by turbinating the water inside the dams. Both natural inflows in water reservoirs and energy prices are random. In all these problems, the number of power units and the number of time steps are usually large (see de [START_REF] De Matos | Improving the performance of stochastic dual dynamic programming[END_REF]).

Standard resolution methods

One classical approach when dealing with stochastic dynamic optimization problems is to discretize the random inputs of the problem using a scenario tree. Such an approach has been widely studied within the stochastic programming community (see [START_REF] Heitsch | Scenario tree modelling for multistage stochastic programs[END_REF], [START_REF] Shapiro | Lectures on Stochastic Programming[END_REF], and used to model and solve energy problems, e.g. by [START_REF] Pflug | Multistage Stochastic Optimization[END_REF]. One of the advantages of such a technique is that, as soon as the scenario tree is drawn, the derived problem can be treated by classical mathematical programming techniques. Thus, a number of decomposition methodologies have been proposed (see for instance [START_REF] Rockafellar | Scenarios and policy aggregation in optimization under uncertainty[END_REF], [START_REF] Carpentier | Stochastic optimization of unit commitment: a new decomposition framework[END_REF], [START_REF] Ruszczyński | Decomposition methods in stochastic programming[END_REF], (Ruszczyński & Shapiro, 2003, Chap. 3) and applied to energy planning problems (see [START_REF] Bacaud | Bundle methods in stochastic optimal power management: A disaggregated approach using preconditioner[END_REF]). Ways to combine the discretization of the expected value together with the discretization of information in a general setting have been presented in [START_REF] Heitsch | Stability of multistage stochastic programs[END_REF], [START_REF] Pflug | Multistage Stochastic Optimization[END_REF] and [START_REF] Carpentier | Stochastic Multi-Stage Optimization volume 75 of Probability Theory and Stochastic Modelling[END_REF]). However, in a multi-stage setting, this methodology suffers from the drawback that arises with scenario trees: as it was pointed out by [START_REF] Shapiro | On complexity of multistage stochastic programs[END_REF], the number of scenarios needed to achieve a given accuracy grows exponentially with the number of time steps of the problem.

The other natural approach to solve SOC problems is to rely on the Dynamic Programming (DP) principle (see [START_REF] Bellman | Dynamic Programming[END_REF], [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF]). The core of the DP approach is the definition of a state variable that is, roughly speaking, the variable that, in conjunction with the time variable, is sufficient to take an optimal decision at every instant. It does not have the drawback of the scenario trees concerning the number of time steps since strategies are, in this context, depending on a state variable whose space dimension does not grow with time (usually linked to the number of power units in the case of power management). However, DP suffers from another drawback which is the so-called curse of dimensionality: the complexity of solving the DP equation grows exponentially with the state space dimension. Hence, solving the DP equation by brute force is generally intractable when the state space dimension goes beyond several units. In [START_REF] Vezolle | Large Scale Experiment and Optimization of a Distributed Stochastic Control Algorithm. Application to Energy Management Problems[END_REF], the authors were able to solve DP on a 10 state variables energy management problem, using parallel computation coupled with adequate data distribution, but the DP limits are around 5 state variables in a straightforward use of the method.

Another popular idea is to represent the value functions (solutions of the DP equation) as a linear combination of a priori chosen basis functions (see [START_REF] Bertsekas | Neuro-Dynamic Programming[END_REF]). This approach, called Approximate Dynamic Programming (ADP) has become very popular and the reader is referred to [START_REF] Powell | Approximate Dynamic Programming: Solving the Curses of Dimensionality[END_REF] and [START_REF] Bertsekas | Dynamic Programming and Optimal Control: Approximate Dynamic Programming[END_REF] for a precise description of ADP. This approximation drastically reduces the complexity of solving the DP equation. However, in order to be practically efficient, such an approach requires some a priori information about the problem, in order to define a well suited functional subspace. Indeed, there is no systematic means to choose the basis functions and several choices have been proposed in the literature (see [START_REF] Tsitsiklis | Feature-based methods for large-scale dynamic programming[END_REF]).

Last but not least is the popular DP-based method called Stochastic Dual Dynamic Programming (SDDP). Starting with the seminal work of [START_REF] Van Slyke | L-shaped linear programs with applications to optimal control and stochastic programming[END_REF], the SDDP method has been designed in [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF]. It has been widely used in the energy management context and lately regained interest in the Stochastic Programming community (see [START_REF] Shapiro | Analysis of Stochastic Dual Dynamic Programming Method[END_REF] and references therein). The idea is to extend Kelley's cutting plane method to the case of multi-stage stochastic problems. Alternatively it can be seen as a multistage Benders (or L-shaped) decomposition method with sampling. It consists of a succession of forward (trajectory computation) and backward (Bellman function refining) passes that ultimately aims at approaching the Bellman function as the supremum of affine hyperplanes (cuts) generated during the backward passes.

Decomposition approach

When dealing with large-scale optimization problems, the decompositioncoordination approach aims at finding a solution to the original problem by iteratively solving subproblems of smaller dimension. In the deterministic case, several types of decomposition have been proposed (e.g. by prices, by quantities or by interaction prediction) and unified in [START_REF] Cohen | Auxiliary Problem Principle and decomposition of optimization problems[END_REF] using a general framework called Auxiliary Problem Principle. In the open-loop stochastic case, i.e. when controls do not rely on any observation, it is proposed in [START_REF] Cohen | Decomposition Coordination Algorithms for Stochastic Optimization[END_REF] to take advantage of both decomposition techniques and stochastic gradient algorithms. The natural extension of these techniques to the closed-loop stochastic case (see [START_REF] Barty | A stochastic gradient type algorithm for closed-loop problems[END_REF]), i.e. when the control is a function of the available observations, fails to provide decomposed state dependent strategies. Indeed, the optimal strategy of a subproblem depends on the state of the whole system, and not only on the local state.

We recently proposed a way to use price decomposition within the closedloop stochastic case. The coupling constraints, namely the constraints preventing the problem from being naturally decomposed, are dualized using a Lagrange multiplier (price). At each iteration, the price decomposition algorithm solves each subproblem using the current price, and then uses the solutions to update it. In the stochastic context, the price is a random process whose dynamics is not available, so the subproblems do not in general fall into the Markovian setting. However, in a specific instance of this problem (see [START_REF] Strugarek | Approches variationnelles et autres contributions en optimisation stochastique[END_REF]), the author exhibited a dynamics for the optimal multiplier and showed that these dynamics were independent from the decision variables. Hence it was possible to come down to the Markovian framework and use DP to solve the subproblems. Following this idea, it is proposed in [START_REF] Barty | Decomposition of largescale stochastic optimal control problems[END_REF] to choose a parameterized dynamics for these multipliers in such a way that solving subproblems using DP becomes possible. While the approach, called Dual Approximate Dynamic Programming (DADP), showed promising results on numerical examples, it suffered from the fact that the induced restrained dual space is non-convex, leading to some numerical instabilities. Moreover, it was not possible to give convergence results for the algorithm. The method has then been improved both from the theoretical and from the practical point of view. The core idea is to replace the current Lagrange multiplier by its conditional expectation with respect to some information process, at every iteration. This information process has to be a priori chosen and adapted to the natural filtration. Moreover, if the information process is driven by a dynamic, the state in each subproblem then consists of the original state augmented by the information process, making the resolution of the subproblem tractable by DP. Interestingly, approximating the multipliers by their conditional expectations is equivalent to solving a relaxed primal problem where the almost-sure coupling constraint has been replaced by its conditional expectation with respect to the information variable, yielding a lower bound of the true optimal cost. Further, the solutions obtained by the DADP algorithm do not necessarily satisfy the initial almost-sure coupling constraint, so we must rely on a heuristic procedure to produce a feasible solution to the original problem.

Contents of the paper

The main contribution of the paper is to give a practical algorithm aiming at solving large scale stochastic optimal control problems and providing closed-loop strategies. The numerous approximations used in the algorithm, and especially the one allowing for feasible strategies, make difficult to theoretically assess the quality of the solution finally adopted. Nevertheless, numerical implementation shows that the method is promising to solve large scale stochastic optimization problems such as those encountered in the field of energy management.

The paper is organized as follows. In §2, we present the hydro valley management problem, the corresponding general SOC formulation and the DP principle. We then focus on spatial decomposition of such a problem and the difficulties of using DP at the subproblem level. In §3, we present the DADP method and give different interpretations. We then propose a way to recover an admissible solution from the DADP results and we briefly discuss the theoretical and practical questions associated to the convergence and implementation of the method. Finally, in §4, we apply the DADP method to the management of hydro valleys. Different examples, corresponding to either academic or realistic valleys, are described. A comparison of the method with SDDP is outlined.

Notations

We will use the following notations, considering a probability space (Ω, A, P):

• i, j is the set of integers between i and j;

• bold letters are used for random variables, normal font for their realizations;

• X F t (resp. X Y ) means that the random variable X is measurable with respect to the σ-algebra F t (resp. with respect to the σ-algebra generated by Y , denoted by σ(Y ));

• x generally stands for the state, u for the control, w for an exogeneous noise;

• f t stands for a dynamics, that is, a transition function modeling the system evolution along time, L t stands for a cost function at time t, K stands for a final cost function;

• V t represents a Bellman's value function at time t;

• the notation X i (resp. U i and Z i ) stands for the discrete time state process (X i 0 , . . . , X i T ) (resp. the two control processes (U i 0 , . . . , U i T -1 ) and (Z i 0 , . . . , Z i T -1 )).

Mathematical formulation

In this section, we present the modeling of a hydro valley and the associated optimization framework.

A generic formulation

We are interested in solving a multistage stochastic optimal control problem over a discrete-time horizon 0, T . In this problem we consider multiple stochastic systems indexed by i ∈ 1, N , that follow independent dynamics but that must satisfy a coupling constraint.

More precisely, we want to address the following problem min

(X i ,U i ) i∈ 1,N E N i=1 T -1 t=0 L i t (X i t , U i t , W t ) + K i (X i T ) , (1a)
s.t. X i t+1 = f i t (X i t , U i t , W t ) , X i 0 given , (1b)
U i t σ(W 0 , . . . , W t ) , (1c)
N i=1 Θ i t (X i t , U i t , W t ) = 0 . (1d)
Constraints (1b) represent the dynamics and constraints (1c) are the nonanticipativity constraints, that is, the fact that each control U i t at time t, considered as a random variable, has to be measurable with respect to the sigma-field σ(W 0 , . . . , W t ) generated by noises up to time t. The last constraints (1d) express the interactions between the production units i. They represent an additive coupling with respect to the different production units, which is termed the "spatial coupling of the problem". Such a general modeling covers other cases than the cascade problem, such that the unit commitment problem, or the problem of exchanging energy on a smart grid.

Dams management problem

We consider a hydro valley constituted of N cascaded dams as represented in Figure 1. The water turbinated at a dam produces energy which is sold on electricity markets, and then enters the nearest downstream dam. The overall goal of the decision maker is to maximize the profit obtained by selling the produced energy on a market. We consider that the hydro valley manager acts as a price follower, in the sense that energy prices are independent of energy produced by the hydro valley. Note that the valley geometry may be more complicated than a pure cascade: see for example the valleys represented at Figure 4.

The representative variables of dam i at time t are u i t for the turbinated water, x i t for the current water volume, a i t for the natural water inflow entering dam i, p i t for the market value of the water at dam i. The randomness is given by w i t = (a i t , p i t ). The modeling of a dam takes into account a possible overflow: the spilled water does not produce electricity, but enters the next downstream dam.

We now cast the problem in the generic framework presented at §2.1, with a slight abuse of notation (U i t stands for (U i t , Z i t ) here). • The dam dynamics (corresponding to Equation (1b)) reads

x 3 t Dam 1 Dam 2 Dam 3 a 1 t x 1 t u 1 t a 2 t u 2 t x 2 t a 3 t u 3 t
x i t+1 = x i t -u i t + a i t + z i t -s i t = f i t (x i t , (u i t , z i t ), w i t ) , (2a)
where s i t is the volume of water spilled by overflowing the dam:

s i t = max 0, x i t -u i t + a i t + z i t -x i . (2b)
The constant value x i stands for the maximal capacity of dam i. The outflow of dam i, that is, the sum of the turbinated water and of the spilled water, is denoted by z i+1 t :

z i+1 t = u i t + s i t = g i t (x i t , (u i t , z i t ), w i t ) . (2c)
This last equation corresponds to Equation (1d) in the general framework. Note that the dynamic equations (2a) are nonlinear because of the max operator in the definition (2b) of the spilled water volume. We assume the Hazard-Decision information structure: the control u i t applied at time t is chosen once the noise w i t at time t has been observed. It is thus possible to ensure that the dam always remains above its minimal admissible volume x i by limiting the control range:

u i ≤ u i t ≤ min u i , x i t + a i t + z i t -x i . Remark 1.
As will be seen in §4, the typical time step length we use is the month (with a time horizon of one year). It is thus reasonable to assume the Hazard-Decision framework, the control applied for a given month being in fact implemented each day taking into account the observed information on a daily basis.

• The objective function of dam i is the sum of different terms.

-The cost at each time t ∈ 0, T -1 is:

L i t (x i t , (u i t , z i t ), w i t ) = -p i t u i t + ε(u i t ) 2 .
The first linear term corresponds to the opposite of the profit when selling the energy produced by the turbinated water on the energy market. The second term ε(u i t ) 2 models the operating cost of the turbine as a quadratic term, and is usually small. This last term ensures the strong convexity of the cost function.

-The final cost at time T is:

K i x i T = α i min{0, x i -x i T } 2 .
It corresponds to a quadratic penalization around a target value x i representing the desired water volume in the dam at the end of the time horizon.

Both functions appear in the cost (1a) in the generic problem formulation.

Dynamic Programming like approaches

In the remainder of the paper, we assume that we are in the so-called white noise setting.

Assumption 1. Noises W 0 , . . . , W T -1 are independent over time.

This assumption can be alleviated, in the case where it is possible to identify a dynamics in the noise process (such as an ARMA model), and by incorporating this new dynamics in the state variables (see [START_REF] Maceira | Use of par(p) model in the stochastic dual dynamic programming optimization scheme used used in the operation planning of the brazilian hydropower system[END_REF] on this topic).

Under Assumption 1, Dynamic Programming (DP) applies to Problem (1): there is no optimality loss to seek each control U i t at time t as a function of both the state and the noise at time t. Then, the Bellman functions V t are obtained by solving the Dynamic Programming equation backwards in time:

V T (x T ) = N i=1 K i (x i T ) , ( 3a
) V t (x t ) = E min u 1 t ,...,u N t N i=1 L i t (x i t , u i t , W t ) + V t+1 f t (x t , u t , W t )) . (3b)
where

x t = (x 1 t , . . . , x N t ), u t = (u 1 t , . . . , u N t ) and f t (x t , u t , W t ) is the collection of new states f i t (x i t , u i t , W t ).
The DP equation is agnostic to whether the state and control variables are continuous or discrete, whether the constraints and the cost functions are convex or not, etc. However, in order to exhaustively solve the DP equation, we need to have discrete state, and to be able to solve each equation to optimality. In practice, the method is subject to the curse of dimensionality and cannot be used for large-scale optimization problems. For example, applying DP to dams management problems is practically untractable for more than five dams (see the results given at §4.3).

Another way to compute the Bellman functions associated to Problem (1) is to use the Stochastic Dual Dynamic Programming (SDDP) method. The method has been first described in [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF], and its convergence has been analyzed in [START_REF] Philpott | On the convergence of stochastic dual dynamic programming and related methods[END_REF] for the linear case and in [START_REF] Girardeau | On the convergence of decomposition methods for multistage stochastic convex programs[END_REF] for the general convex case. SDDP recursively constructs an approximation of each Bellman function as the supremum of a number of affine functions, thus exploiting the convexity of the Bellman functions (arising from the convexity of the cost and constraint functions). SDDP has been used for a long time for solving large-scale hydrothermal problems (see [START_REF] De Matos | Improving the performance of stochastic dual dynamic programming[END_REF] and the references therein) and allows to push the limits of DP in terms of state dimension (see the results given at §4.4).

Spatial coupling and approach by duality

A standard way to tackle large-scale optimization problems is to use Lagrange relaxation in order to split the original problem into a collection of smaller subproblems by dualizing coupling constraints. As far as Problem ( 1) is concerned, we have in mind to use DP for solving the subproblems and thus want to dualize the spatial coupling constraints (1d) in order to formulate subproblems, each incorporating a single dam. The associated Lagrangian L is accordingly

L X, U , λ = E N i=1 T -1 t=0 L i t (X i t , U i t , W t ) + K i (X i T ) + T -1 t=0 λ t • Θ i t (X i t , U i t , W t ) ,
where the multiplier λ t associated to Constraint (1d) is a random variable. From the measurability of the variables X i t , U i t and W t , we can assume without loss of optimality that the multipliers λ t are σ(W 0 , . . . , W t )-measurable random variables.

In order to be able to apply duality theory to the problem (which is mandatory for algorithmic resolution), we make the two following assumptions.

Assumption 2. A saddle point of the Lagrangian L exists.

Assumption 3. The Uzawa algorithm applies to compute a saddle-point of L (see [START_REF] Ekeland | Convex Analysis and Variational Problems volume 28 of Classics in Applied Mathematics[END_REF], Chap. VII) for a complete presentation).

Assumption 2 corresponds to a Constraint Qualification condition and ensures the existence of an optimal multiplier. Assumption 3 allows to use a (dual) gradient ascent algorithm to compute the optimal multiplier. An important question in order to be able to satisfy these two assumptions is the choice of the spaces where the various random variables of the problem are living in. Duality theory and associated algorithms have been extensively studied in the framework of Hilbert spaces (see [START_REF] Ekeland | Convex Analysis and Variational Problems volume 28 of Classics in Applied Mathematics[END_REF]), but the transition to the framework of stochastic optimal control poses difficult challenges [START_REF] Rockafellar | Integrals which are convex functionals[END_REF][START_REF] Rockafellar | Integrals which are convex functionals, ii[END_REF]), which will be briefly presented at §3.4. One way to get rid of these difficulties is to assume that the space Ω is finite, assumption also required for convergence of SDDP. 1 When using the Uzawa algorithm to compute a saddle-point of the Lagrangian, the minimization step with respect to (X i , U i ) i∈ 1,N splits in N independent subproblems each depending on a single pair (X i , U i ), and therefore allows for a dam by dam decomposition. More precisely, the k-th iteration of the Uzawa algorithm consists of the two following steps.

1. Solve Subproblem i, i ∈ 1, N , with fixed λ (k) : min X i ,U i E T -1 t=0 L i t (X i t , U i t , W t ) + λ (k) t • Θ i t (X i t , U i t , W t ) + K i (X i T ) (4a) s.t. X i t+1 = f i t (X i t , U i t , W t ) , X i 0 given (4b) U i t σ(W 0 , . . . , W t ) , ( 4c
)
whose solution is denoted U i,(k) , X i,(k) . 2. Use a gradient step to update the multipliers λ t :

λ (k+1) t = λ (k) t + ρ t N i=1 Θ i t X i,(k) t , U i,(k) t , W t . (5)
Note that even if Subproblem (4) only involves the "physical" state variable X i t and the control variable U i t , a situation which seems favorable to DP, it also involves two exogenous random processes, namely W and λ (k) . The white noise Assumption 1 applies for the first process W , but not for the second one λ (k) , so that the state of the system cannot be summarized by the physical state X i t ! Moreover if we just use the fact that λ (k) t is measurable with respect to the past noises, the state of the system must incorporate all noises prior to time t, that is, (W 0 , . . . , W t ). The state size of the subproblem increases with time. Without some additional knowledge on the process λ (k) , DP cannot be applied in a straightforward manner: something has to be compressed in order to use Dynamic Programming.

Dual Approximate Dynamic Programming

In [START_REF] Strugarek | Approches variationnelles et autres contributions en optimisation stochastique[END_REF], for a very specific instance of Problem (1), the author exhibited the dynamics of the optimal multiplier of the coupling constraint (1d). Hence it was possible to come down to the Markovian framework and to use DP to solve the subproblems (4) with an augmented space, namely the "physical" state X i t and the state associated to the multiplier's dynamics. Following this idea for a general Problem (1), [START_REF] Barty | Decomposition of largescale stochastic optimal control problems[END_REF] proposed to choose a parameterized dynamics for the multiplier: then solving the subproblems using DP 1 Recall that the aim of the present paper is mainly to present numerical results. The reader is referred to [START_REF] Leclère | Contributions aux méthodes de décomposition en optimisation stochastique[END_REF] for these difficult theoretical questions.

became possible, the parameters defining the multiplier dynamics being updated at each iteration of the Uzawa algorithm. This new approach, called Dual Approximate Dynamic Programming (DADP), has then largely improved through a series of PhD theses [START_REF] Girardeau | Résolution de grands problèmes en optimisation stochastique dynamique[END_REF], [START_REF] Alais | Risque et optimisation pour le management d'énergies[END_REF] and [START_REF] Leclère | Contributions aux méthodes de décomposition en optimisation stochastique[END_REF]) both from the theoretical and from the practical point of view. We give here a brief overview of the current DADP method.

DADP core idea and associated algorithm

In order to overcome the obstacle explained at §2.4 concerning the measurability of random variables λ (k) t , we choose a random variable Y t at each time t, each Y t being measurable with respect to the noises W 0 , . . . , W t up to time t. We call Y = Y 0 , . . . , Y T -1 the information process associated to Problem (1).

Method foundation

The core idea of DADP is to replace the multiplier λ

(k) t by its conditional expectation E[λ (k) t | Y t ]
with respect to Y t . From an intuitive point of view, the resulting optimization problem will be a good approximation of the original one if Y t is close to the random variable λ (k) t . Note that we require that the information process is not influenced by controls because introducing a dependency of the conditioning term with respect to the control would lead to very serious difficulties for optimization.

Using this core idea, we replace Subproblem (4) by: min

X i ,U i E T -1 t=0 L i t (X i t , U i t , W t ) + K i (X i T ) + E[λ (k) t | Y t ] • Θ i t (X i t , U i t , W t ) , (6a) s
.t. X i t+1 = f i t (X i t , U i t , W t ) , X i 0 given , (6b)
U i t σ(W 0 , . . . , W t ) . (6c)
According to the Doob property (Dellacherie & Meyer, 1975, Chapter 1, p. 18), the Y t -measurable random variable E[λ

(k) t | Y t ] can be represented by a mea- surable mapping µ (k) t , that is, µ (k) t (y) = E λ (k) t Y t = y , (7)
so that Subproblem (6) in fact involves the two fixed random processes W and Y . If the process Y follows a non-controlled Markovian dynamics driven by the noise process W , i.e. if there exist functions h t such that Y t+1 = h t (Y t , W t ) then (X i t , Y t ) is a valid state for the subproblem and DP applies.

DADP algorithm

Assume that the information process Y follows the dynamics Y t+1 = h t (Y t , W t ).

• The first step of the DADP algorithm at iteration k consists of solving all the subproblems (6) with λ

(k) t fixed, that is, with µ (k) t (•)
given. It is done by solving the Bellman functions associated to each subproblem i, that is,

V i,(k) T (x i , y) = K i (x) , V i,(k) t (x i , y) = E Q i,(k) t (x i , y, W t ) , where Q i,(k) t (x i , y, w t ) is the value of min u i L i t (x i , u i , w t ) + µ (k) t (y) • Θ i t (x i , u i , w t ) + V i,(k) t+1 x i t+1 , y i t+1 s.t. x i t+1 = f i t (x i , u i , w t ) , y t+1 = h t (y, w t ) .
Storing the argmin obtained during the Bellman resolution, we obtain the optimal feedback laws γ i,(k) t as functions of both the state (x i , y) and the noise w t at time t. These functions allow to compute the optimal state and control processes U i,(k) , X i,(k) of subproblem i at iteration k. Starting from X i,(k) 0 = X i 0 the optimal control and state variables are obtained by applying the optimal feedback laws from t = 0 up to T -1:

U i,(k) t = γ i,(k) t (X i,(k) t , Y t , W t ) , X i,(k) t+1 = f i t (X i,(k) t , U i,(k) t , W t ) .
• The second step of the DADP algorithm consists of updating the multiplier process λ (k) . Instead of updating the multipliers themselves by the standard gradient formula

λ (k+1) t = λ (k) t + ρ t N i=1 Θ i t X i,(k) t , U i,(k) t , W t , (8)
it is sufficient to deal with their conditional expectations w.r.t. Y t . Using the optimal processes X i,(k) and U i,(k) obtained at the previous step of the algorithm for all subproblems, the conditional deviation from the coupling constraint is represented by a measurable mapping ∆

(k) t : ∆ (k) t (y t ) = E N i=1 Θ i t X i,(k) t , U i,(k) t , W t Y t = y t . (9)
Gathering the functional representations ( 7) and ( 9)of the conditional multiplier and of the conditional deviation, the gradient update reduces to the following functional expression:

µ (k+1) t (•) = µ (k) t (•) + ρ t ∆ (k) t (•) . ( 10
)
This last equation is equivalent to the multipliers conditional expectation update:

E λ (k+1) t Y t = E λ (k) t Y t + ρ t E N i=1 Θ i t X i,(k) t , U i,(k) t , W t Y t . (11)
From a practical point of view, computing the gradients using Formula (9), instead of (8) opens the way to important numerical improvements in the DADP algorithm. Indeed, instead of a gradient formula in a large space, we can use more sophisticated direction descent algorithms: as a matter of fact, if the support of the random variable Y t is finite, it becomes possible to efficiently implement a quasi-Newton method, thus obtaining a much faster convergence than the one of the standard gradient ascent method (see §4.3.2 for details).

DADP algorithm is depicted in Figure 2.

Multiplier function µ (k) t (y) • • • Solving subproblem 1: DP on (X 1 t , Y t ) Solving subproblem N : DP on (X N t , Y t ) E N i=1 Θ i t • Y t = y ∆ (k) t (y) = 0 ? µ (k+1) t (•) = µ (k) t (•) + ρ t ∆ (k) t (•) Θ i t X i,(k) t , U i,(k) t , W t Information Process Y t+1 = h t (Y t , W t )

DADP interpretations

The DADP method, as it has been presented up to now, makes use of an approximation of the optimal multiplier, that is, the multiplier λ t is replaced by its conditional expectation E λ t Y t . Such an approximation is equivalent to a decision-rule approach for the dual problem (see also [START_REF] Kuhn | Primal and dual linear decision rules in stochastic and robust optimization[END_REF]), obtained by imposing that the dual variables λ t is measurable with respect to Y t .

DADP may also be viewed as a relaxation of the constraints in the primal problem. More precisely, we replace the almost sure coupling constraint (1d) by the following conditional expectation constraint,

E N i=1 Θ i t (X i t , U i t , W t ) Y t = 0 . ( 12
)
Proposition 1. Assume that the Lagrangian associated with this relaxed problem has a saddle point. Then the DADP algorithm on Problem (1) can be interpreted as the Uzawa algorithm applied to the relaxed Problem.

Proof. Consider the duality term E E[λ

(k) t | Y t ] • Θ i t (X i t , U i t , W t
) which appears in the cost function of subproblem i in DADP. This term can be written equivalently E λ

(k) t • E[Θ i t (X i t , U i t , W t ) | Y t ] ,
which corresponds to the dualization of the coupling constraint handled in the relaxed problem. DADP thus consists of replacing an almost-sure constraint by its conditional expectation w.r.t. the information variable Y t . From this interpretation, we deduce that the optimal value provided by DADP is a guaranteed lower bound of the optimal value of Problem (1).

Admissibility recovery

Solving the relaxed problem, that is Problem (1), where constraints (1d) is replaced by the less binding constraints (12), does not necessarily yield a solution admissible for Problem (1). Nevertheless it produces at each time t a set of N local Bellman functions V i,∞ t , each depending on the extended state (x i t , y t ). We use these functions to produce an approximation V ∞ t of the "true" Bellman function V t of the global state x 1 t , . . . , x N t by simply summing the local Bellman functions:

V ∞ t x 1 t , . . . , x N t , y t = N i=1 V i,∞ t x i t , y t .
We then obtain an admissible feedback policy for Problem (1): for any value of the state x 1 t , . . . , x N t , any value of the information y t and any value of the noise w t at time t, the control value is obtained by solving the following one-step DP problem min

(u 1 t ,...,u N t ) N i=1 L i t x i t , u i t , w i t + V ∞ t+1 x 1 t+1 , . . . , x N t+1 , y t+1 , s.t. x i t+1 = f i t x i t , u i t , w i t , i ∈ 1, N , y t+1 = h t y i t , w i t , N i=1 Θ i t (x i t , u i t , w t ) = 0 .
In this framework, DADP can be viewed as a tool allowing to compute approximated Bellman functions for Problem (1) which in turns yields an online admissible feedback policy for Problem (1). Applying this online feedback policy along a bunch of noises scenarios allows to compute a Monte Carlo approximation of the cost, which is accordingly a stochastic upper bound of the optimal value of Problem (1).

Theoretical and practical questions

The theoretical questions linked to DADP are addressed in [START_REF] Leclère | Contributions aux méthodes de décomposition en optimisation stochastique[END_REF], and the practical ones in [START_REF] Girardeau | Résolution de grands problèmes en optimisation stochastique dynamique[END_REF] and [START_REF] Alais | Risque et optimisation pour le management d'énergies[END_REF].

Theoretical questions

In the DADP approach, we treat the coupling constraints of a stochastic optimization problem by duality methods and solve it using the Uzawa algorithm. The Uzawa algorithm is a dual ascent method which is usually described in an Hilbert space such as L 2 (Ω, A, P, R n ), but we cannot guarantee the existence of an optimal multiplier in such a space. To overcome the difficulty, the approach consists of extending the setting to the non-reflexive Banach space L ∞ (Ω, A, P, R n ), to give conditions for the existence of an optimal multiplier in L 1 Ω, A, P; R n (rather than in the dual space of L ∞ ) and to study the Uzawa algorithm convergence in this space.

Practical questions

An important practical question is the choice of the information variables Y t . We present here some possibilities.

Perfect memory: Y

t = W 0 , . . . , W t .
From the measurability properties of λ

(k) t , we have E[λ (k) t | Y t ] = λ (k)
t , that is, there is no approximation! Indeed a valid state for each subproblem is X t , W 0 , . . . , W t : the state is growing with time. 2. Minimal information: Y t = 0.

Here λ

(k) t is approximated by its expectation E[λ (k) t ]
. The information variable does not deliver any online information, and a valid state for subproblem i is X i t .

Dynamic information: Y

t+1 = h t Y t , W t+1 .
This choice corresponds to a number of possibilities, as mimicking the state of another unit, or adding a hidden dynamics. A valid state for subproblem i is X i t , Y t . The question of accelerating the DADP algorithm by using a more sophisticated method than the simple gradient ascent method in the multiplier update step has been discussed at the end of §3.1.2. Numerical experiments have shown that it has a great impact on the convergence speed of the method (see §4.3.2). Another improvement would be to replace the standard Lagrangian by an augmented Lagrangian.

Numerical experiments

In this section, we present numerical results obtained on a large selection of hydro valleys. Some of these valleys (see Figure 4) correspond to academic examples, in the sense that their characteristics (size of dams, range of controls, inflows values) do not rely on existing valleys. These examples allow us to quantify the performance of different optimization methods (DP, DADP and SDDP) on problems of increasing size, from a valley incorporating 4 dams, and thus solvable by DP, up to a valley with 30 dams, and thus facing the curse of dimensionality ( §4.3 and §4.4). We also present two instances corresponding to more realistic hydro valleys, where the models respect the orders of magnitude of the dam sizes of existing valleys ( §4.5).

All the results presented here have been obtained using a 3.4GHz, 4 cores -8 threads Intel Xeon E3 based computer.

Application of DADP to a hydro valley

We go back to the problem formulation presented at §2.2. In order to implement the DADP algorithm, we dualize the coupling constraints

Z i+1 t -g i t (X i t , U i t , W i t , Z i t ) = 0 , (13)
and we denote by λ i+1 t the associated multiplier (random variable). When minimizing the dual problem at iteration k of the algorithm, the product of ( 13) with a given multiplier by λ i+1,(k) t is additive with respect to the dams, that is, the term λ i+1,(k) t

• g i t X i t , U i t , W i t , Z i t pertains to dam i subprob- lem, whereas the term λ i+1,(k) t
• Z i+1 t pertains to dam i+1 subproblem, hence leading to a dam by dam decomposition for the dual problem maximization in (X, U , Z) at λ i+1,(k) t fixed.

DADP implementation

The DADP method consists of choosing a multiplier process Y and then replacing the coupling constraints by their conditional expectations with respect to Y t . Here we adopt the choice Y t = 0 (minimal information), so that Constraints ( 13) are replaced in the approximated problem by their expectations:

Dam i x i t a i t s i t z i t u i t z i+1 t Dam i + 1
E Z i+1 t -g i t (X i t , U i t , W i t , Z i t ) = 0 . ( 14
)
The expression of Subproblem ( 6) attached to dam i reads

min U i ,Z i ,X i E T -1 t=0 L i t X i t , U i t , W i t , Z i t + E λ i,(k) t • Z i t (15a) -E λ i+1,(k) t • g i t X i t , U i t , W i t , Z i t + K i X i T , (15b)
s.t. X i t+1 = f i t (X i t , U i t , W t ) , X i 0 given (15c) U i t σ(W 0 , . . . , W t ) . (15d)
Because of the crude relaxation due to a constant Y i t , the multipliers λ i,(k) t appear only in the subproblems by means of their expectations E[λ i,(

], so that all subproblems involve a 1-dimensional state variable, that is, the dam stock X i t , and hence are easily solvable by Dynamic Programming. We denote by U i,(k) , Z i,(k) , X i,(k) the optimal solution of each subproblem i, and by V i,(k) t (x i ) the Bellman function obtained for each dam i at time t.

With the choice of constant information variables Y i t , the coordination update step (11) reduces to

E λ i,(k+1) t = E λ i,(k) t + ρ t E Z i+1,(k) t -g i t X i,(k) t , U i,(k) t , W i t , Z i,(k) t , (16)
that is, a collection of deterministic equations involving the expectation of (13) which is easily estimated by a Monte Carlo approach. Assume that DADP converges, leading to optimal Bellman functions V i,∞ t . We know that the initial almost-sure coupling constraints are not satisfied. To recover admissibility, we use the heuristic rule proposed at §3.3, solving the following deterministic one-step DP problem:

min (u 1 ,...,u N ) N i=1 L i t x i , u i , w i t , z i + V ∞ t+1 x 1 t+1 , . . . , x N t+1 , (17a)
s.t. x i t+1 = f i t x i , u i , w i t , z i ∀i , (17b)
z i+1 = g i t (x i , u i , w i t , z i ) ∀i . (17c)

Complete process

We can summarize the whole process as follows. In the optimization stage we first compute the local Bellman functions V i,∞ t , and form the approximate global Bellman functions V ∞ t by summing the local ones. In the simulation stage, we evaluate by Monte-Carlo the strategy induced by V ∞ t . We draw a large number of noise scenario, and compute the admissible control values along each scenario by solving Problem (17), from t = 0 to t = T -1, and storing payoffs.

SDDP implementation

As explained in §4.3, the controls of the original problem are discrete, which is a difficulty for SDDP implementation though recent extension has been proposed in [START_REF] Zou | Stochastic dual dynamic integer programming[END_REF]. In the optimization stage we relax the integrity constraints to obtain relaxed Bellman value functions V ∞ t . Then, in the simulation stage, we use these relaxed Bellman value functions to design policies taking into account the discrete controls by solving problems akin to Problem (17). Furthermore, we consider that the spillage is a control variable, so as to render the dynamics linear, which is the convex-costs SDDP framework.

The whole process of SDDP is as follows. In the optimization stage, lower approximations of Bellman functions V t are built iteratively. At iteration k, the procedure consists of two passes.

• During the forward pass, we sample a scenario of noise. We then simulate a stock trajectory by using the current approximation of the Bellman functions. This is done by successively solving one-step DP problem, akin to Problem (17), where V ∞ t+1 is replaced by its current piecewise linear outer-approximation, to determine the next stock value. Note that each of these one-step DP problem is a continuous quadratic programming (QP) problem.

• In the backward pass, duality theory allows to find subgradient of lower approximations of the Bellman functions. This subgradients are computed along the trajectory obtained during the forward pass, and used to construct valid cuts, that is hyperplanes that are lower than the Bellman functions. Those cuts are then added to the current outer-approximations of the Bellman functions.

In order to assess the convergence of the SDDP algorithm, we compute (say every 20 iterations of SDDP) a Monte Carlo approximation of the expected cost value with its associated 95% confidence interval, and compare the upper value of the confidence interval with the lower bound provided by SDDP up to a given threshold in order to stop the algorithm (see [START_REF] Shapiro | Analysis of Stochastic Dual Dynamic Programming Method[END_REF]). In our experiments, the Monte Carlo simulation has been made using 10,000 scenarios, and the relative convergence threshold was around 0.5%. The simulation stage is identical to the one described at §4.1.2 using the global Bellman's value function obtained by SDDP .

We have used a version of SDDP implemented in Julia (StochDynamicProgramming package2 ) built on top of the JuMP package used as a modeler (see [START_REF] Dellacherie | JuMP: A modeling language for mathematical optimization[END_REF]). The QP problems are solved using CPLEX 12.5. Every 10 iterations, redundant cuts are removed thanks to the limited memory level-1 heuristic described in [START_REF] Guigues | Dual dynamic programing with cut selection: Convergence proof and numerical experiments[END_REF]. Indeed, without cuts removal, the resolution of each QP becomes too slow as the number of cuts increases along iterations.

Results obtained for academic valleys

We model a first collection of hydro valleys including from 4 to 12 dams, with arborescent geometries (see Figure 4). The optimization problem is stated on a time horizon of one year, with a monthly time step (T = 12). All the dams have more or less the same maximal volume. The maximal amount of turbinated water for each dam varies with the location of the dam in the valley (more capacity for a downstream dam than for an upstream dam), as well as the random inflows in a dam (more inflow for an upstream dam than for a downstream dam). We assume discrete probability laws with finite supports for the inflows, and deterministic market prices. We also assume that the available turbine controls are discrete, so that each dam is in fact modeled using a discrete Markov chain. These valleys do not correspond to realistic valleys, in the sense that a true valley incorporates dams with very heterogeneous sizes.

SDDP convergence

We first illustrate the convergence of the SDDP algorithm for the 8-Dams valley on Figure 5 (note that most of the valleys display a similar convergence pattern). As explained at §4.2, the exact lower bound given by SDDP (black curve) increases along the iterations, and the gap between this lower bound and the upper value of the confidence interval (red curve) is less than 0.5% at iteration 140. corresponding to the fact that the lower bound provided by SDDP becomes greater than the lower value of the confidence interval.

DADP convergence

Let us first detail the method used for the update of the multipliers involved by DADP. Thanks to the choice of constant information variables, the gradient expression involved in the update formula ( 16) is an expectation, that can be approximated by a Monte Carlo approach. We draw a collection of statistically independent scenarios of {W t } and then compute at iteration k of DADP the optimal solutions X i,(k) t

, U i,(k) t , Z i,(k) t
of Subproblem (15) along each scenario. One has to note that this collection of scenarios is independent of the one 1. the standard Dynamic Programming method (DP), when possible, 2. the SDDP presented at §4.2, 3. the DADP method.

All these methods produce Bellman functions (optimization stage described at §4.1.2), whose quality is evaluated by the simulation stage of §4.1.2. The obtained results are given in Table 1. The lines "CPU time" correspond to the time (in minute) needed to compute the Bellman functions (optimization stage), whereas the lines "value" indicate the cost obtained by Monte Carlo on the initial model (simulation stage, performed using a 100,000 scenarios sample, except for the 12-Dams valley for which a smaller sample set was used to reduce the computational load). The comparisons between the different cost values for the same valley are thus relevant. For both SDDP and DADP, we also give the lower bound corresponding to the Bellman value obtained at the end of the optimization stage. We first note that a direct use of DP is only possible for the 4-Dams valley: it corresponds to the well-known curse of dimensionality inherent to DP. The value given by DP is the true optimal cost value for the 4-Dams valley and can be used as the reference value. The SDDP method, although relying on the integrity constraints relaxation in the optimization stage (hence a not so tight lower bound), gives excellent results for the 4-Dams valley: we thus elect SDDP as the reference method in order to evaluate the DADP method. Note that the CPU time remains reasonable, the optimization problems inside SDDP corresponding to a continuous linear-quadratic formulation (here solved using the CPLEX commercial solver).

Remark 2. Note however that all the methods we are comparing face the curse of dimensionality associated to the combinatorics of the control during the simulation stage, as the controls associated to the whole valley have to be enumerated at each time t along each scenario. This is the reason why the value obtained for the 12-Dams valley have been computed using 1, 000 scenarios (100, 000 for the others valleys) and hence are not so accurate.

We now turn to the DADP method. We first notice that the lower bound given by the method is rather bad (as a consequence of solving a problem with relaxed coupling constraints in the optimization stage), but the values obtained in the simulation stage are reasonable compared to the ones given by SDDP (as indicated by the last line of Table 1). The most noticeable point is that the CPU time needed for the optimization stage seems to grow more slowly for DADP than for SDDP. This aspect will be highlighted in §4.4.

Let us finally materialize more finely the difference in the results between SDDP and DADP. Beyond average values given in Table 1, Figure 7 represents the payoff empirical probability laws (optimal cost over the time horizon), obtained by the simulation stage using 100,000 scenarios, for both SDDP and DADP. We observe that, although the expectations are fairly close, the shapes of the two distributions differ significantly.

Challenging the curse of dimensionality

The experiments made in §4.3 seem to indicate that DADP is less sensitive to the size of the valley than the SDDP method. In order to validate this observation, we design a new collection of academic hydro valleys incorporating from 14 up to 30 dams. It is of course no longer possible to perform the simulation stage for these instances: the combinatorics induced by the set of possible values of the controls is too large to allow simulation of the valley behavior along a large set of scenarios. We thus limit ourselves to the computation of the Bellman functions (optimization stage). The corresponding results are reported in Table 2.

It appears that the CPU time required for the DADP method grows linearly with the number of dams, while the growth rate of SDDP is more or less exponential. Figure 8 shows how the CPU time varies for the three methods. As expected, DP is only implementable for small instances, say up to 5 dams. Eventually, the limits of SDDP and DADP have not really be reached, but DADP displays a near linear rate of CPU time allowing to tackle instances of even greater size.

Results for two realistic valleys

We finally model two hydro valleys corresponding to existing systems in France, namely the Vicdessos valley and the Dordogne river (see Figure 9).

The optimization problem is stated again on a one year horizon, with a monthly time step. What mainly differ here from the academic examples used at §4.3 are the characteristics of the dams. For example, the Dordogne river valley encompasses large dams (as "Bort" whose capacity is say 400) and small dams (as "Mareges" the capacity of which is equal to 35, that is, ten times smaller). This heterogeneity induces numerical difficulties, for example the requirement to have a wide range of possible controls for the small downstream reservoirs, or the need to use a fine discretization for the state grid in DP-like methods. We again assume discrete probability laws with finite support for the inflows, and we also assume that the available turbine controls are discrete.

The comparison results of SDDP and DADP are given in Table 3. As for the academic examples, SDDP displays the best results and is therefore used as the reference. The large number of possible discrete controls penalizes the DADP method, although the gap between SDDP and DADP remains limited.

Conclusion

In this article, we have depicted a method called DADP which allows to tackle large-scale stochastic optimal control problems in discrete time, such as the ones found in the field of energy management. We have presented the practical aspects of the method, without deepening in the theoretical issues arising in the foundations of the method. Lots of numerical experiments have been presented on hydro valley problems ("chained models"), which complements the ones already made on unit commitment problems ("flower models") [START_REF] Barty | Decomposition of largescale stochastic optimal control problems[END_REF]. The main conclusions are that DADP converges fast and gives near-optimal results even when using a "crude" relaxation (here a constant information process Y ). More precisely, DADP allows to deal with optimization problems that are out of the scope of standard Dynamic Programming, and beats SDDP for very large hydro valleys in terms of CPU time. We thus hope to be able to implement DADP for very large stochastic optimal control problems such as the ones encountered in smart management of urban districts, involving hundreds of houses and thus hundred of states variables. Such problems are formulated on a short-term time scale (typically a one day horizon with 15 minutes time steps), and incorporate on/off devices. In that new context, controls will have to be modeled using discrete variables (whereas this assumption was not mandatory for the study presented in this paper). Moreover, on a short-term time scale, the randomness of the markets prices plays an important role, and it will thus be necessary to take them into account as a noise in the problem.

We plan to extend this study in two directions. First to implement the DADP method for general spatial structures (not only "flower models" or "chain model", but "smart-grid models" involving a generic graph). The second direction is to implement more sophisticated decomposition methods than price decomposition. On the one hand we want to use decomposition schemes such that resource allocation or interaction prediction principle [START_REF] Cohen | Optimization by Decomposition and Coordination: A Unified Approach[END_REF]. On the other hand we want to use augmented Lagrangian based methods such as alternating direction method of multiplier (ADMM) and proximal decomposition algorithm (PDA) for decomposition in order to obtain the nice convergence properties of this kind of methods (see [START_REF] Lenoir | A survey on operator splitting and decomposition of convex programs[END_REF] for a survey).

Finally, let us mention that a theoretical work has begun in order to provide foundations of the method [START_REF] Leclère | Contributions aux méthodes de décomposition en optimisation stochastique[END_REF]). It includes conditions for existence of a multiplier in the L 1 space when the optimization problem is posed in L ∞ and conditions for convergence of the Uzawa algorithm in L ∞ . A lot of work remains to be done on these questions, mainly to relax the continuity assumption in order to be able to deal with extended functions, and to obtain more general assumptions ensuring the convergence of the Uzawa algorithm.

Figure 1 :

 1 Figure 1: Operating scheme of a hydro valley with 3 dams.

Figure 2 :

 2 Figure 2: DADP flowchart.

Figure 3 :

 3 Figure 3: Decomposition by dam.

Figure 4 :

 4 Figure 4: Some academic examples of hydro valleys.

Figure 5 :

 5 Figure 5: Convergence of SDDP for the 8-Dams valley

Figure 7 :

 7 Figure 7: 4-Dams payoff SDDP (left) -DADP (right)

Table 1 :

 1 Results obtained by DP, SDDP and DADP

	Valley	4-Dams	6-Dams	8-Dams	10-Dams	12-Dams
	DP CPU time	1600 '	∼ 10 8 '	∼ ∞	∼ ∞	∼ ∞
	DP value	-3743	N.A.	N.A.	N.A.	N.A.
	SDDP CPU time	6 '	10 '	13 '	50 '	97 '
	SDDP value	-3742	-7027	-11830	-17070	∼ -17000
	SDDP lower bound	-3754	-7050	-11960	-17260	-19490
	DADP CPU time	7 '	12 '	18 '	24 '	22 '
	DADP value	-3667	-6816	-11570	-16760	∼ -17000
	DADP lower bound	-3996	-7522	-12450	-17930	-20480
	Gap DADP/SDDP	2.0%	3.0%	2.2%	1.8%	?

Table 3 :

 3 Results obtained by SDDP and DADP

See the github link https://github.com/JuliaOpt/StochDynamicProgramming.jl.

Acknowledgments. This research benefited from the support of the FMJH "Program Gaspard Monge for optimization and operations research", and from the support from EDF.

used during the simulation stage of the complete process described at §4.1.2. We thus obtain realizations of Z i+1,(k) t

) , whose arithmetic mean gives the (approximated) gradient component at time t for the coupling between dam i and dam i + 1. This gradient can be used either in the standard steepest ascent method such as in ( 16), or in a more sophisticated algorithm such as the conjugate gradient or the quasi-Newton method. We use in our numerical experiments a solver (limited memory BFGS) of the MOD-ULOPT library from INRIA by [START_REF] Gilbert | Libopt-an environment for testing solvers on heterogeneous collections of problems[END_REF]. For all the valleys we studied, the convergence was fast (around 100 iterations regardless of the problem size). Figure 6 represents the evolution of the multipliers λ i t for the 8-Dams valley along the iterations of the algorithm. The order of magnitude of the optimal multipliers decreases with the geographical position of the link in the hydro valley. Nevertheless, the convergence rate is very similar for all links: this practical consideration remains true for almost all valleys, and it explains why the number of iterations required for the DADP convergence does not vary too much with the size of the valley.

Methods comparison

We solve Problem (1) for the first collection of academic valleys by 3 different methods: