Wim Van Ackooij

Nicolas Lebbe

Jérôme Malick

Regularized decomposition of large scale block-structured robust optimization problems

Keywords: Large scale block-structured problems, robust optimization, cutting-plane methods, bundle methods, unit-commitment Mathematics Subject Classification (2000) MSC 90C15, MSC 90C25, MSC 49M27

We consider a general robust block-structured optimization problem, coming from applications in network and energy optimization. We propose and study an iterative cutting-plane algorithm, generic for a large class of uncertainty sets, able to tackle large-scale instances by leveraging on their specific structure. This algorithm combines known techniques (cutting-planes, proximal stabilizations, efficient heuristics, warm-started bundle methods) in an original way for better practical efficiency. We provide a theoretical analysis of the algorithm and connections to existing literature. We present numerical illustrations on real-life problems of electricity generation under uncertainty. These clearly show the advantage of the proposed regularized algorithm over classic cutting plane approaches. We therefore advocate that regularized cutting plane methods deserve more attention in robust optimization.

f i (x i ) + ψ(x) (1)
such that x i ∈ X i , for all i = 1, . . . , m, W. van Ackooij EDF R&D, OSIRIS, 7 Boulevard Gaspard Monge, 91120 Palaiseau, France E-mail: wim.van-ackooij@edf.fr N. Lebbe CEA LETI 17 Avenue des Martyrs, 38000 Grenoble, France E-mail: lebbe.nicolas@gmail.com J. Malick CNRS, LJK Bâtiment IMAG, Université de Grenoble, 38000 Grenoble, France E-mail: jerome.malick@univ-grenoble-alpes.fr

where ψ : R n → R ∪ {+∞} and f i : R ni → R are given convex mappings and X i ⊆ R ni arbitrary compact sets. The complicating function ψ has a coupling effect on the specific block-structure: whenever ψ ≡ 0, problem [START_REF] Bandi | Tractable stochastic analysis in high dimensions via robust optimization[END_REF] can in principle be solved by solving the m smaller problems

min xi∈Xi f i (x i ). (2)
Such block-structured optimization problems appear in energy optimization (see e.g., [START_REF] Tahanan | Large-scale unit commitment under uncertainty: a literature survey[END_REF] on unitcommitment problems) and in network optimization (see e.g., [START_REF] Chatzipanagiotis | An augmented lagrangian method for distributed optimization[END_REF] with applications in network utility maximization and linear network flow problems). We now consider that the coupling function ψ is further subject to uncertainty, and that we want to take optimal decisions over the whole uncertainty set. More precisely we assume that ψ depends on a parameter d lying in a given uncertainty set D ⊂ R p which leads to the robust optimization problem

such that x i ∈ X i , for all i = 1, . . . , m

where the function ψ : R n × R p → R ∪ {+∞} is convex with respect to the first variable. In section 2, we give several examples of situations where such robust block-structured optimization problems appear. Note also that (3) can have an interpretation in the terminology of stochastic optimization with ψ viewed as a full second stage (e.g., [68, eq. ( 1)]), and as a Benders formulation (e.g., [START_REF] Geoffrion | Generalized benders decomposition[END_REF]). The goal of this paper is to design a generic efficient method for solving this robust optimization problem in large scale context by leveraging on its specific block-structure.

However, the practical difficulty of solving the two problems (3) and [START_REF] Bandi | Tractable stochastic analysis in high dimensions via robust optimization[END_REF] strongly depends on the model of uncertainty D. In some particular cases, the two problems have comparable computational costs. For example, when the set D is a finite set with small cardinality, computing Ψ in (4) consists of computing ψ(x, d) finitely many times. This is the common situation in stochastic programming and does not necessarily cause modelling issues if we allow the set of scenarios D to be sufficiently representative within a large expressive set of "realizations" (typically called scenarios). It is possible indeed to select representative scenarios by quantization (see e.g., [START_REF] Pagès | Optimal quadratic quantization for numerics: the gaussian case[END_REF]), methods from scenario reduction (see e.g., [START_REF] Heitsch | A note on scenario reduction for two-stage stochastic programs[END_REF]), or adaptive partitioning schemes (see e.g., [START_REF] Song | An adaptive partition-based approach for solving two-stage stochastic programs with fixed recourse[END_REF][START_REF] Van Ackooij | An adaptive partition-based level decomposition for solving two-stage stochastic programs with fixed recourse[END_REF]). These methods are well understood for objective functions involving expectations, but not for worst case objective functions. Their applications in our robust optimization context is unclear as it would raise many theoretical questions. These questions are related to how to produce a set D ′ , close to D, but significantly smaller than D or such that (4) is cheaper to compute. This raises questions on distances between the optimal solutions (and associated optimal solution sets) in this context, which seems to be a new research area, e.g., [START_REF] Guo | Stability analysis for mathematical programs with distributionally robust chance constraint[END_REF].

For better modelling of uncertainty in our robust optimization situation, we would like to use directly richer uncertainty sets D, such as finite sets with an exponential number of scenarios (see e.g., [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF]) or elementary geometrical shapes (see e.g., [START_REF] Ben-Tal | Robust Optimization[END_REF][START_REF] Ben-Tal | Robust solutions of uncertain linear programs[END_REF][START_REF] Ghaoui | Robust solutions to least-squares problems with uncertain data[END_REF][START_REF] Ghaoui | Robust solutions to uncertain semidefinite programs[END_REF]) as

D = d ∈ R T : M (d -D) p ≤ κ , (5)
where p ≥ 1 refers to a choice of norm, M is an appropriate matrix and D the "center" of the shape. This shape includes famous models as hyper-rectangles (p = ∞ and M = I) and hyper-ellipsoids (p = 2). For these choices of uncertainty in general, solving the robust problem (3) is much more complicated than solving [START_REF] Bandi | Tractable stochastic analysis in high dimensions via robust optimization[END_REF], because computing Ψ is usually expensive (and can be NP-hard). We refer the reader to [START_REF] Ben-Tal | Robust convex optimization[END_REF][START_REF] Bertsimas | Theory and applications of robust optimization[END_REF] for general presentation and to [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF] for an extensive discussion on NP-hardness and tractability of the computation of Ψ . For a concrete construction of D in another way, we refer the reader to [START_REF] Bandi | Tractable stochastic analysis in high dimensions via robust optimization[END_REF].

For rich and complicated sets D, it may be difficult to reformulate the problem into a tractable form. In general the robust counterpart to an arbitrary convex optimization problem is intractable, see e.g., [START_REF] Ben-Tal | Robust convex optimization[END_REF] (even if some approximation results for robust convex problems with a conic structure are discussed in [START_REF] Bertsimas | Tractable approximations to robust conic optimization problems[END_REF]). However we can always trivially generate a linearization or a cutting-plane which lies below the function Ψ . Indeed, for any given arbitrary d ∈ D, convexity of ψ(•, d) yields

ψ(x, d) + g, z -x ≤ ψ(z, d) ≤ Ψ (z), for any g ∈ ∂ x ψ(x, d), (6)
so we get a cutting plane of Ψ given by a value Ψ x = ψ(x, d) and a vector g x = g. Here ∂ x ψ(x, d) denotes the sub-differential of convex analysis of the mapping x → ψ(x, d). In the terminology of nonsmooth optimization (see e.g., [START_REF] Malick | Uncontrolled inexact information within bundle methods[END_REF][START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF]), we have a so-called (uncontrolled) inexact oracle. When

d maximizes ψ(x, •) over D (that is Ψ x = ψ(x, d) = Ψ (x)
), this linearization touches the graph of the function Ψ , and the oracle is said to be exact. Calling this oracle at several points x 1 , ..., x k in the domain of Ψ , we can set up the so-called cutting-plane model:

Ψk (x) := max j=1,...,k Ψ x j + g x j , x -x j ≤ Ψ (x), (7)
lying below the function by convexity. Such models are often used in robust optimization by cutting-plane algorithms; see e.g., the recent article [START_REF] Bertsimas | Reformulation versus cutting-planes for robust optimization[END_REF] for a numerical comparison in robust optimization between the two main approaches (reformulation to a deterministic problem vs iterative cutting plane method). In our context with general uncertainty set D, cutting plane algorithms should be methods of choice. However in a direct application, we cannot exploit the block-structure of the problem. This prevents the conventional use of cutting-plane methods in the large-scale setting of this paper, because the master optimization problem can not be handled directly. This is developed in section 3.

Contributions and outline of the article

This paper presents an efficient iterative algorithm for solving large scale robust optimization problem (3), which has the three following features:

it is generic for any uncertainty set D, -it exploits the block-structure of (3), -it incorporates numerical improvement (stabilization, warm-starting...) to tackle large-scale instances.

The only assumption on ψ and D is that we can (approximately) solve ( 4). Thus cutting-plane approximations of Ψ can be built using ( 6) and [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF]. A straightforward use of an iterative cutting-plane algorithm would make the "master" problems computing the next iterates hard to solve, due to their large-scale.

The first algorithmic specificity of our approach is to leverage on the block-structure of (3) by decomposing the master problem by duality into parallel subproblems of the form (2). We make no specific structural assumptions on the compact sets X i which can be "hard to handle" (i.e., a set with nonlinearities, nonconvexities, mixed-integer variables or probabilistic constraints). The practical assumption is that solving (2) is a relatively "easy" task compared to solving (1) or (3).

The second key algorithmic feature of the algorithm is the quadratic regularization of the cutting-plane approach. Regularized cutting-plane methods have proved to be highly efficient in nonsmooth optimization [START_REF] Ben Amor | On the Choice of Explicit Stabilizing Terms in Column Generation[END_REF][START_REF] Borghetti | Lagrangian heuristics based on disaggregated bundle methods for hydrothermal unit commitment[END_REF], stochastic optimization (see e.g., [START_REF] Oliveira | Inexact bundle methods for two-stage stochastic programming[END_REF] and [START_REF] Zverovich | A computational study of a solver system for processing twostage stochastic lps with enhanced benders decomposition[END_REF]), chance constrained optimization (e.g., [START_REF] Van Ackooij | Level bundle methods for constrained convex optimization with various oracles[END_REF][START_REF] Van Ackooij | Constrained bundle methods for upper inexact oracles with application to joint chance constrained energy problems[END_REF]), energy management ( [START_REF] Feltenmark | Dual applications of proximal bundle methods, including lagrangian relaxation of nonconvex problems[END_REF]) and combinatorial optimization [START_REF] Briant | Comparison of bundle and classical column generation[END_REF]. However they have not received much attention in the context of robust optimization. In contrast, the recent article [START_REF] Bertsimas | Reformulation versus cutting-planes for robust optimization[END_REF] reports that stabilizations of cutting-plane methods "did not yield significant improvements in computation time in our experiments" on a variety of problem instances for both robust linear optimization and robust mixedinteger optimization problems using polyhedral and ellipsoidal uncertainty sets. We will show here the exact opposite: our numerical experiments on real-life robust problems in energy optimization show a clear domination of the regularized cutting-plane method.

We design here a (regularized) cutting-plane algorithm adapted to (3), inspired but different from recent developments in nonsmooth optimization (see e.g., [START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF] and [START_REF] Hare | A proximal bundle method for nonconvex functions with inexact oracles[END_REF]). The main difference is that our master problem cannot be solved directly, which has practical implications (including the use of a second layer of cutting-plane) and theoretical impacts (including the fact that the analysis is not covered by recent studies of inexact bundle methods). We then provide an original analysis of the convergence and of the quality of the obtained final solution.

The paper is structured as follows. In section 2, we discuss how the structure of problem (3) encompasses several problems already considered in the literature. We also mention our motivating example and some of its features arising from unit-commitment. The cutting-plane approach designed for (3) is presented in Section 3, together with the ingredients for practical efficiency and a theoretical study. Numerical experiments are presented in section 4, showing the advantage of regularization over the simpler variant. The conclusion summarizes this work and advocates that regularized cutting-plane algorithms deserve more attention in robust optimization. The paper ends with appendices to recall notions and terminology used in this paper of nonsmooth optimization and operations research.

Examples of robust block-structured problem

We present in this section several situations where the robust two-stage optimization problem (3) appears.

We partially focus on the robust unit-commitment problem that we use in our numerical illustrations. We also briefly list other applications in stochastic programming.

Motivating application: large scale unit-commitment

Unit-commitment is one of the basic problems in energy optimization; see the recent review [START_REF] Tahanan | Large-scale unit commitment under uncertainty: a literature survey[END_REF]. We assume that there are m units (thermal and hydro valleys) over the time interval {1, ..., T }. Each individual unit has its own cost f i : R ni → R and production constraints X i ⊂ R ni . We consider an abstract form of this problem with the decision vector x = (x 1 , ..., x m ) min

x∈R n m i=1 f i (x i ), (8)
such that x ∈ X ∩ C.

Here X is the set of technically feasible production schedules, n = m i=1 n i and C the set of system wide constraints. We consider, for instance, the power balance constraints, as the set of system-wide constraints, that is:

C = {x ∈ R n : Ax ≥ d} , (9)
where d ∈ R T is the (net) customer/system load and Ax the total amount of produced power. We remark that realistic modelling of the French electrical system gives birth to a large scale unit-commitment problem of this form, e.g., the example given in the numerical section.

When renewable generation has overall high installed capacity, uncertainty should be accounted for in the model. We choose to consider an uncertain customer load and we further assume that complete recourse decisions exist with respect to market conditions (following [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF]). Typically the commitment variables, i.e., the integer variables, are restricted to the first stage and assumed to be fixed ahead of time, unable to respond to changes in d. A second stage deals with solving an "economic dispatch problem". The latter problem deals with finding appropriate power output levels of power plants in order to meet the actual customer load. We consider a general piecewise affine penalization function ψ t : R × R → R given by

ψ t (y, d) = max i=1,...,nt {a i,t (d -y) + b i,t } . ( 10
)
for coefficients a i,t , b i,t (e.g., see Figure 1 below). This suggestion covers for instance the penalization function of [START_REF] Dubost | A primal-proximal heuristic applied to french unitcommitment problem[END_REF]Fig 2.2]. It also covers the situation considered in [START_REF] Ben-Salem | Gestion Robuste de la production électrique à horizon court-terme[END_REF] and [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF] wherein

ψ t (y, d) = max {α(d -y), β(d -y)} . (11)
With this model, we further consider the worst case outcomes, so we are interested in solving the optimization problem (3) with this penalization function

Ψ (x) = sup d∈D T t=1 ψ t (d t , (Ax) t ), (12)
where (Ax) t refers to total generation at instant t.

Other examples

Adaptive robust optimization A popular set of approaches for a block-structured problem (3) is based on 2-stage robust optimization, also called adaptive robust optimization; see for example its application to unit-commitment in [START_REF] Bertsimas | Adaptive robust optimization for the security constrained unit commitment problem[END_REF]. In this situation the mapping ψ(x, d) is of the following form:

ψ(x, d) =    min y∈R p b(d), y such that A(d)x + B(d)y ≤ h(d) y ∈ Y, (13)
where Y can be an arbitrary convex set and A(d), B(d), h(d), b(d) are vectors and matrices of appropriate dimension depending on d.

Worst case penalization When a coupling constraint of the type g(x) ≤ b is present in (1), an approach is to penalize this constraint with a mapping attributing a value ϕ(g(x) -b) to the violation of the constraint; see e.g., [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF]. The mapping ϕ should be increasing and satisfy ϕ(z) = 0 for z ≤ 0. Then when realizing that g and/or b are subject to uncertainty, this leads to the situation described above with ψ(x, d)

:= ϕ(g d (x) -b(d))
. This is also a special case of simple recourse since

ψ(x, d) = min z {ϕ(z) s.t. g d (x) -b(d) ≤ z} . ( 14
)
Semi-infinite, Soyster like problems Problem (1) can involve constraints of the type

A(d)x ≤ b(d) for all d ∈ D, (15)
as in [START_REF] Bertsimas | Reformulation versus cutting-planes for robust optimization[END_REF]. This case fits into the general framework by defining

ψ(x, d) = min y≤0 { 0, y s.t. A(d)x + y ≤ b(d)} .
The value function is then used to define feasibility cuts. Whenever x / ∈ dom Ψ (x), one can typically also obtain some first order information: a feasibility cut. Then by solving the auxiliary problems ,d) is a convex mapping, and an optimal dual multiplier to the equation

ψ f (x, d) := {min y y : s.t. A(d)x + y ≤ b(d)}, it becomes clear that dom Ψ = d∈D x ∈ R n : ψ f (x, d) ≤ 0 . Moreover ψ f (.
A(d)x + y ≤ b(d), satisfies -A(d) T λ ∈ ∂ψ f (x, d).
From this discussion it becomes clear that feasibility cuts are but another cutting-plane model for taking into account a convex constraint of the type

Ψ f (x) ≤ 0, where Ψ f (x) = sup d∈D ψ f (x, d).
Finite small cardinality D This situation considered in [START_REF] Van Ackooij | Decomposition algorithm for large-scale two-stage unit-commitment[END_REF] and [START_REF] Van Ackooij | A comparison of four approaches from stochastic programming for large-scale unit-commitment[END_REF] has a mapping that itself is subject to complicating non-convex and/or combinatorial constraints. In this situation ψ(x, d) can be written as:

ψ(x, d) =        min (y1,...,ym) m i=1 g i (y i ) such that m i=1 θ i (y i ) ≤ b(d) y i ∈ Y i P y = P x, (16)
Although ψ(x, d) is not convex in x because of the potentially non-convex structure of some of the Y i , an implicitly computable convexification is available using the Lagrangian dual of ( 16) with respect to the coupling constraints m i=1 θ i (y i ) ≤ b(d) and P x = P y. An extensive discussion of the convexification effect of such a scheme as well as an efficient algorithm for dealing with these problems is discussed in a similar case in [START_REF] Van Ackooij | Decomposition algorithm for large-scale two-stage unit-commitment[END_REF].

(Regularized) cutting-plane methods

In this section we present two cutting-plane algorithms for solving the block-structured robust problems (3). We start in section 3.1 with recalling the cutting-plane approximation in general and how to deal with them by duality for our block-structured problems. In section 3.2 we discuss an elementary cuttingplane method for solving (3) as it helps to show the ideas and to formalize the notation. Then we present in section 3.3 the regularized cutting-plane method, that we use for the numerical illustrations of the next section. We use the classical terminology in convex optimization and operation research (see the foundational references [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF][START_REF] Geoffrion | Generalized benders decomposition[END_REF][START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I[END_REF][START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms II[END_REF], and the discussion in appendices A and B if necessary).

Cutting-plane models

Cutting-plane methods are one of the major approaches in robust optimization, and more generally in operation research. As discussed in the introduction, we have, in our situation, an oracle producing, for a given x:

either an optimality cut: a lower linearization of Ψ (if x ∈ dom Ψ ), -or a feasibility cut: a cutting-plane separating x and dom Ψ (if x ∈ dom Ψ ).

After calling the oracle on several iterates x 1 , ..., x k , we denote with O k the set of iterates at which an optimality cut was derived and F k the set of iterates at which a feasibility cut was obtained. Then we define the optimality cutting-plane model Ψk as [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] with respect to the set of iterates O k only; and similarly we introduce the feasibility cutting-plane model Ψ f k defined with respect to F k only. The idea of the cutting-plane method is to use the above-defined cutting-plane models to replace the unknown Ψ and dom Ψ in problem [START_REF] Ben-Salem | Gestion Robuste de la production électrique à horizon court-terme[END_REF]. More specifically, we replace the original convex objective Ψ by the convex piecewise affine function Ψk and the convex implicit constraint x ∈ dom Ψ by the convex polyhedral constraint Ψ f k (•) ≤ 0. Thus we consider the following approximation of problem ( 3)

min (x,r)∈R n+1 m i=1 f i (x i ) + r (17) such that x ∈ X, Ψk (x) ≤ r, Ψ f k (x) ≤ 0,
which is again equivalent to:

min (x,r)∈R n+1 m i=1 f i (x i ) + r such that x ∈ X, Ψ (x ℓ ) + g ℓ , x -x ℓ ≤ r, ℓ ∈ O k Ψ f k (x ℓ ) + g f ℓ , x -x ℓ ≤ 0, ℓ ∈ F k .
This problem, called master problem at iteration k, is still a difficult optimization problem because of the large scale and potential arbitrarily bad structure of the compact sets X i . To leverage on the block-structure of this problem, our idea is to move to the Lagrangian dual with respect to the coupling constraints. Note that from a dual view, the distinction between feasibility and optimality cuts becomes anecdotic. For all k ≥ 1,

|O k | + |F k | = k, let us define: ( Ψ (x ℓ ), ḡℓ ) := (Ψ (x ℓ ), g ℓ ) if ℓ ∈ O k (Ψ f (x ℓ ), g f ℓ ) if ℓ ∈ F k (18)
Lemma 1 (Oracle for the k-th dual function) The concave dual function Θ k associated with [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF] can be expressed, for all µ ∈ R k + such that ℓ∈O k µ ℓ = 1, as

Θ k (µ) := min (x1,...,xm)∈ m i=1 Xi m i=1 f i (x i ) + k ℓ=1 µ ℓ ḡℓ , x + k ℓ=1 µ ℓ Ψ (x ℓ ) -ḡℓ , x ℓ . ( 19
)
Thus solving the m subproblems

min xi∈Xi f i (x i ) + k ℓ=1 µ ℓ ḡℓ i , x i
gives the value Θ k (µ), a sub-gradient

Ψ (x 1 ) + ḡ1 , x(µ) -x 1 , • • • , Ψ (x k ) + ḡk , x(µ) -x k T ∈ ∂Θ k (µ)
and an associated primal solution x(µ) = (x 1 (µ), . . . , x m (µ)) ∈ X.

Proof Starting from [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF] we introduce a Lagrangian multiplier µ ∈ R k + and define the dual function by

Θ k (µ) := min xi∈Xi r∈R m i=1 f i (x i ) + r + k ℓ=1 µ ℓ Ψ (x ℓ ) + ḡℓ , x -x ℓ - ℓ∈O k µ ℓ r.
As soon as we choose r = 0 and any x i , we can pick µ = (t, 0, . . . , 0) with |t| → +∞ which leads to Θ k (µ) → -∞. Thereby we see that we can restrict µ to the subset of r on which we have r = ℓ∈O k µ ℓ r, which is exactly the set {µ : ℓ∈O k µ ℓ = 1}. Using this restriction the variable r vanishes and we find equation [START_REF] Oliveira | A bundle method for exploiting additive structure in difficult optimization problems[END_REF] for the dual function. Finally, applying a standard result of Lagrangian duality (see e.g., [33, Prop. XII.2.2.2]), the remaining assertions follow.

⊓ ⊔

This lemma shows that we can leverage on the block-structure of the problem to set up an efficient oracle of the dual function associated with [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF] by decomposing the computation into m smaller subproblems.

To maximize this concave dual function, the so-called bundle methods (see e.g., [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms II[END_REF][START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF]) are then the methods of choice. We will use these methods as "inner engines" within the algorithms presented in the next two sections. For the convenience of non-expert readers, we make a short presentation of proximal bundle methods in appendix A.

Cutting-plane method

We propose to solve (3) with a standard cutting-plane algorithm where the novelty is to use a duality to exploit the block-structure of the problem. From the developments of the previous section, this approach gives the following algorithm (inspired from the cutting-plane method of [START_REF] Van Ackooij | Decomposition algorithm for large-scale two-stage unit-commitment[END_REF] proposed in another context).

-(Step 0 : Initialization) : Generate several initial points in X to build an initial model for Ψ . Set k and the stopping tolerance δ stop .

-(Step 1 : Lagrangian dual) : Use a bundle method to maximize the dual function [START_REF] Oliveira | A bundle method for exploiting additive structure in difficult optimization problems[END_REF] with the additional constraint ℓ∈O k µ ℓ = 1. Get an approximate dual solution µ k+1 .

-(Step 2 : Primal Recovery) : Exploit the dual information generated in Step 1 to recover a solution of good quality x k+1 ∈ X. Estimate the duality gap

∆ A k := f (x k+1 ) + Ψk (x k+1 ) -Θ k (µ k+1
Ψ (x k+1 ) -Ψk (x k+1 ) ≤ δ stop . (20)
Otherwise let k = k + 1 and return in Step 1.

This algorithm uses a standard cutting-plane machinery with two specific working horses: the bundle method of Step 1 and the primal recovery of Step 2. We discuss these key computational ingredients later in section 3.4 and study here the convergence of the overall algorithm. We emphasize that the accuracy of the final iterate, though estimated, is not controlled explicitly by the algorithm due to the potential non-convexities in X. The meaning of the stopping test is thus the following.

Lemma 2 (Approximate optimality) If at iteration k the cutting-plane algorithm terminates then x k+1 is an (∆ A k + δ stop )-optimal solution to problem (3).

Proof Letting x * ∈ X 1 be the optimal solution to problem (3), we derive

f (x * ) + Ψ (x * ) ≥ f (x * ) + Ψk (x * ) ≥ f (x k+1 ) + Ψk (x k+1 ) -∆ A k ≥ f (x k+1 ) + Ψ (x k+1 ) -∆ A k -δ stop ,
when using convexity of Ψ , ∆ A k optimality of x k+1 for [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF], and the stopping criteria respectively.

⊓ ⊔

We then have the following convergence result under two assumptions on finite feasibility iterations and boundedness of subgradients. These assumptions are due to the fact that Ψ might not be defined everywhere and they are standard in similar contexts (for example, they are already present in the seminal work [START_REF] Geoffrion | Generalized benders decomposition[END_REF] on Benders decomposition).

Proposition 1 (Convergence of the cutting plane method) Assume that the algorithm does not generate infinitely many consecutive feasibility cuts, that all the used subgradients of Ψ are uniformly bounded. Let δ stop > 0, then the algorithm terminates for some finite k and x k+1 is a (∆ A k + δ stop )optimal solution to problem [START_REF] Ben-Salem | Gestion Robuste de la production électrique à horizon court-terme[END_REF].

Proof Let us establish convergence by contradiction; we assume that there is a δ stop > 0 such that for all k either x k+1 ∈ dom Ψ or

Ψ (x k+1 ) -Ψk (x k+1 ) > δ stop .
Let us denote by K the set of indexes of the iterations where the above inequality is satisfied. This set is infinite by assumption. For any k ≥ 1 and ℓ ≤ k, we have by [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF],

δ stop < Ψ (x k+1 ) -Ψk (x k+1 ) ≤ Ψ (x k+1 ) -Ψ (x ℓ ) -g ℓ , x k+1 -x ℓ .
Using the Cauchy-Schwarz inequality yields:

δ stop < Ψ (x k+1 ) -Ψ (x ℓ ) -g ℓ , x k+1 -x ℓ < Ψ (x k+1 ) -Ψ (x ℓ ) + g ℓ x k+1 -x ℓ . (21)
By assumption, all the g ℓ are bounded by a certain C ≥ 0 and by compactness of X the sequence (x k ) k admits a cluster point. Therefore, using also the fact that Ψ is continuous, the right hand side in [START_REF] Oliveira | Inexact bundle methods for two-stage stochastic programming[END_REF] tends to 0 along an appropriate subsequence, yet this immediately leads to a contradiction. Lemma 2 allows us to conclude.

⊓ ⊔

In theory the boundedness of the subgradients can be ensured on the interior of the domain of Ψ , but not on its boundary. In practice this is not a problem, since computer codes naturally require "bounded" subgradients in order to ensure well-posedness of the linear systems involved in the solution process.

When dom Ψ = R n as in our application, the convergence is guaranteed, as formalized in the following easy corollary.

Corollary 1 When dom Ψ = R n , the assumptions of Proposition 1 hold.

Proof The set X is compact, and therefore so is Co(X) by [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms I[END_REF]III.1.4.3], where Co(X) denotes the convex hull. Moreover, Ψ is convex on Co(X) so in particular locally Lipschitz continuous at each x ∈ Co(X). By compactness, Ψ is Lipschitz continuous on all Co(X) and using [16, proposition 2.1.2] there is a constant

C > 0 such that g ℓ ≤ C for all ℓ ∈ {1, . . . , L} because g ℓ ∈ ∂Ψ (x ℓ ) with x ℓ ∈ Co(X). ⊓ ⊔
We finish this section by two remarks. First the above convergence results provides the information that only the final precision matters: in early iterations one could therefore solve [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF] with a very imprecise method, as long as, when nearing convergence, problem ( 17) is solved with more precision. Second, although our interest originates from problems wherein (3) is large scale, the above algorithm can also be employed for smaller scale problems. For some problems, it might be possible to short-cut the bundle methods of Step 1 and solve directly [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF]. Proposition 1 and the preceding remark on imprecise solution still hold in this case.

Regularized cutting plane method

The cutting plane method in practice shows the well known oscillation effect and slow convergence, as illustrated in the numerical experiments of the next section (see Figure 3). To improve efficiency of the cutting plane method, we use the standard technique for adding proximal stabilization (see e.g., the textbook [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms II[END_REF]), that usually significantly reduces the overall computational effort.

The regularized cutting-plane method uses the same approach and the same numerical components as for the previous algorithm. The main difference is that we add a "proximal" quadratic regularization term to the master problem. Instead of solving [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF], we solve min

(x1,...,xm),r m i=1 f i (x i ) + r + 1 2t k x -xk 2 2 (22)
such that x ∈ X, Ψk (x) ≤ r, Ψ f k (x) ≤ 0, where xk ∈ X is a current 'best' feasible iterate. Notice that other types of quadratic regularizations of cutting-plane algorithms could be considered, including stabilization by levels or trust-regions; see e.g., [START_REF] Lemaréchal | New variants of bundle methods[END_REF]. Proximal stabilizations that we use here have not been considered yet in a context of the Benders decomposition, contrary to level and trust-region stabilizations (see [START_REF] Van Ackooij | Inexact stabilized Benders' decomposition approaches: with application to chance-constrained problems with finite support[END_REF][START_REF] Zaourar | Optimisation convexe non-différentiable et méthodes de décomposition en recherche opérationnelle[END_REF]). However, the proximal stabilization is the most popular in operation research in general [START_REF] Briant | Comparison of bundle and classical column generation[END_REF]. As in the previous section, an efficient version of a (proximal) regularized cutting-plane algorithm has to reveal the block-structure of the master problem [START_REF] Van Dinter | The unit commitment model with concave emissions costs: a hybrid benders decomposition with nonconvex master problems[END_REF]. Despite the quadratic term, the dual of problem [START_REF] Van Dinter | The unit commitment model with concave emissions costs: a hybrid benders decomposition with nonconvex master problems[END_REF] still preserves its decomposability, as formalized in the next lemma. The proof of the lemma follows the exact same lines as the one of Lemma 1 and is therefore omitted.

Lemma 3 (Oracle for the

k-th dual function) The dual function Θ R k : [0, 1] k → R is defined for all µ ∈ R k such that ℓ∈O k µ ℓ = 1 by Θ R k (µ) := min x∈X f (x) + 1 2t k x -xk 2 2 + k ℓ=1 µ ℓ ḡℓ , x + k ℓ=1 µ ℓ ( Ψ (x ℓ ) -ḡℓ , x ℓ ) . ( 23
)
Thus we have to solve the m subproblems

min xi∈Xi f i (x i ) + k ℓ=1 µ ℓ ḡℓ i - 1 t k xk i , x i + x i , 1 2t k I x i ,
which gives:

Ψ (x 1 ) + ḡ1 1 , (x(µ) -x 1 ) , • • • , Ψ (x k ) + ḡk k , x(µ) -x k T ∈ ∂Θ R k (µ)
This gives the following regularized cutting plane method:

-(Step 0 : Initialization) : Generate one or several initial elements x ∈ X and use these to build an initial model for Ψ . Set k = 1, the stopping tolerance δ stop , t 1 > 0, and let x1 = x 1 be one of these initial iterates.

-(Step 1 : Lagrangian dual) Use a bundle method to maximize [START_REF] Dubost | A primal-proximal heuristic applied to french unitcommitment problem[END_REF] with the additional constraint

ℓ∈O k µ ℓ = 1. -(
Step 2 : Primal Recovery) : Exploit the dual information generated in Step 1 to recover a solution of good quality x k+1 ∈ X. Estimate the duality gap

∆ A k := f (x k+1 ) + Ψk (x k+1 ) + 1 2t k x k+1 -xk 2 2 - Θ R k (µ * ). Compute the predicted decrease δ k := f (x k ) + Ψ (x k ) -f (x k+1 ) -Ψk (x k+1
). -(Step 3 : Oracle call) : Call the oracle for Ψ to compute a value and subgradient at x k+1 , enriching the cutting-plane model by an optimality or a feasibility cut.

-(Step 4 : Descent condition) Check if the observed decrease is at least a fraction m > 0 of predicted decrease δ k :

f (x k+1 ) + Ψ (x k+1 ) ≤ f (x k ) + Ψ (x k ) -mδ k . (24)
If ( 24) holds, let xk+1 = x k+1 , else xk+1 = xk

-(Step 5 : Stopping test) : Check if δ k ≤ δ stop (25)
and stop the algorithm if this holds. Let k = k + 1, choose t k+1 arbitrarily and return in Step 1.

This algorithm does not fit in the recent development of inexact nonsmooth optimization methods (see e.g., [START_REF] Oliveira | Level bundle methods for oracles with on demand accuracy[END_REF] and references therein) where the inexactness is the one of the oracle. Here the main inexactness comes from the fact that the approximating subproblem ( 17) is solved inexactly. The convergence of the algorithm is thus not covered by existing results. However, we can adapt standard rationale to get the following convergence result, under the same assumptions as the ones of Proposition 1.

Proposition 2 (Convergence of regularized cutting plane method) Assume that the algorithm does not generate infinitely many consecutive feasibility cuts and that all used subgradients of Ψ are uniformly bounded. Let δ stop > 0, then the regularized algorithm terminates for some finite k and x k+1 is a (∆ A k + δ stop )-optimal solution to problem [START_REF] Ben-Salem | Gestion Robuste de la production électrique à horizon court-terme[END_REF].

Proof We prove this result by contradiction; we assume that

δ k = (f + Ψ ) (x k ) -f + Ψk (x k+1 ) > δ stop for all k ≥ 1.
If there is an infinite number of iterations where ( 24) holds (called serious iterations, in the bundle terminology), then for all such iterations, we have

(f + Ψ )(x k+1 ) ≤ (f + Ψ )(x k ) -mδ k ≤ (f + Ψ )(x k ) -mδ stop .
Consequently, the fact that the objective function decreases infinitely many times by a fraction of at least mδ stop contradicts the fact that f + Ψ is bounded from below on the compact set X.

Consequently, there exists some k such that x := xk = xk for all k ≥ k. By assumption there exist some C > 0 such that g ≤ C for all g ∈ ∂Ψ (x) and for all x ∈ X. Hence,

δ stop ≤ (f + Ψ ) (x) -f + Ψk (x k+1 ) ≤ (f + Ψ ) (x) -f (x k+1 ) -Ψ (x ℓ ) -g ℓ , x k+1 -x ℓ ≤ (f + Ψ ) (x) -f (x k+1 ) -Ψ (x ℓ ) + C x k+1 -x ℓ (26)
by [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] and the Cauchy-Schwarz inequality. The set X is compact, so by moving to a subsequence if needed, the sequence x k k can be assumed to admit a limit x. Now [START_REF] Geoffrion | Generalized benders decomposition[END_REF] allows us to move to the lim inf and observe:

δ stop ≤ lim inf ℓ≤k,(ℓ,k)→∞ (f + Ψ ) (x) -f (x k+1 ) -Ψ (x ℓ ) + C x k+1 -x ℓ = (f + Ψ ) (x) -lim sup ℓ≤k,(ℓ,k)→∞ f (x k+1 ) -Ψ (x ℓ ) + C x k+1 -x ℓ = (f + Ψ ) (x) -(f + Ψ ) (x), (27)
by continuity of f and Ψ on R n . For each k ≥ k, x k+1 is a null-iterate and therefore:

(f + Ψ ) (x) -(f + Ψ ) (x k+1 ) ≤ m (f + Ψ ) (x) -f + Ψk (x k+1 )
Moving to the lim inf on both sides of the inequality, we establish:

(f + Ψ ) (x) -(f + Ψ ) (x) ≤ m (f + Ψ ) (x) -lim sup k f + Ψk (x k+1 ) , which yields (1 -m) ((f + Ψ ) (x) -(f + Ψ ) (x)) ≤ 0 hence contradicting (27) since m ∈ (0, 1). ⊓ ⊔
The contents of the discussion following Proposition 1 (about the assumptions, inexact resolution of the master in practice, and direct resolution of the master for small problems) still hold for this algorithm.

Algorithmic specificities

Two ingredients of the algorithms play a special role in the practical efficiency of the approach. First, the primal recovery heuristic (Step 2) should have a good balance between computational cost and quality of feasible solutions computed. Second, the bundle method for maximizing the dual of the approximate problem (Step 3) should not be a computational burden. We discuss these two points here.

Primal recovery

Over the course of the dual bundle method, many elements of X are generated. A very simple primal recovery procedure is thus to simply retain the iterate that does best regarding the primal objective function. A more advanced primal recovery would use the convexity of certain subproblems (constrained by X i ), as follows.

Step 1 of the algorithm also provides simplicial multipliers (see e.g., [START_REF] Dubost | A primal-proximal heuristic applied to french unitcommitment problem[END_REF], ) which can be used to form a convex combination of iterates. If X i is convex, this combination will be feasible. An approach for recovery, called convexified best iterate, consists then in considering the best iterate after the convexification. Alternatively one can also solve the Dantzig-Wolfe-like optimization problem (e.g., [START_REF] Briant | Comparison of bundle and classical column generation[END_REF])

min z,x p j=1 m i=1 z j i f i (y j i ) + Ψk (x) (28)
s.t. x i = p j=1 z j i y j i , i = 1, ..., m, p j=1 z j i = 1, i = 1, ..., m Ψ f k (x) ≤ 0, z j i ∈ {0, 1} , i = 1, ..., m, j = 1, ..., p,
where y j , j = 1, ..., p are the solutions obtained while maximizing the Lagrangian dual. Note that, in the above problem [START_REF] Ghaoui | Robust solutions to uncertain semidefinite programs[END_REF], f i (y j i ) are fixed and x is an auxiliary variable, so that the problem ( 28) is a mixed-integer linear program for which efficient greedy heuristics exist. For example, this procedure was exploited in [START_REF] Takriti | Using integer programming to refine lagrangian-based unit commitment solutions[END_REF] in the context of unit-commitment.

Warm-starting the master problem resolution

Maximizing the Lagrangian dual function [START_REF] Oliveira | A bundle method for exploiting additive structure in difficult optimization problems[END_REF] or [START_REF] Dubost | A primal-proximal heuristic applied to french unitcommitment problem[END_REF] requires in turn building a cutting plane model of them. These cutting plane models are the essential building block of a bundle method (see e.g., the textbook [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms II[END_REF] or appendix A). Since most of the computational cost resides in solving the m subproblems, it is useful to carry as much cutting-plane information as possible when moving from iteration k to iteration k + 1. This is the idea behind the following warm-starting procedure for the master problem. Observe that we have two identifications. First, for any k ≥ 1 and any

µ ∈ R k Θ k+1 ((µ, 0)) = Θ k (µ).
Second, for any iterate µ ∈ R k with associated x(µ) ∈ X encountered while maximizing Θ R k , we have

Θ R k+1 ((µ, 0)) = f (x(µ)) + 1 2t k+1 x(µ) -xk+1 2 2 + k ℓ=1 µ ℓ g ℓ , x(µ) + k ℓ=1 µ ℓ Ψ (x ℓ ) -g ℓ , x ℓ = Θ R k (µ)
if t k+1 = t k and xk+1 = xk . When t k+1 and xk changed, we can still interpret x(µ) as an inexact solution to the optimization problem [START_REF] Dubost | A primal-proximal heuristic applied to french unitcommitment problem[END_REF] at dual signal (µ, 0) for iteration k + 1. This gives rise to the usual rationale of an "upper"-oracle (see [START_REF] Oliveira | A bundle method for exploiting additive structure in difficult optimization problems[END_REF]Lemma 2.1] and [START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF]) since we are working in the concave dual.

In practice, we just need to save past primal iterates with their true function values to set up an initial cutting plane model for Θ k+1 or Θ R k+1 following Lemmas 1 and 3. This strategy proves highly effective, as the computational cost of maximizing the Lagrangian dual for a given k is highly reduced. This is illustrated in Figure 4.

Possible extensions to inexact oracles for Ψ

We finally note that our scheme can be combined with inexact computations of the slave problem, which is particularly useful if the set D is hard to maximize over. Indeed, we have assumed so far that we can solve sup d∈D ψ(x, d) for a given x; in other words, that the oracle for Ψ is exact. As mentioned in the introduction, performing this maximization is not mandatory, since any d can provide us with valid linearizations for Ψ by [START_REF] Ben-Tal | Robust solutions of uncertain linear programs[END_REF]. This inexactness (i.e., computing Ψ inexactly) can be handled by bundle methods. Following [START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF]Lemma 5.1], we need to append the test δ k < 1 2t k x k+1 -xk 2 in between steps 2 and 3 of the algorithm Whenever this condition holds true, the presence of noise is suspected and we suggest to apply the usual remedy [START_REF] Kiwiel | A proximal bundle method with approximate subgradient linearizations[END_REF] which consists in increasing t k strictly and then return to Step 1 to produce a new iterate. The analysis of the effect of using inexact oracles of [START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF]. reveals that only when nearing convergence, it is of importance to request some precision when computing Ψ to get a good final accuracy. Requesting more precision here means maximizing the mappings ψ(x, d) sufficiently well over D. We refer the reader to [START_REF] Oliveira | Level bundle methods for oracles with on demand accuracy[END_REF][START_REF] Van Ackooij | Probabilistic optimization via approximate p-efficient points and bundle methods[END_REF] for further details and alternatives. We note that in our numerical experiments computing Ψ is fast so that we have not explored this point any further.

Numerical experiments

The goal of the numerical experiments is twofold: first, illustrating that the approach of the paper can be deployed on a real-life application; second, comparing the cutting plane method versus the stabilized cutting plane method, showing the interests of the latter.

Experimental setting

Our numerical experiments are conducted on robust large scale unit-commitment problems as presented in section 2.1, similar to those that EDF, faces everyday. We briefly present here our unit-commitment model and data.

Unit-commitment model and data

We use standard choices in the unit-commitment literature for unit modeling, i.e., a convex model for hydro valleys (e.g., [START_REF] Van Ackooij | Joint chance constrained programming for hydro reservoir management[END_REF]) and a standard non-convex model for thermal units (e.g., [START_REF] Frangioni | Sequential Lagrangian-MILP Approaches for Unit Commitment Problems[END_REF]). The overall model comes from a part of the set of French electricity generation assets. The total number of variables describing the set X is around 50000 continuous variables, 27000 binary variables and we need around 815000 constraints. The data is taken from the operational data sets of 2013 of EDF: we have T = 96 half hourly time steps, and 6 data sets for 6 days of 2013, given by Table 1. Moreover we generate a variant of each data set by perturbing the load data with a multiplicative factor along the lines of [60, Section 5.1] with three dispersion load scenarios (low, medium, and large). We hence obtain a total of 18 instances defined by its date and load factor.

Robust unit-commitment model In this above-described large scale unit-commitment setting, we consider the robust 2-stage model of the form (3) as described in section 2.1. The coupling penalization function Ψ is defined from the piecewise affine functions ψ t of (10) with EDF industrial choices for parameters a i,t and b i,t (see [START_REF] Dubost | A primal-proximal heuristic applied to french unitcommitment problem[END_REF]Fig 2.2]).

For the model of uncertainty, our approach could handle any set D for which we can easily solve [START_REF] Borghetti | Lagrangian heuristics based on disaggregated bundle methods for hydrothermal unit commitment[END_REF]. This includes in particular finite sets with small cardinality and Cartesian products of intervals. However these two models have a limited expressivity: on one hand the finite set is obviously limited; and on the other hand, the product of intervals is too large and does not express the dependencies between two successive time-steps that we have in our application. Then we choose the so-called "state-space represented" model as suggested in [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF]. This model of uncertainty is finite but of very high cardinality, expressing temporal dependencies while preserving a very fast computability by dynamic programming. This model is illustrated and shortly described in appendix C; for a complete analysis, we refer to [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF].

Lagrangian approaches vs monolithic resolution approaches

In our context, the deterministic problem (i.e., without the robust penalization function ψ but with load-demand balance equations) is so large with the data instances described above, that it cannot be solved (to any decent precision) by a monolithic solver (e.g., CPLEX). On top of this, the state-space model for D contains an exponential number of scenarios, in our case 1.35 × 10 81 , so that an explicit reformulation involving all these scenarios is simply not possible. The Lagrangian dual approach is therefore the only viable alternative to solving the master problem. In general, as described in [START_REF] Tahanan | Large-scale unit commitment under uncertainty: a literature survey[END_REF], Lagrangian dualization is one of the essential methods for solving large scale unit-commitment problems with extensive modeling detail. The Lagrangian approach produces lower-bounds and (infeasible) primal solutions of very good quality exploited by efficient heuristics. In related contexts using appropriate methods (e.g., [START_REF] Borghetti | Lagrangian heuristics based on disaggregated bundle methods for hydrothermal unit commitment[END_REF][START_REF] Feltenmark | Dual applications of proximal bundle methods, including lagrangian relaxation of nonconvex problems[END_REF][START_REF] Takriti | Using integer programming to refine lagrangian-based unit commitment solutions[END_REF]) solutions of very good quality, near 0.5 % optimality, can be obtained (e.g., [START_REF] Frangioni | Sequential Lagrangian-MILP Approaches for Unit Commitment Problems[END_REF]).

Fig. 2 Comparison of the three generic heuristics on a typical run of the bundle method within the solving process Fig. 3 Illustration of the oscillation effect of the cutting plane algorithm versus the stabilized one on instance number 6

Algorithmic choices We use rudimentary choices of parameters after some preliminary tests. In particular we use a fixed prox-parameter for all data sets (which is theoretically justified, see the discussion in the appendix on bundle methods). We empirically chose the value t k = 1000 since it tends to "balance off" the terms in the objective function of problem [START_REF] Van Dinter | The unit commitment model with concave emissions costs: a hybrid benders decomposition with nonconvex master problems[END_REF]. We also use the Dantzig-Wolfe heuristic (see section 3.4) as recovery step in the two algorithms, since it tends to produce much better feasible solutions. We illustrate this on a typical example in Figure 2.

Numerical results

We compare the stabilized cutting plane method of section 3.3 against the cutting plane method of section 3.2. All computations have been carried out using CPLEX 12.6 on an Intel Xeon X5670 westmere computer cluster with 8 Gb of reserved memory without using parallelization.

We consider large instances: the cost of the master problems are negligible in the total computing amount and the main computational cost in the solution process consists of solving the m subproblems (2). So we count the cost of the algorithms in terms of the total number of sub-problem resolutions, and the comparison is based on this cost. Figure 3 shows a typical behaviour of the two algorithms on one instance. It immediately becomes apparent that the cutting plane method suffers significantly from oscillation and that the regularized variant shows a better behaviour as this aspect is concerned. Tables 2 and3 report the complete comparison on all the instances with two different stopping thresholds. The value -corresponds to a computation that did not reach convergence within the time limit. The tables clearly show that stabilization reduces the number of (outer-)iterations, and most importantly the number of oracle calls. Even instances that were not solved by the basic cutting plane method, can be solved by the stabilized method. The impressive gain of the stabilized algorithm could be explained by two factors. First, the stabilized algorithm produces better iterates, which is the raison-d'être of this method, but not always easy to obtain in a context of mixed-integer programming ( [START_REF] Zaourar | Optimisation convexe non-différentiable et méthodes de décomposition en recherche opérationnelle[END_REF], where typically the method is not outperformed uniformly on all instances). Second, the strategy of the warm-starting presented in section 3.4, though highly effective for both algorithms, is even more efficient for the stabilized algorithm. We illustrate the effect of warm-starting on a typical example on Figures 4 and5. We see on Figure 4 that the computational cost of maximizing the Lagrangian dual for a given k is highly reduced for the typical example. Figure 5 then illustrates that the warm start is more efficient for the stabilized algorithm which, roughly speaking, explores the solutions in a more continuous way. The shown efficiency of warm-starting (under regularization) seem to be observations of some sort of general mechanism. Indeed, the phenomenon was already reported in [START_REF] Van Ackooij | Decomposition approaches for block-structured chance-constrained programs with application to hydro-thermal unit commitment[END_REF] but also in [START_REF] Zaourar | Optimisation convexe non-différentiable et méthodes de décomposition en recherche opérationnelle[END_REF] (less feasibility cuts).

Conclusions

We have proposed and studied a generic decomposition scheme for robust block-structured problems. The main algorithm is a regularized cutting-plane method intertwining two layers of bundle methods that scales up to large instances. Its practical performance relies on the combination of efficient numerical features, including warm-starting and heuristics. Its theoretical analysis is not covered by existing results because of the technical difficulty of the inexact resolution of the master problem. We have provided a convergence analysis using basic arguments. We have illustrated the approach with numerical experiments on large scale instances coming from unit-commitment. Though cutting-plane methods are basic methods in robust optimization, regularized cutting-plane algorithms are not often presented or studied in this context; in contrast, it is sometimes reported to be not useful. Our numerical experiments show clearly that for the robust block-structured optimization problems we consider here, the decomposition scheme is much more efficient with regularization, in terms of computing times and stability. We also underline that the numerical features of the two cutting-plane algorithms are essentially the same, and we advocate that the jump from the standard to the regularized one is not complicated, neither conceptually nor implementation-wise. The main contribution of this paper is thus to promote the combination of cutting-plane and regularization in the context of robust optimization. To support this message we release1 the code implementing the regularized algorithm.

A Brief overview of bundle methods in nonsmooth optimization

This appendix provides a brief description of bundle methods using standard notation and terminology, which are also the ones of this paper. For an in-depth description, we refer to the textbook [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms II[END_REF] and the recent article [START_REF] Oliveira | Convex proximal bundle methods in depth: a unified analysis for inexact oracles[END_REF].

We are interested in solving minx f (x), where f : R n → R is a nonsmooth convex mapping known through an "oracle" returning, for a given x ∈ R n as an input, the value f (x) and a subgradient g ∈ ∂f (x) (or approximations of them). Given a set of iterates x 1 , ..., x k for which the oracle has been called, we can set up the cutting plane model fk for f by defining fk

(x) := max ℓ=1,...,k f (x ℓ ) + g ℓ , x -x ℓ ,
where g ℓ ∈ ∂f (x ℓ ). The cutting-plane algorithm, which can be traced back to [START_REF] Kelley | The cutting-plane method for solving convex programs[END_REF], consists of defining the next iterate by solving minx fk (x), which can be formulated as a linear problem.

The cutting-plane method suffers from various drawbacks (including stability and memory issues) and an improved scheme is available: regularized cutting-plane methods or so-called bundle methods, dating back to the seminal work of [START_REF]An extension of davidon methods to nondifferentiable problems[END_REF]. These methods use a quadratic stabilization to ensure that the next iterate is not far from the best current iterate, called stability center x. The next iterate is the solution of

min x fk (x) + 1 2t k x -xk 2 , (29)
which is a convex quadratic program (easy to solve by specialized methods [START_REF] Kiwiel | A method for solving certain quadratic programming problems arising in non-smooth optimization[END_REF][START_REF] Kiwiel | A cholesky dual method for proximal piecewise linear programming[END_REF]). The predicted decrease δ k := f (x k )fk (x k+1 ) is then used to update the stability center: if enough decrease is observed, i.e. the inequality

f (x k+1 ) ≤ f (x k ) -mδ k . (30)
holds, the new iterate x k+1 becomes the new stability center xk+1 ; Otherwise the information from the oracle is just used to improve the cutting plane model. From a variational analysis perspective, bundle methods are thus damped proximal algorithms where the inner proximal iterations are stopped whenever sufficient decrease is observed on f .

Note that the proximal parameter t k is basically free: for exact oracles, convergence is guaranteed if t k stays within bounds (see [START_REF] Correa | Convergence of some algorithms for convex minimization[END_REF]). In practice, t k can be arbitrary if the stability center was changed; t k should not increase if the stability center was not changed; t k should strictly increase if the oracle is subject to imprecise calculations during noise detection steps (see e.g., [61, eq. (2.16)]). A complete description of asymptotical rules can be found in [START_REF] Van Ackooij | Incremental bundle methods using upper models[END_REF]).

On top of stabilizing iterates, bundle methods also come up with the possibility of aggregating cutting-plane information into an aggregate cut. In theory, we can delete all but two cutting planes (the aggregate cut and the cut belonging to the current stability center).

B Benders decomposition: terminology and references

The so-called Benders decomposition hinges on the observation that fixing some variables makes certain problems much easier. This observation, originally made for problems with an underlying linear structure [START_REF] Benders | Partitioning procedures for solving mixed-variables programming problems[END_REF] potentially with stochastic aspects [START_REF] Van Slyke | L-shaped linear programs with applications to optimal control and stochastic programming[END_REF], was generalized to problems with some underlying convexity in the seminal work [START_REF] Geoffrion | Generalized benders decomposition[END_REF]. We recall in this appendix the terminology and some important references about this classical method in operation research in a broad sense.

The slave problem consists in the resulting problem with the fixed variables. The master problem is related to finding the optimal allocation of the previously fixed variables. Since the domain of the value function need not be the whole space, the master problem is augmented with the so called feasibility cuts (i.e., an outer approximation of the convex domain of the value function). The master problem is also enriched with optimality cuts, as a matter of fact, a cutting-plane model of the value function. Benders algorithm is thus as a variant of the cutting-plane method, presented previously, and therefore subject to the well-known oscillation effect and slow convergence. Regularized Benders algorithms have been recently proposed [START_REF] Zaourar | Quadratic stabilization of benders decomposition[END_REF] and [START_REF] Van Ackooij | Inexact stabilized Benders' decomposition approaches: with application to chance-constrained problems with finite support[END_REF].

Many articles, including [START_REF] Van Dinter | The unit commitment model with concave emissions costs: a hybrid benders decomposition with nonconvex master problems[END_REF][START_REF] Sahiridis | Accelerating benders method using covering cut bundle generation[END_REF][START_REF] Sherali | On generating maximal nondominated benders cuts[END_REF][START_REF] Wentges | Accelerating benders' decomposition for the capacitated facility location problem[END_REF][START_REF] Yang | A tighter cut generation strategy for acceleration of benders decomposition[END_REF] concretely deal with strategies for generating good (optimality) cuts, occasionally with a problem-dependent flavour. The authors of [START_REF] Codato | Combinatorial benders' cuts for mixed-integer linear programming[END_REF] are concerned with generating strong feasibility cuts and provide an important contribution in this view. Further improvements to the general scheme are concerned with a relaxation of the convexity assumption of the slave problem by allowing for integer variables [START_REF] Caroe | L-shaped decomposition of two-stage stochastic programs with integer recourse[END_REF] or the use of generalized "logic" duality (called inference duality) in [START_REF] Hooker | Logic-based benders decomposition[END_REF]. In [START_REF] Zakeri | Inexact cuts in benders decomposition[END_REF], it is moreover suggested to inexactly solve the slave problem to compute some inexact but "cheap" cuts. In general, an appropriate choice of cuts (approximating the value function) is crucial, as illustrated in [START_REF] Kolokolov | Analysis of decomposition algorithms with benders cuts for p-median problem[END_REF].

Let us finally illustrate optimality and feasibility cuts more specifically on the following parametrized linear program: ψ(x) := min y q T y : W y ≤ b -T x .

The inequality system W y ≤ b -T x need not admit a solution for all x and that consequently ψ(x) = ∞ may happen for some x. Minimizing the convex mapping ψ then amounts to handling also the hidden constraints x ∈ dom ψ. By strong duality, we have ψ(x) = max λ≥0 λ T (b -T x) : W T λ ≥ q , so that any optimal dual multiplier λ * satisfies -T T λ * ∈ ∂ψ(x).

When W y ≤ b -T x does not admit a feasible solution, by Farkas Lemma an unbounded ray can be identified and used to define an outer linearization for dom ψ. Note finally that, in the general non-linear situation, first order information about ψ can also be identified with tools from nonsmooth analysis (e.g., [START_REF] Rockafellar | of Grundlehren der mathematischen Wissenschaften[END_REF]Theorem 10.13]), but feasibility cuts will need to be derived by solving an auxiliary problem typically involving the minimization of slack variables.

C A rich but computationally cheap model of uncertainty

For our numerical experiments, we choose to use the uncertainty set D introduced by [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF]. This so-called "state-space represented" model has a rich expressivity while being computationally easy to manipulate. This model has an exponential number of scenarios expressing temporal dependencies, while preserving a very fast computability by basic dynamic programming. We refer to [START_REF] Minoux | Two-stage robust optimization, state-space representable uncertainty and applications[END_REF] more details and a complete analysis. This model is illustrated in Figure 6. Yet, for the convenience of the reader, we describe here the main ideas briefly: At each time step t ∈ τ , we set up a set of nodes Et containing both a value for dt and a weight wt. This set of nodes is assumed to be sorted according to increasing values for dt. Using this set of nodes we set up a graph (G, V ), with G = t∈τ Et, such that V connects all nodes in Et to those in E t+1 that are not further apart than H, i.e., node i ∈ Et is connected to node j ∈ E t+1 if and only if their local labels |L(i) -L(j)| ≤ H. At least one of the nodes in Et is assumed to have a zero weight and H is assumed to be such that each

1 Introduction 1 . 1

 111 Robust block-structured optimization We consider block-structured optimization problems of the following form in R n (with n = n 1 + • • • + n m ) min x=(x1,...,xm) m i=1

1. 2

 2 Models of uncertainty, reformulation, and cutting-planes Formally, the robust optimization problem (3) has the same structure as the initial problem (1) by just replacing ψ with the convex function Ψ : R n → R defined as Ψ (x) := sup d∈D ψ(x, d).

Fig. 1

 1 Fig. 1 Example of a penalization function ψt

  ). -(Step 3 : Oracle call) : Call the oracle for Ψ at x k+1 to enrichen the model for Ψ by adding a (feasibility or optimality) cut. -(Step 4 : Stopping test) : In the case where x k+1 lies in dom Ψ , stop the algorithm if

Fig. 4

 4 Fig.4Comparison of the number of inner iterations (within Step 1) for each outer iteration, with and without warm-starting

Fig. 5

 5 Fig. 5 Less subproblem resolutions for quadratic algorithm using warm-start

Fig. 6

 6 Fig. 6 Illustration of the uncertainty model used in the numerical experiments

Table 1

 1 Instance label with corresponding date and numbers of units; the data sets contain 96 half hourly time steps.

	Instance	date	# Thermal units # Hydro valleys
	1	15/01/2013	106	41
	2	20/03/2013	106	32
	3	13/05/2013	104	35
	4	22/08/2013	104	34
	5	25/10/2013	104	37
	6	10/12/2013	104	40

Table 2

 2 Comparison of the cutting plane method (column cp) and stabilized cutting plane method (column stab cp) for a stopping criteria δstop = 5%.

			cp			stab. cp	improvement
	Instance	load factor	# It. Oracle Calls	# It. Oracle Calls	# It.	Oracle Calls
		low	64	711	43	449	149%	158%
	15/01/2013	medium	45	760	44	495	102%	154%
		high	41	536	26	317	158%	169%
		low	2045	2281	251	416	815%	548%
	20/03/2013	medium	848	1082	157	317	540%	341%
		high	393	613	103	245	382%	250%
		low	3074	3862	439	964	700%	401%
	13/05/2013	medium	1608	2404	266	720	605%	334%
		high	764	2315	319	974	239%	238%
		low	-	-	684	1098	∞%	∞%
	22/08/2013	medium	-	-	662	1152	∞%	∞%
		high	2648	3452	356	790	744%	347%
		low	-	-	297	815	∞%	∞%
	25/10/2013	medium	2923	3894	369	1044	792%	373%
		high	724	2163	209	1147	346%	189%
		low	236	916	96	398	246%	230%
	10/12/2013	medium	136	661	55	321	247%	206%
		high	78	578	34	229	229%	252%
	average	1042	1749	245	661	420 %	285 %

Table 3

 3 Comparison of the cutting plane method (column cp) and stabilized cutting plane method (column stab cp) for a stopping criteria δstop = 1%.

https://github.com/NLebbe/rbsopb

Acknowledgements Part of this work was done during the second author's internship at EDF R&D. The first and the third author gratefully acknowledge the support of PGMO project "advanced nonsmooth optimization methods for stochastic programming".