Christophe Duhamel
email: christophe.duhamel@isima.fr

Philippe Mahey
email: philippe.mahey@isima.fr

Multicommodity flow problems with a bounded number of paths : a flow deviation approach

Keywords: Multicommodity flow, routing, flow deviation

Multicommodity flow problems with a bounded number of paths : a flow deviation approach

Christophe Duhamel, Philippe

Introduction

We will consider a network routing problem with multiple pairs of origins and destinations where we want to minimize the maximal relative congestion on the arcs of the network under the restriction that the number of paths used to carry the traffic is bounded. The problem of minimizing the flow on the most congested link is of valuable interest for the design of data communication networks (see [START_REF] Bertsekas | Data networks[END_REF]). This basic problem is known to be hard to solve even if it can be written as a linear program. When additional constraints are present in the model, like path restrictions as considered here, it will result in very difficult problems for which exact methods are likely to be useless (see [START_REF] Leblanc | Packet routing in telecommunications networks with path and flow restrictions[END_REF] for instance).

The particular path restriction we consider here is that the number of paths which support any feasible commodity flow is bounded by a given number (possibly depending on the commodity). Of course, when that number is equal to 1, we get the classical and difficult non bifurcated routing problem of unsplittable flows as denoted by Kleinberg [START_REF] Kleinberg | Approximation algorithms for disjoint path problems[END_REF], and when the number is very large, we obtain the relatively easier routing problem with any flow splitting allowed. The intermediate situation considered here is of practical interest when dealing with label switched paths (LSP) in modern broadband communication networks. Indeed, in multiple protocol labelswitched networks (MPLS) like all-optical IP backbones, different LSP are allowed to support the traffic for a given pair of nodes, but too many LSPs will deteriorate the performance of the protocols. One should thus try either to minimize the total number of LSPs (thus minimizing the number of wavelength conversions) with additional delay bounds to avoid congestion (see [START_REF] Banerjee | Wavelength-routed optical networks : Linear formulation, resource-budgeting tradeoffs and a reconfiguration study[END_REF]) , or to minimize the congestion with bounds on the number of supporting LSPs. Observe that the first choice turns to be much more intricated as it includes the problem of searching a minimal set of paths to carry a given feasible flow, which is purely combinatorial. This is why we focus here on the second choice.

We will analyze here the adaptation of the Flow Deviation algorithm for Network Routing (see [START_REF] Bertsekas | Data networks[END_REF]) to this constrained minimal congestion problem. As the Flow Deviation method was designed to solve multicommodity flow problems with separable convex costs, we will show first how the maximal congestion cost function can be approximated by these nonlinear functions, following earlier approaches based on polynomial approximation schemes and potential functions (see Bienstock's book for a comprehensive state-of-the art, [START_REF] Bienstock | Potential Function Methods for Approximately Solving Linear Programming Problems : Theory and Practice[END_REF]).

The Flow Deviation will be briefly surveyed in section 3, with emphasis on the path flow updates. In section 4, the routing model with a bounded number of paths will be established and different heuristic procedures will be proposed. These algorithmic issues will be validated in the last section.

Minimal congestion problems and potential functions

In the remainder, we will call active a path carrying a positive path flow.

Let G = (V, E) be a directed graph such that V is the set of nodes with |V | = n, and E is the set of arcs with |E| = m. Each arc is assigned a capacity C e , e ∈ E, and T is a traffic requirement matrix such that, for each pair (i, j) of nodes, t i j represents the amount of traffic required from node i to node j. Each pair of nodes such that t i j > 0 will be referred to as a commodity and the index k will be associated with a commodity, i.e. a pair of origin and destination nodes (respectively o k and d k ), and a traffic requirement t k = t o k d k . The total flow on a given arc e ∈ E will be denoted by x e = ∑ k x k e where x k e is the amount of commodity k routed on arc e. In the nodearc formulation of the multicommodity flow problem, we will need the node-arc incidence matrix A, (where a ie = +1, a je = -1, if e = (i, j) ∈ E) to express the individual flow constraints as

F k = {x k ∈ IR m | Ax k = b k , x k ≥ 0}, where b k i =    +t k if i = o k -t k if i = d k 0 otherwise
The basic problem of minimizing the most congested arc in the routing of a multicommodity flow consists in minimizing the piecewise linear convex function f (x) = max e x e c e . The problem may nevertheless be modelled as a linear program by adding an additional variable z as described below :

(MINCONG) min z s.t. ∑ k x k e -C e z ≤ 0, ∀e ∈ E x k ∈ F k , k = 1, . . ., K
Even if it is a linear multicommodity flow problem, thus an LP which can be solved by standard decomposition techniques exploiting the underlying flow structure, it is generally considered a hard problem. The first reason why this occurs is that the objective function f is convex piecewise linear and not separable with respect to arcs. The second reason is that it produces optimal solutions with a large number of active paths, being in that sense equivalent to the Maximum Concurrent Flow problem as shown below. Bienstock has reported in [START_REF] Bienstock | Potential Function Methods for Approximately Solving Linear Programming Problems : Theory and Practice[END_REF] a set of numerical experiments on very large maximum concurrent flow problems (with up to 4 • 10 5 rows and 2 • 10 6 columns) exhibiting abnormal cubic growth of the cpu time to solve them with the CPLEX dual code.

Observe that (MINCONG) is formulated without explicit capacity constraints on the total arc flows. This means that, besides its natural applications to congestion control in data networks, the routing problem (MINCONG) may be considered as an optimization formulation of the multicommodity flow feasibility problem, as we have the following relations : (MINCONG) is feasible if and only if there exists at least one path linking each origin to its destinations. Let z * be an optimal value for (MINCONG); there exists a feasible multicommodity flow if and only if z * ≤ 1.

The path structure of an optimal solution of (MINCONG) can be analyzed by considering its arc-path formulation. Let P k be the set of paths linking origin o k with destination d k and x kp be the path flow flowing on path p ∈ P k , i.e. x kp is the portion of the demand t k routed on path p. The arc-path version of (MINCONG) is then :

(PCONG) min z s.t.        ∑ k ∑ p;e∈p x kp -C e z ≤ 0 e ∈ E ∑ p∈P k x kp = t k k = 1, . . ., K x kp ≥ 0 k = 1, . . ., K, p ∈ P k
We recall that a path is active when it carries some positive flow; an arc e will be called critical when it corresponds to the maximal congestion, i.e. when x e = C e z.

An active path containing critical edges will be called a critical path. Suppose now that some commodity is routed on a critical path p at the optimal solution and that there exists a second active path p supporting that commodity which is not critical. Both paths define a cycle, so that one can modify the solution deviating a small quantity from p to p . A new basic optimal solution should be obtained when p turns to be critical. That situation can thus only occur when there are multiple optimal solutions, else :

Proposition 1 Suppose (PCONG) has a unique optimal solution; then, if any commodity is routed on a critical path in an optimal solution of (PCONG), then all paths used to route that commodity are critical.

An optimal solution to (PCONG) will contain at most K + σ active paths (and at least K), where σ is the number of critical edges.

As observed in [START_REF] Shahrokhi | The maximum concurrent flow problem[END_REF], (MINCONG) is also strongly related to the Maximum Concurrent Flow problem (MAXTHRU), where one wants to maximize the throughput of the network for a given set of capacities. The throughput Z of the network is a load factor that multiplies the traffic matrix. In the model shown below, X k represents again the k-th commodity flow :

(MAX T HRU ) max Z s.t. ∑ k X k e ≤ C e X k ∈ F k (Z) k = 1, . . . , K where F k (Z) = {X k ∈ IR m | AX k = Zb k , X k ≥ 0}.
One can easily verify that, for any optimal solution x * of (MINCONG), with optimal value z * , there exists an optimal solution X * of (MAXTHRU), with optimal value Z * , such that x * = z * X * and z * = 1/Z * . Moreover, the crucial fact is that both problem solutions share the same set of active paths.

Both problems (MINCONG) and (MAXTHRU) have received a lot of attention in the past decade, namely since the seminal paper by Shahrokhi and Matula [START_REF] Shahrokhi | The maximum concurrent flow problem[END_REF] who first proposed a fully polynomial approximation scheme to solve (MAXTHRU) with uniform capacities. They showed that the minimization of a separable exponential penalty function on the arcs yields a flow with a nearly maximal throughput. They chose the following penalty function :

φ e (x e ) = exp( 2m 2 ε x e )
where ε is a positive parameter which defines the approximation. They found a complexity bound of O(nm 7 /ε 5 ) and showed that the number of active paths is in O(m 3 /ε 2 ). Faster algorithms based on refinements of that exponential penalty function have been proposed later and extended to non uniform capacities and to other packing and covering problems (see [START_REF] Leighton | Fast approximation algorithms for multicommodity flow problems[END_REF], [START_REF] Grigoriadis | Fast approximation schemes for convex programs with many blocks and coupling constraints[END_REF], [START_REF] Garg | Faster and simpler algorithms for multicommodity flow and other fractional packing problems[END_REF]). Following Bienstock in [START_REF] Bienstock | Potential Function Methods for Approximately Solving Linear Programming Problems : Theory and Practice[END_REF], we will denote these separable penalty functions by the term of potential functions.

The flow deviation method

The Flow Deviation method (FD) is an adaptation of the classical linearization algorithm of Frank and Wolfe [START_REF] Frank | An algorithm for quadratic programming[END_REF] to solve multicommodity flow problems with convex costs. It has been first proposed by Fratta, Gerla and Kleinrock [START_REF] Fratta | The flow deviation method: an approach to store-and-forward computer-communication network design[END_REF] in the context of designing packet-switched networks and has been widely used by the transportation community (see [START_REF] Leblanc | An efficient approach to solving the road network equilibrium traffic assignment problem[END_REF]) and in telecommunications networks (see [START_REF] Bertsekas | Data networks[END_REF]). We recall below the formal ideas behind Frank-Wolfe's algorithm for general nonlinear programs with linear constraints :

Minimize Φ(x) s.t. Ax = b x ≥ 0
The method proceeds by successive linearization solving LP subproblems at each iteration t where the gradient ∇Φ(x t ) has been computed. Let xt be the optimal solution of the subproblem at iteration t

xt = Argmin x∈P ∇Φ(x t ) • x where P = {x ∈ IR n | Ax = b, x ≥ 0}
is the polyhedron of feasible solutions supposed bounded. Thus, we can assume that xt is an extreme point of P. The direction d t = xtx t is a descent direction in the sense that the directional derivative ∇Φ(x t ) • d t is strictly negative. Then the new iterate is obtained by carrying out a line search on the segment [x t , xt ] with the non linear function Φ. Convergence results have been obtained in the strictly convex case (see [START_REF] Bertsekas | Nonlinear Programming[END_REF] for example) but the method suffers from very slow convergence tail which is evidenced by the following fact : the solution of the non linear program is in general not a vertex of the feasible polyhedron so that, as the sequence of feasible solutions converges towards a point on a face of the polyhedron, the angle between successive descent directions tends to 180 • (see Fig. 1). It is shown in [START_REF] Bertsekas | Nonlinear Programming[END_REF] how typical sublinear convergence rate can be exhibited in very simple and low-dimensional situations. This drawback is somehow compensated by the observed fact that nearly optimal solutions -typically within 1% of optimalitycan be obtained very quickly. On the other hand, the positive aspect of the implementation of Frank-Wolfe's method to multicommodity flow problems is the simplicity of the subproblems resolutions which reduce to shortest-path computation and the solution update which consists in fairly deviating flows from the active paths towards the new ones without explicitly computing the individual commodity flows as we can see below. Now, to apply the Flow Deviation method to (MINCONG), we need to approximate the non smooth congestion function f (x) = max e x e /C e by a smooth convex function. This can be done by using a separable potential function as shown in the former section. An interesting link between the Flow Deviation algorithm and the resolution of (MINCONG) has been analyzed recently by Bienstock and Raskina in [START_REF] Bienstock | Asymptotic analysis of the flow deviation method for the maximum concurrent flow problem[END_REF]. They proposed an algorithm to solve (MINCONG)which alternates between a magnification step where the throughput is increased for a fixed congestion and a potential reduction step where congestion is minimized for a fixed throughput. That latter step is realized by performing flow deviation inner steps with the Kleinrock's congestion function used in [START_REF] Fratta | The flow deviation method: an approach to store-and-forward computer-communication network design[END_REF]. The interesting fact is that this function has an interpretation as a measure of the Quality of Service (QoS) of the network. It is indeed the average delay suffered by a data packet on arc e. In [START_REF] Fratta | The flow deviation method: an approach to store-and-forward computer-communication network design[END_REF], independence assumptions with Poisson arrivals on the queuing network lead to Kleinrock's delay function which is proportional to Φ e (x e ) = x e C e -x e . The total cost is simply Φ(x) = ∑ e Φ e (x e ). Observe that the objective function acts as a barrier function on the capacity constraints which can thus be ignored in the model :

x 2 x 1 x 0 x opt x *
Minimize Φ(x) s.t. x k ∈ F k , k = 1, . . ., K
where

F k = {x k ∈ IR |E| | Ax k = b k , x k ≥ 0} is the set of flows for commodity k.
Assume that we get a strictly feasible multicommodity flow at iteration k, i.e. such that x e < C e , ∀e. Then, the linearization step of (FD) reduces to separate shortest-path computations for each commodity k with arc lengths ∇Φ e (x t e ). The solution xt corresponds to the situation where all demands t k are routed separately on these new paths. Note that xt may not satisfy the capacity constraints. The new solution x t+1 is obtained by carrying out a line search on the segment [x t , xt ]. Note again that the new solution will be forced to strict feasibility by the barrier function.

The main point in the procedure above is the fact that all computations are performed on the total flow variables. However, in some situations, one may be interested in computing the optimal path flows. An arc-path formulation is then necessary and the (FD) iterations can be carried out on the path flow variables as explained in [START_REF] Bertsekas | Data networks[END_REF]. The first step of the procedure is unchanged so that, if p is the shortest path with the first-derivative arc lengths for a given commodity k, then we set xk p= t k and xkp = 0 for p = p, so that the flow deviation step can be computed in the following way :

x t+1 kp = x t kp + θ t ( xt kp -x t kp ), ∀p ∈ P k
where θ t is the optimal step size which minimizes the objective function over all θ ∈ [0, 1]. As mentioned before, the method is a descent algorithm so that the values Φ(x t ) decrease monotonically. On the other hand, we can obtain easily a lower bound of the optimal value at each iteration by computing LB t = Φ(x t ) + ∇Φ(x t ) • ( xtx t ); these succesive values need not increase monotonically, so we must keep the highest value as the current lower bound

LB t = max h=1,...,t LB h
In the original (FD) method, an equal fraction of the flow on the nonshortest paths is shifted to the shortest path. One may observe that this update of path flows will tend to increase monotonically the number of active paths. Variants of the basic Flow Deviation method can be built by modifying one of the inner steps of the algorithm, i.e.

• either by modifying the search direction, for instance by substituting the shortestpath calculations by min-cost flow subproblems for each commodity,

• or by deviating non uniform flow proportions on the newly generated paths.

To understand the first strategy, one can add the redundant constraints x k e ≤ C e , ∀e, k to the model (MINCONG); after the linearization, the direction-finding subproblem of (FD) decomposes now in K minimum-cost flow problems. The computational overhead of these subproblems compared to the original shortest-path calculations is compensated by the choice of feasible paths (considering one commodity at a time). But the drawback is that more paths are likely to be generated which is exactly what we do not wish, and the path flow updates are not so straightforward as in the original method. Our testing with that variant has confirmed these difficulties and we had to decide not to implement it in our algorithm. Observe nevertheless that some authors have chosen to implement min-cost flow computations instead of shortest-path calculations to improve the worst-case behaviour of some approximation algorithms (see [START_REF] Grigoriadis | Fast approximation schemes for convex programs with many blocks and coupling constraints[END_REF] or [START_REF] Plotkin | Fast approximation algorithms for fractional packing and covering problems[END_REF]).

There are many different strategies to modify the search direction which result in non uniform deviations from the current active paths to the shortest one. Secondorder information can be used to yield Newton or Quasi-Newton directions as suggested in [START_REF] Bertsekas | Data networks[END_REF]. Conjugate gradient strategies have also been tested in earlier works (see [START_REF] Leblanc | Improved efficiency of the frankwolfe algorithm for convex network programs[END_REF]). Again, the expected gain in convergence rate is overtaken by the extra work of performing each iteration and reconstituting the path support.

We will discuss in the last section heuristic procedures to carry out (FD) steps with a bounded number of paths.

Bounding the number of paths

We will finally consider a more complex but realistic situation where some path constraints are added to the (MINCONG) model. A way to find an interesting compromise between low congestion and unsplittable flow is to add some constraints on the number of active paths for each commodity. The general case has been scarcely considered in the literature, the research focussing mainly on the very special cases where either a single path is forced for each commodity or two disjoint paths are required for security reasons (see Kleinberg's thesis [START_REF] Kleinberg | Approximation algorithms for disjoint path problems[END_REF] for example). Another interesting type of constraints is the case of k-splittable flows where each commodity flow may be splitted uniformly among k routes [START_REF] Baier | On the k-splittable flow problem. Algorithmica[END_REF]. Besides these studies, approximation and heuristic techniques have been devoted to the reduction of the number of supporting paths with a linear or non linear objective function (see [START_REF] Burns | Path selection and bandwidth allocation in mpls networks : a non linear programming approach[END_REF], [START_REF] Martens | Flows on few paths : algorithms and lower bounds[END_REF], [START_REF] Mirrokni | A simple polynomial time framework for reduced-path decomposition in multi-path routing[END_REF]).

Let R k , k = 1, . . ., K be the maximum number of paths on which commodity k can be routed. Thus a feasible routing is a set of R k paths chosen in the set P k for each k such that ∑ R k p=1 x kp = t k , k = 1, . . . , K. It can be referred to as a (non uniform) k-splittable flow, following Skutella and others. In [START_REF] Martens | Flows on few paths : algorithms and lower bounds[END_REF], additional constraints on the size of the trunks (path capacities) are considered. We will not consider these constraints here but instead use the congestion objective function discussed in the previous sections. These problems are all known to be NP-complete as soon as R k is lower than the number of paths supporting the optimal solution of (PCONG). A mixed-integer programming formulation can be formalized to model the bound on the number of paths, introducing 0-1 variables associated with each path in (PCONG) :

(BPCONG) min z s.t.                      ∑ k ∑ p;e∈p x kp -C e z ≤ 0 e ∈ E ∑ p∈P k x kp = t k k ∈ K x kp -t k y kp ≤ 0 k ∈ K , p ∈ P k ∑ p∈P k y kp ≤ R k k ∈ K x kp ≥ 0 k ∈ K , p ∈ P k y kp ∈ {0, 1} k ∈ K , p ∈ P k
In the following, the number of active paths to support a given feasible flow will be called the flow width. If one solves a capacitated multicommodity flow problem with linear costs, the total number of active paths at the solution is at most K + m, and, supposing K = O(m), this results in an average of 2 active paths per commodity.

In practice, unfortunately, many commodities will be routed on a single path and a few of them will spread their flow on a large number of paths. Now considering an uncapacitated model with a nonlinear barrier function as above, the flow width at the optimal solution can be bounded thanks to Caratheodory's theorem by m + 1 for each commodity. This is still a too large number and, even worse, the behaviour of the (FD) procedure will tend to add more paths than necessary in the construction of the current solution. Indeed, supposing one new path is generated at each iteration, as the procedure deviates flow in equal proportion from all active paths to these new ones, the flow width will monotonically increase with the iteration count. This nasty fact could be avoided at the cost of eliminating redundant paths in the definition of the current solution as a convex combination of path flows. Typically, pivoting steps are needed to compute the smallest representation in term of path flows.

The question we address here is how to maintain a limited number R of paths in the (FD) process. The difficulties to achieve this goal are :

• Global optimality will in general be far out of reach because most situations will result in NP-hard problems;

• the arc-path model being implicit, we can only reroute the flows issued from the cancelled paths towards the already generated paths;

• feasibility issues turn to be crucial as one does not know if the network will be able to support the traffic on a smaller number of paths.

The central idea to approximately solve (BPCONG) is to adapt the Flow Deviation (FD) to the case of k-splittable flows. One can observe that the same idea has been used very early in the seminal paper by Fratta, Gerla and Kleinrock [START_REF] Fratta | The flow deviation method: an approach to store-and-forward computer-communication network design[END_REF] who applied it to the particular case of unsplittable flows (R k = 1, ∀k). The heuristic procedure simply tries to deviate all the demand on the newly generated shortest paths without violating any capacity constraint. This can be performed commodity per commodity until no improvement in the objective function is observed. It seems clear that that simple method can stop too early with a very poor feasible solution unless the traffic load is very low. Besides that, the major problem with the adaptation of (FD) to k-splittable flows is the accumulation of active paths in the iterative process. Indeed, the flow is deviated in equal proportions from all current paths to load the new one and, either all previous path flows are set to zero (θ t = 1) or all are decreased (0 < θ t < 1) and the flow width increased by one unit, the latter situation being the most likely to occur. Now observe that the general problem to find a set of paths supporting a given feasible flow has in general many solutions. Indeed, if M k is the arc-path incidence matrix associated with the set P t+1 k of paths for the k-th commodity after updating the flows at iteration t and x t+1 (resp. f t+1 ) is the path flow vector for the new active paths (resp. the arc flow vector), we look for a solution of the following linear system :

         M 1 • • • M k • • • M K e T . . . e T . . . e T          x t+1 =          f t+1 d 1 . . . d k . . . d K         
where e is a vector of 1. This system always possess at least one solution with non negative components (the one which corresponds to the classical uniform deviation and used all paths in the current representation), but we should be able to find a better solution using less paths for the same arc flow values. Unfortunately, the problem to find a minimum supporting path set for a given flow, i.e. to compute the minimal flow width, is NP-hard (see [START_REF] Vatinlen | Minimizing congestion with a bounded number of paths[END_REF]). Moreover, the indetermination in the definition of the path support is more intricate in the case of multicommodity flows as, for a given feasible multicommodity flow, there are in general an infinity of individual arc flow solutions which satisfy the multicommodity flow constraints ∑ k x k e = x e , ∀e ∈ E. For a single positive flow on a directed graph, the flow decomposition theorem (see [START_REF] Ahuja | Network Flows : Theory and Algorithms[END_REF]) guarantees there exists a support with at most n + m paths and cycles (and at most m cycles). When the only directed cycles in the network are the ones which use the return arc from the sink back to the source, it can be easily seen that the minimum number of active paths is bounded above by mn + 2 (using paths with linearly independent incidence vectors). But this number will in general be too large to be of practical interest as network protocols do not allow for a splitting in more than, say, 10 paths, even in large dense networks.

The basic ingredients of our procedure are :

Cleaning the poorest routes: the cleaning procedure is based on the idea that some routes belonging to the initial solution may not belong to the optimal solution. Then, the flow deviation procedure will decrease the amount of flow but will still maintain a positive quantity of flow on these routes. Thus a first strategy to reduce the number of paths consists in canceling the paths which support a quantity of flow which is less than a given relative threshold. The solution is then updated, either by optimizing the cost on the corresponding reduced set of paths before generating any new path, or by fairly sharing the canceled quantity of traffic on the remaining paths. Observe that the new objective function value can have increased during that phase, but in most situations, the canceled paths will not be active any more in the future solutions, including the optimal one.

Iterative loading of commodities: a second strategy, already mentioned in the case of unsplittable flows, consists in updating the flow one commodity at a time. This can be interpreted as a Gauss-Seidel like version of the Flow Deviation method, where the line-search is performed on a single commodity flow. As a consequence, the stepsize θ wil be greater and more path flows will decrease significantly with a higher probability to be canceled.

Avoiding unprofitable shortest path calculations: in the original (FD) algorithm, a shortest path is computed for each commodity at each iteration. However, computational experiments show that, most of the time, the new shortest path already belongs to the current solution set of active paths. Only few brand new improving paths are identified during the optimization (typically, less than 100 paths are stored at the end of several thousands FD iterations). Then, it may be interesting to try to look for an improving path among the active paths instead of computing a shortest path. We have modified the improving path selection in the following way : during a given number τ of iterations, the improving path is computed among the active paths. At the begining, τ ← 10. Then, at the end of the minor iterations, a shortest path is computed. If this path is also an active one, then τ ← 2τ otherwise τ ← 10.

Controlling the flow width: there are two ways to control the flow width for each commodity. The first one (external) is based on relaxing first the width constraint, then trying to satisfy it progressively. The second one (internal) keeps the width constraint during the flow deviation procedure. In the external procedure, the width constraint is restored using the same idea as in the cleaning procedure. At each iteration of the restoration, the active path with the smallest amount of flow is fairly rerouted among the remaining ones, then a flow deviation with τ ← ∞ is applied to locally optimize the flow structure. A path control is added at each step in the internal procedure : once a flow reaches its width limit, τ ← ∞, if its width drops below the limit, then τ ← 10.

Numerical experiments

The (FD) algorithm and its variants adapted to problem (BPCONG) have been coded in C and compiled through gcc 3.2 on the gnu/linux system. The numerical tests were performed on an Intel P3 800 MHz computer with 256 Mb RAM.

Test networks:

The test networks used for the numerical experimentation are of two types : 1) a realistic core network instance with 12 nodes, 88 arcs and 74 commodities; 2) a testbed of 20 random instances; those instances are built on planar graphs (obtained through Delaunay triangulation) whose size ranges from 10 to 40 nodes. For each size, 5 different sets of demands have been randomly defined. The number of commodities has been set to 2n, 4n, 6n, n(n -1)/2 and n(n -1) where n is the number of nodes.

Initialization procedure: flow deviation needs an initial solution. Thus, one has to provide a way to compute a feasible solution. This is achieved through the following heuristic, denoted by INITFD. It first consists in cutting all demands into equal-sized packets. Then each packet is routed, one at a time, on the shortest path where arc costs are the load first derivative. Arc loads are updated after each packet has been sent. Packets are randomly chosen from the remaining ones at each iteration to prevent bias and ill behaviour. However, this is not sufficient to guaranty feasibility of the solution since wrong routing decisions may be made. In practice, when the network is not close to saturation, those problems are not likely to appear.

In Table 1, results are shown for varying values of the throughput γ (a common multiplying factor for all commodities demand). INITFD, FD and FD+ rows respectively give informations for the heuristic, the classical flow deviation method and the flow deviation method using the cleaning procedure and the improved path generation. For each method and each γ, CPU time in seconds, value of the optimal solution are reported. The last seven columns give the number of commodities in the solution with the given number of active paths. The cleaning procedure (FD+ method) helps a lot reducing the number of active paths. Combined with the modified path generation, they also provide a strong reduction in CPU time. The reason is that the cleaning procedure removes paths that are not likely to be in the optimal γ meth. time (s) solution. Thus, removing those suboptimal paths avoid a lot of iterations that would otherwise be needed to reduce their amount of flow down to zero.

z * #1 #2 #3 #4 #
Table 2 shows the impact of the required final gap (column "gap") on the FD and FD+ methods. The meaning of the columns is the same as before. It can be clearly seen that the cleaning procedure and the internal path generation help keeping low CPU time as well as providing solutions with a low number of path per commodity. This also suggests the fact that FD+ as a better experimental convergency towards the optimal solution. Table 3 summarizes the behaviour of FD and FD+ for the whole set of instances. Again, FD+ compares well against FD. BLABLABLA à ajouter deviation with internal path management (FD+int) requires less CPU time than flow deviation with external path management (FD+ext). There are two reasons: first, when FD+int reaches some commodity width limit, it is quite difficult to find an alternate path to exchange with an active one. Second, removing an active path in FD+ext does a lot of perturbation. The quality of the solution for both methods is nearly the same. No accurate gap can be produced since the only valid lower bound is unconstrained flow deviation lower bound.

Finaly, the last table shows flow deviation behaviour for the unsplittable case that is, when R = 1. We compare flow deviation with internal path management (FD+int) against two heuristics. H1 (resp. H2) sorts the commodities by decresing (resp. increasing) demands and then route each one on the shortest path, where the cost are the delay first derivative. Those costs are updated after each routing. It can be seen H1 is quite close to FD for a much lesser CPU time. This shows the limits of our approach, when the width constraint is very low.

Conclusion

The Flow Deviation algorithm remains a very versatile and prolific tool to solve routing and design problems in networks since the pioneer work of M. Gerla in the seventies. Several variants and algorithmic improvements have been proposed and validated here to take in consideration some bounds on the number of active paths to support the less congested solutions. As the resulting model leads to NPhard problems for most practical situations, heuristic have been first proposed to show how the simplicity of the Flow Deviation updates is very useful to control the number of paths for each commodity. Of course, extreme situations like the unsplittable situation associated with heavy loads are unlikely to yield nice results with these heuristics as the theoretical approximation bounds are known to be quite weak for these hard problems. In these cases, branching procedures and the efficient generation of valid cuts should be added to the basic (FD) iteration to get good solutions. These studies are currently under studies with a limited computational experience on small to medium-size networks.

Figure 1 :

 1 Figure 1: Zigzagging convergence of Frank and Wolfe's method

Table 1 :

 1 cleaning procedure and improved path generation (γ influence)

	5 #6 #7

Table 4

 4 illustrates the path control procedure on the optimal solution. The first column shows the width upper bound R for every commodities. In all cases, flow

	gap meth. time (s)	z * #1 #2 #3 #4 #5 #6 #7
	-INITFD 9.64 28.014 39 16 10 5 0 2 2
	10 -1 FD	0.00 28.012 38 17 10 5 0 2 2
	FD+	0.00 28.012 38 17 10 5 2 1 1
	10 -5 FD	0.06 27.629 38 14 10 8 0 2 2
	FD+	0.03 27.629 38 15 11 6 2 2 0
	10 -3 FD	0.95 27.441 38 14 10 8 0 0 4
	FD+	0.50 27.441 38 14 18 4 0 0 0
	10 -4 FD	17.61 27.420 38 14 10 8 0 0 4
	FD+	4.73 27.418 38 22 12 2 0 0 0
	10 -5 FD	84.63 27.418 38 14 10 8 0 0 4
	FD+	5.01 27.418 38 23 11 2 0 0 0

Table 2 :

 2 cleaning procedure and internal path generation (gap influence)

Table 3 :

 3 results for the random instances

	R meth.	time (s)	z * #1 #2 #3 #4
	no FD	5.03 27.418 38 23 11 2
	3 FD+int	1.90 27.438 38 26 10 0
	FD+ext	8.53 27.418 38 29 7 0
	2 FD+int	9.47 27.478 40 34 0 0
	FD+ext	33.38 27.478 39 35 0 0
	1 FD+int	12.05 39.608 74 0 0 0
	FD+ext	15.41 39.608 74 0 0 0
		Table 4: path control

Table 5 :

 5 unsplittable flow


Philippe Mahey

Jonas Koko

Arnaud Lenoir

Decomposition methods for a spatial model for long-term energy pricing problem

come

Decomposition methods for a spatial model for long-term energy pricing problem

Introduction

In order to forecast electricity prices for one or several coupled markets, one can simulate the markets behavior. This is usually done by modeling the supply (production units) and the demand (electricity demand forecast) and assuming some market rules. For instance, we will assume in this paper a perfect market situation. From the classical micro-economic theory, this is equivalent to one lead optimizer who is responsible for minimizing the overall cost of the network.

The present model includes energy production planning on a multiperiod horizon projected into a far future (2025) and a network of production zones, interconnected by energy transactions, like in the European Market for electricity. It results in a large-scale multistage optimization problem. In order to simulate the behavior of the production network and to estimate the interzonal prices of electricity in a far future, we need to determine strategies for the management of the yearly cycle of European water dams. Such a strategy can theoretically be computed via dynamic programming where the zonal demands, water inflows and variable renewable generation are random processes. The state variable is the vector of aggregated zonal dams levels. Its size typically lies between 10 and 20 making the computation of the optimal Bellman functions hopeless. To avoid the curse of dimensionality, a relaxed version of the problem involving a simplification of the information exchanged between zones can be solved by coordination of zonal dynamic programs giving rise to zonal strategies. Indeed, it is well known that the computational cost of Dynamic Programming depends on the number of possible states at each stage (for example, for a 20-dimensional state, if we discretize each component into 50 values, we will need to consider 50 20 choices !). Thus, decomposing by zone will reduce this huge number to a tractable number of states (no more than 3 states, including the independent noise variables, in our model). On the other hand, if we decide to decompose w.r.t. scenarios, we will face the exponential growth of the scenario tree w.r.t. the number of stages (365 days in our case study). In other words, we aim at solving the dynamic but deterministic spatial model using a fast and decentralized solution strategy which will be embedded in a huge scenario-based model in a future work.

Then we have to think about approximate methods or decomposition schemes in order to solve such problems. Decomposition techniques have been widely used to cope with the inherent difficulties of multistage optimal control problems. In the present model, we are faced with the necessity to decompose the problem into decentralized zonal subproblems for which optimal solutions some efficient Dynamic Programming routine can be used (see [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]). As the aggregate model is in general convex, piecewise linear or quadratic, regularized coordination like in Monotone Operator Splitting techniques is a natural candidate to approximate the global optimum.

Monotone operator splitting methods are designed to solve monotone inclusions including the sum of two monotone operators and the algorithms aim at alternating forward (subgradient) steps and backward (proximal) steps on each operator separately, inducing at the same time decomposition techniques and Augmented Lagrangian techniques. Theoretical results and state-of-the-art of the main splitting algorithms can be found in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Lenoir | Monotone operator splitting methods and decomposition of convex programs[END_REF]. In this paper, we will compare the behavior of two classical splitting methods, namely Alternate Direction Method of Multipliers (ADMM) and Proximal Decomposition Algorithm (PDA) applied to a large-scale multistage optimal control problem.

It is well-known that ADMM and PDA are very close algorithms [START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF], but they may differ thanks to different reformulations of the coupling relations and variables between both operators. For PDA, we will decouple the incidence arcs of the production networks and, for ADMM, production costs and transmission costs are decoupled, decomposing the calculus of the former into zonal subproblems. Both algorithms have produced many studies, variants with inexact proximal calculations, multidimensional and adaptive scaling, as well as applications to diverse fields like nonlinear mechanics [START_REF] Glowinski | Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics[END_REF], stochastic multistage optimization problems (like the Progressive Hedging method of Rockafellar and Wets [START_REF] Rockafellar | Scenarios and policy aggregation in optimization under uncertainty[END_REF]), multicommodity network flows ( [START_REF] Fukushima | Application of the alternating direction method of multipliers to separable convex programming[END_REF], [START_REF] Ouorou | A survey of algorithms for convex multicommodity flow problems[END_REF]), image reconstruction [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] or classification [START_REF] Boyd | Distributed optimization and statistical learning with the alternating direction method of multipliers[END_REF].

We compare two different separable formulations, one solved by PDA (indeed similar to Spingarn's Partial Inverse algorithm, see [START_REF] Spingarn | Applications of the method of partial inverses to convex programming:decomposition[END_REF]), and the other by ADMM, both closely related to Douglas-Rachford's splitting scheme (see [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF]). After presenting the dynamic model in section 3, numerical results on realistic simulations in the last section show similarities and contrasted performance of both algorithms. It will be shown that ADMM is better on smaller and sparse networks, while PDA is better on large and dense networks.

We first begin by presenting, in Section 2, the network based optimal control problem we want to solve. Then we propose decomposition methods, in Section 3, in order to treat this problem. In Section 4, some numerical results are presented to illustrate the behavior of the decomposition methods proposed.

The multizonal optimization model

We consider a set of geographical zones Z where there is some demand for a commodity (electricity or gas). The time line is split in T time intervals indexed by t ∈ {0, . . . , T -1} and we denote d z,t the total demand for zone z ∈ Z at step t. There are several sources for supplying the demand:

1. Each zone has local production units which deliver a quantity p z,t at time step t.

The corresponding aggregated cost is c z (p z,t ).

2. Commodity can be imported from a neighboring zone z through transportation line e = (z, z ). Interconnections between zones define a directed graph G = (Z, E ⊂ Z × Z) with |Z| = n, |E| = m, vertices and edges of which are zones and transportation lines. The quantity transported through line e is denoted f e,t and the associated cost is l e,t (f e,t ).

3. Finally, some part of the production can be stored in a local storage and reused later. The level of stored commodity in zone z at the beginning of time step t is denoted by x z,t . The initial level x z,0 is given and the final level is denoted x z,T . It obeys the dynamics

x z,t+1 = x z,t -u z,t + i z,t (1)
where u z,t is the usage of the storage (positive for withdrawal assuming for the sake of clarity that no pumping is allowed in the present model). Quantity i z,t is an additional input in the storage, assumed to be known : we can think about water inflows in dams.

The network flow balance can be stated for each period t to satisfy the demand d z,t of each zone. In practice, a failure can occur if the balance equation is not satisfied. Then a load shed η z,t ≥ 0 is introduced and penalized by a high quadratic cost h z (η z,t ) = c f ail η 2 z,t . Thus, the balance between production and demand can be written by the following equations :

p z,t + u z,t + e∈z + f e,t + η z,t = d z,t + e∈z - f e,t , ∀z ∈ Z, ∀t = 0, . . . , T -1.
(2)

Thermal production

The thermal production cost is a piecewise affine and convex function of the production levels p z,t . To avoid increasing the number of unknowns, it is replaced by a quadratic approximation (using least-square regression) as shown in Figure 1. Then, in the sequel, we consider the production cost g zt as a quadratic function.

Hydroelectric production

The level of stored commodity in zone z at the beginning of time period t is denoted by x z,t . It will take the role of the state variable. The initial level x z,0 is given and the final level is denoted x z,T . It obeys the dynamics equation [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Minimal and maximal levels are fixed for each zone and each period, so that X min z,t ≤ x z,t ≤ X max z,t , ∀z, t. We neglect here the hydroelectric production cost. On the other hand, we add a cost on the final state x zT → Ψ(x zT ) to penalize the excess of water reserves at the end of the horizon :

Ψ(x z,T ) = c f inal z max{0, x z,0 -x z,T }. (3)
Thus we impose a linear cost on the reduction in storage between the start and the end of the horizon, but we do not include rewards for an increase in the reserves.

Interzonal transfer costs

For an arc e = (z, z ) ∈ E interconnecting two zones z and z , the flow transfer during period t is the variable f e,t which is bounded by 0 ≤ f e,t ≤ κ e,t . The transfer cost is linear and denoted by l et (f et ) = c inter e f et .

Dynamic quadratic model

The model consists in minimizing the sum of the production costs under offer-demand equilibrium constraints in each zone including the dynamic equations on the state and control variables.

min

(p,u,f,x,η) T -1 t=0 z∈Z [g zt (p zt ) + h z (η zt )] + e∈E l et (f et ) + z∈Z ψ z (x zT ) (4)
Equations ( 1), (2

) ∀z ∈ Z, t ∈ [0, T -1] X min zt ≤ x zt ≤ X max zt , 0 ≤ u zt ≤ U max z δ h , ∀z ∈ Z, t ∈ [0, T -1] 0 ≤ f et ≤ κ et , ∀e ∈ E, t ∈ [0, T -1]
This large-scale model will now be reformulated to focus on the interzonal coupling. To simplify the notations, we will use the 4-dimensional vector q z,t = (p z,t , u z,t , η z,t , x z,t ) in order to build the concatenated vectors q z = (q z,t , t = 1, . . . , T -1) and q t = (q z,t , z ∈ Z). Similarly, we define the vectors d z , d t , i z , i t . For the flow variables f e,t , we define f e = (f e,t , t = 1, . . . , T -1) and f t = (f e,t , e ∈ E). Finally, the whole concatenation of production and flow values on the horizon will be naturally denoted by q and f (the bold face notation is thus reserved for subvectors with one or more running indexes).

The graph G will be represented by its node-arc (n × m) incidence matrix A, so that Af t = e∈z + f e,t -e∈z -f e,t . We will also use the (n × 4n) matrix B = I n ⊗ [1 1 1 0], so that the demand satisfaction equations can be globally rewritten as :

Bq t + Af t = d t , ∀t.
Thus, q and f are coupled through the incidence matrix A at each time period t and letting

G z (q z ) = T -1 t=0 [g zt (p z,t ) + h z (η z,t )] + ψ z (x z,T ), L e (f e ) = T -1 t=0 l e,t (f e,t ),
we can reformulate the model (4) as :

min (q,f ) z∈Z G z (q z ) + e∈E L e (f e ) (5)
Bq t + Af t = d t , t = 1, . . . , T -1 (6)
x z,t+1 = x zt -u zt + i zt ∀z, t

q z,t ∈ P z,t , f e,t ∈ F e,t ∀z, t, e where

P z,t = 0 ≤ p z,t ≤ P max z,t , 0 ≤ u z,t ≤ U max z,t , X min z,t ≤ x z,t ≤ X max z,t F e,t = {0 ≤ f e,t ≤ κ et } .
Observe that the dynamic equations ( 7) are separable with respect to zonal indexes. Later, they will be kept in the zonal subproblems where they are supposed to be treated by Dynamic Programming. The coupling constraints (6) could be easily dualized to induce zonal and network subproblems, but we propose here to use ADMM as a reference for further comparisons with other splitting schemes. ADMM will be compared with PDA as described below, and both algorithms alternate resolutions of the subproblems associated with separable Augmented Lagrangian functions. The difference between both splitting resides in the treatment of the coupling constraints. While ADMM will directly consider the Augmented Lagrangian associated with the coupling constraints (6) in model ( 5), PDA will decouple the terms Af t by creating local copies of the ingoing flows (f e , e ∈ z -) as explained in the next section.

Decomposition Algorithms

We now study decomposition strategies for our dynamic quadratic programming model [START_REF] Dallagi | A stochastic zonal decomposition algorithm for network energy management problems[END_REF].

We can reorganize the columns (arcs) of the node-arc incidence matrix A by ordering the zones and the outgoing arcs in each zone

A =     1 1 1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1 -1 -1 -1 -1 1 1 1    
Observe in the above example with 4 zones that the arc ordering gives rise to diagonal blocks (arrays of three 1) denoted hereafter by A z . We define f z,t to denote the vector of outgoing flows f e,t , e ∈ z + and we introduce copies φ z,t of the ingoing flows for each zone z associated with each -1 in row z of the incidence matrix to obtain the following equivalent reformulation :

Bq z,t + A z f z,t - e∈z - φ ez,t = d z,t , ∀z, t (8)
φ ez,t = f ez ,t , for e = (z , z) ∈ z + ∩ z -, ∀z, z ∈ Z, ∀t (9)
Observe that the demand equation ( 8) is now completely decentralized with additional variables φ ez . The set of equations ( 9) represents a coupling subspace between zones in R 2m(T -1) , denoted hereafter by A. By setting

Q z (f z , φ z ) = min qz,xz t g z,t (p z,t ) + h z (η z,t ) + e∈z + l e,t (f ez,t ) + ψ z (x z,T ) s.t. Bq z,t + A z f z,t - e∈z - φ ez,t = d z,t , ∀t
x z,t+1 = x zt -u zt + i zt , ∀t q z,t ∈ P z,t , f ez,t ∈ F e,t , φ ez,t ∈ F e,t , ∀e, t

Finally, we can put problem (5) in the following separable form :

min q,f ,φ z Q z (f z , φ z ) (f , φ) ∈ A
Each Q z is a convex function defined on a polyhedron with a maximal monotone subdifferential operator, allowing the application of PDA as seen below.

Proximal decomposition algorithm (PDA)

The Proximal Decomposition algorithm is a generalized version of Spingarn's Partial Inverse method ( [START_REF] Spingarn | Applications of the method of partial inverses to convex programming:decomposition[END_REF]) analyzed by Mahey et al [START_REF] Mahey | Proximal decomposition on the graph of a maximal monotone operator[END_REF]. It is specially tailored to minimize separable convex functions on a coupling subspace.

The key ingredient is that we need a pair of primal and dual variables lying respectively in subspace A and its orthogonal subspace A ⊥ . We denote the dual variables by (u, v) so that we have the following dual relations :

X = f φ = I I α = Eα, for some α ∈ R n(n-1) (10)
W = v w = I -I β = Dβ, for some β ∈ R n(n-1) (11)
Algorithm 1 Proximal decomposition algorithm

Step 1. Zonal subproblems: Compute (q k+1 z , f

k+1/2 z , φ k+1/2 z
) as solution of min (qz,fz,φz)

Q z (q z , f z ) + 1 2λ f z -f k z -λv k z 2 + 1 2λ φ z -φ k z -λw k z 2
Step 2. Dual update:

v k+1/2 z = 1 λ (f k+1/2 z + λv k z -f k z ) w k+1/2 z = 1 λ (φ k+1/2 z + λw k z -φ k z )
Step 3. Projection step: Compute X k+1 = (f k+1 , φ k+1 ) and W k+1 = (v k+1 , w k+1 ) as follows.

f k+1 ez = φ k+1 ez = 1 2 (f k+1/2 ez + φ k+1/2 ez ) v k+1 ez = -w k+1 ez = v k+1/2 ez - 1 2 (v k+1/2 ez + w k+1/2 ez )
The method alternates proximal steps on the primal variables (the optimal production of each zone), the convex subproblem with equilibrium guaranteed by each zone independently, and projection steps on each subspace. The cost function of the zonal subproblems take the following form at iterations k min (qz,fz,φz)

Q z (q z , f z ) + 1 2λ X z -X k z -λW k z 2 which yields the intermediate solution (q k+1 z , X k+1/2 z
). The dual update are such that X k+1/2 +λW k+1/2 = X k +λW k so that (X k+1/2 , W k+1/2 ) is the projection of (X k , W k ) onto the graph of the separable maximal monotone operator

∂[Q 1 • • • Q n ].
The different steps of the proximal algorithm are presented in Algorithm 1. In Step 2, we use the same ordering of the variables X and W (f ez at the same column index as φ ez , and w ez at the same column index as v ez ).

Observe that X k ∈ A and W k ∈ A ⊥ , ∀k, so that f k is a feasible flow for graph G. On the other hand, f k+1/2 is not feasible but it will converge asymptotically.

It is well-known that numerical enhancements are necessary to take full profit of the numerical stability and linear rate of convergence of Algorithm 1, mainly :

• Inexact proximal steps;

• Adjustable scaling parameter λ k ;

• Additional relaxation parameter extending the Douglas-Rachford splitting.

Alternating direction method of multiplier (ADMM)

Alternating direction method of multiplier (also known as Uzawa block relaxation method) has been used in various domains as an operator-splitting method: nonlinear mechanics (e.g., [START_REF] Fortin | Augmented Lagrangian Methods: Application to the Numerical Solution of Boundary-Value Problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics[END_REF][START_REF] Koko | Uzawa block relaxation for the unilateral contact problem[END_REF]), image processing (e.g., [START_REF] Koko | An augmented lagrangian method for T V g +L 1 -norm minimization[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]).

Algorithm 2 Alternating direction method of multipliers

Step 1. Zonal production subproblems: Compute (q k+1 z ) as solution of min Step 2. Multiplier update

(qz∈Pz) T -1 t=0 g zt (p z,t ) + h z (η z,t ) + y k z,
y k+1 t = y k t + 1 λ (Bq k+1 t + Af k+1 t -d t ), ∀t
Back to model [START_REF] Dallagi | A stochastic zonal decomposition algorithm for network energy management problems[END_REF], we now apply directly ADMM, building the corresponding Augmented Lagrangian, using dual multipliers y t with components y z,t associated with the flow equations (6) :

L λ (q, f , y) = z∈Z G z (q z )+ e∈E L e (f e )+ t y t (Bq t + Af t -d t ) + 1 2λ Bq t + Af t -d t 2 ,
where λ is the penalty parameter. Starting with (q 0 , f 0 , y 0 ), ADMM alternates minimization steps on the Augmented Lagrangian w.r.t. primal production variables q which decompose into zonal subproblems, and to primal flow variables f which decompose into single stage static subproblems :

q k+1 = arg min q∈P L λ (q, f k , y k ) ( 12
)
f k+1 = arg min f ∈F L λ (q k+1 , f , y k ) ( 13
)
y k+1 t = y k t + 1 λ (Bq k+1 t + Af k+1 t -d t ), ∀t. ( 14
)
where the set P represents constraints [START_REF] Fukushima | Application of the alternating direction method of multipliers to separable convex programming[END_REF] for every zone z and period t with corresponding bounds, and the set F represents the bound constraints on the flow variables.

Observe that the production subproblems ( 12) decompose in n zonal subproblems which are quadratic dynamic programs with 4(T -1) variables. On the other hand, the network flow subproblems are quadratic minimizations with simple bounds which decompose in T -1 static subproblems with m variables.

The alternating direction method of multiplier algorithm for problem ( 5) is outlined in Algorithm 2.

As already observed, similar zonal subproblems and static flow subproblems are solved in Algorithms 1 and 2. The main differences lie in the fact that zonal and flow subproblems are solved in parallel in Algorithm 1 whereas they are solved sequentially in Algorithm 2 (like Jacobi's and Gauss-Seidel's algorithms in iterative methods. Additionally, projection steps on the coupling subspaces are performed in Algorithm 1 to define feasible subsequences.

Numerical experiments

We compare, in this section, the decomposition algorithms outlined in the previous sections. The computations have been carried out on a Dell Precision T3500 computer, equipped with Intel Xeon 2.67GHz processor with 12GB RAM, using Matlab version 7. We use the MATLAB function quadprog to solve the (quadratic programming) zonal problems. quadprog implements an interior-point algorithm designed for large-scale convex quadratic programs exploiting sparsity. All algorithms are stopped if the gradient of the Lagrangian function and the relative residual of the production/demand equilibrium constraint

Bq k + Af k -d ≤ ε
where ε = 10 -4 is a tolerance defined by the user. We will first compare both algorithms with the direct application of quadprog to the centralized model on a dense network.

Dense network

We consider a complete transportation graph with m = (n -1)n. We generate randomly the additional uncontrolled water input i z,t and the demand d zt as autocorrelated processes. We also generate randomly, U max z , x z,0 , x max z,t , P j zt and c j zt . We set

• interzonal transfer cost c inter e = 1, ∀e;

• failure cost c f ail = 10 3 ;

• final state cost c f inal z = 10 6 , ∀z;

• arc capacity κ et = 5, ∀t, e.

After some tests, we set λ = 10 -5 in Algorithm 1-2 (observe that update formulae for the scaling parameter have been proposed in the literature, see [START_REF] Lenoir | Monotone operator splitting methods and decomposition of convex programs[END_REF] for instance, but where not used here). Table 1 shows he comparative performance of Algorithm 1-2 with MATLAB function quadprog for different instances.

quadprog is faster but we should observe that significant speedup could have been obtained by solving subproblems in parallel in the splitting algorithms. One can notice

Sparse network

We consider a sparse transportation graph G = (Z, E) with |Z| = n and E such that each node has no more than 4 connections. We generate data randomly as in the previous subsection. For c inter e , c f ail , c f inal z , κ et and λ, we use the same values as in the previous subsection.

Table 2 shows the performance of our decomposition algorithms on a sparse transportation network. One can again notice that both algorithms are almost comparable in terms of the number of iterations. Due to the sparsity pattern of the network (no more than 4 connections for each node), Algorithm 2 is faster as noticed in Section 4.1. MAT-LAB function quadprog is again faster for these instances, but higher speedup is expected with a parallel implementation of our decomposition Algorithm 1-2.

A realistic network

The concrete energy planning problem considered here is a one-year simulation of a standard year in a far future (say in 2025). The granularity of the dynamic model is one day, thus ignoring intra-days fluctuations. The objective is to simulate import-export marginal prices between the zones (and not to solve the detailed operation plan) for a set of adapted scenarios.

We mean by 'realistic' here a network with 8 to 20 zones representing the European countries but the data used to test the decomposition algorithms do not correspond to real-life data but to randomly generated scenarios of daily variations of demand and water input for each country during one year simulations.

For illustration purpose, we consider an energy network consisting of 8 zones: Belgium (BE), Spain (ES), France (FR), Germany (GE), Italy (IT), Portugal (PT), Switzerland (SW), United Kingdom (UK). The network is depicted in Figure 2 with the arcs capacity. The hydroelectric production data are given in Table 3 while Table 4 summarize As in the previous subsections, the additional uncontrolled (water) input i z,t and the demand d z,t are randomly generated as autocorrelated processes. For example, the annual variation of energy consumption is modelled as the following sinusoid :

µ z,t = α z + β z cos( 2πt T ), t = 1, . . . , T -1
Then the zonal demand is generated recursively by independent centered random sampling (ξ z,t , t = 1, . . . , T -1) :

d z,0 = [µ z,0 + σξ z,0 ]δ h d z,t = [µ z,t + e -ω (d z,t-1 -µ z,t-1 ) + σξ z,t ]δ h
where σ, ω are adhoc parameters and δ h = 24.

Similarly, water inflows are generated with sinusoidal fluctuations, but the latter are shifted in phase to model the fact that consumption peaks appear at the beginning of the year and water inflows will be higher during spring.

We set T = 365 (one year). Consequently, the size of the original problem is 35048. In Algorithm 1, the size of zonal subproblems varies from 3286 to 6936. In Algorithm 2 the size of zonal production subproblems is 2556 while the size of flow subproblems is 7300 (20 for each time step).

Table 5 shows the performances of the algorithms with respect to the parameter λ. We can notice that Algorithm 1 is more sensitive to the choice of the parameter λ. Overall, since the network is sparse, Algorithm 2 is more efficient. • FR, GE and PT import energy;

• ES, IT and SW export energy.

Due to the final storage (high) cost, the hydroelectric generation is privileged as expected.

in Figure 10-11, two distinct scenarios are applied applied to illustrate (on a given zone, here UK) that the zones can change their status of exporting or importing.

Conclusion

We have applied two different splitting method to a dynamic long-term multizonal energy planning problem. Both methods lead to zonal subproblems which adjust the production levels of each zone in each period. In the first method, PDA, the subproblems include all ingoing and outgoing arcs of the zone by the use of the arc flow copies. In the second method, ADMM, zonal are decoupled from the network subproblems, the latter being decomposed by periods. In both cases, convergence is guaranteed by the convexity assumptions on the local cost functions and the performance of the resulting algorithms depends heavily on the choice of the scaling parameter λ.

Further study is underway to enhance the model. Indeed, the demand d z and the additional water inflows i z are random processes. To simulate the behavior of the production network and to estimate the interzonal prices of electricity in a far future, we need to define a stochastic long-term planning model where the zonal demands and the water inflows are random processes, and the cost function is the total expected cost on the whole horizon (see [START_REF] Dallagi | A stochastic zonal decomposition algorithm for network energy management problems[END_REF] for details about the stochastic model). This stochastic problem could be solved by Stochastic Dual Dynamic Programming techniques (see [START_REF] Pereira | Multistage stochastic optimization applied to energy planning[END_REF]) but it will be limited by the huge dimension of the problem even with a small number of scenarios. We have thus shown in this paper that an appropriate decomposition techniques like the Proximal Decomposition method is able to decentralize the corresponding computation by solving local dynamic programs of relatively small dimensions.

Figure 1 :

 1 Figure 1: Quadratic approximation of a piecewise linear cost function

Figure 2 :

 2 Figure 2: Energy network with arc capacity (in GW)

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 34567891011 Figure 3: Computed hydroelectric and thermic production for Belgium

  t p z,t + u z,t + η z,t +

									f k e,t -	f e,t -d z,t
									e∈z +	e∈z -
		+	1 2λ	p z,t + u z,t + η z,t +	e∈z +	f k e,t -	e∈z -	f k e,t -d z,t	2 + ψ z (x z,T )
	s.t.	x z,t+1 = x zt -u zt + i zt , ∀t
	Flow subproblem: Compute f k+1 t	as solution of
		min ft∈Ft		e	l e,t (f e,t ) + (y k t ) (p k+1 t	+ u k+1 t	+ η k+1 t	+ Af t -d t ) +
					1 2λ	p k+1 t	+ u k+1 t	+ η k+1

t + Af t -d t 2

Table 1 :

 1 Centralized GQP Vs Algorithms 1, 2 on dense networks: N number of unknowns; N c number of coupling unknowns; M c number of coupling constraints. that both algorithms are almost equivalent in terms of the number of iterations required for convergence. For small problems (n = 4, 8, 16), Algorithm 2 is up to twice as fast as Algorithm 1. For the largest problem (i.e. n = 64), Algorithm 1 stops after 137.783 seconds while Algorithm 2 dose not stop after 1000 seconds. This is due to the large number of coupling variables since Algorithm 2 has to solve a global network subproblem (Step 2.)

		Problem data		Algorithm 1		Algorithm 2		quadprog
	n	N	N c /M c	iter CPU (Sec.) iter CPU (Sec.) iter CPU (Sec.)
	4	524	240/40	2	0.662	2	0.467	11	0.053
	8	1688	1120/80	2	1.911	2	0.805	14	0.251
	16 5936	4800/160	3	5.131	2	2.417	13	0.743
	32 22112 19840/320	3	17.445	2	24.98	13	3.034
	64 85184 80640/640	5	137.783	15	> 1000	13	19.841

Table 2 :

 2 Performances of Algorithms 1-2 on a sparse transportation network: N number of unknowns; N c number of coupling unknowns; M c number of coupling constraints.

Table 3 :

 3 Hydroelectric maximal production U max

	z	, initial level (x z,0 ), minimal and maxi-

Table 4 :

 4 Thermal generation data g z (p zt ) = 1 2 a z p 2 zt + b z p zt

Table 5 :

 5

Number of iterations Versus λ for Algorithm 1-2 Figures 3-10 show the annual production for each zone. One can notice that • BE and UK are self sufficient;

Acknowledgements

The authors acknowledge the financial support of PGMO program (Programme Gaspard Monge pour l'Optimisation). The authors are very grateful to the anonymous referees for their helpful and valuable suggestions and remarks, which greatly improved the earlier version of this paper.