Machine Learning-Guided Dual Heuristics and New Lower Bounds for the Refueling and Maintenance Planning Problem of Nuclear Power Plants

: This paper studies the hybridization of Mixed Integer Programming (MIP) with dual heuristics and machine learning techniques, to provide dual bounds for a large scale optimization problem from an industrial application. The case study is the EURO/ROADEF Challenge 2010, to optimize the refueling and maintenance planning of nuclear power plants. Several MIP relaxations are presented to provide dual bounds computing smaller MIPs than the original problem. It is proven how to get dual bounds with scenario decomposition in the different 2-stage programming MILP formulations, with a selection of scenario guided by machine learning techniques. Several sets of dual bounds are computable, improving signiﬁcantly the former best dual bounds of the literature and justifying the quality of the best primal solution known.


Introduction
Hybridizing heuristics with mathematical programming and Machine Learning (ML) techniques is a prominent research field to design optimization algorithms, these techniques having complementary advantages and drawbacks [1]. Matheuristics, hybridizing mathematical programming and heuristics, use the modeling facilities of mathematical approaches for industrial optimization problems, possibly highly-constrained ones which are usually difficult for meta-heuristics, with a good scalability to face large size instances, which is a characteristic of meta-heuristics and the bottleneck for exponential methods like Branch&Bound (B&B) tree search [2].
Usually, hybrid algorithms compute primal solutions of optimization problems as in [3], rare works use such hybridization to improve dual bounds. Using dual heuristics for Mixed Integer Linear Programming (MILP) is mentioned in [4]: "In linear programming based branch-and-bound algorithms, many heuristics have been developed to improve primal solutions, but on the dual side we rely solely on cutting planes to improve dual bounds". Until the 2000s, many efforts on improving the quality of dual bounds relied on improving cut generation to speed up B&B approaches [5] or using Lagrangian bounds like the Danzig-Wolfe reformulation [6]. In the 2000s, primal heuristics inside B&B search improved exact MILP solving, having earlier good primal solutions in the B&B search

Problem Statement
Planning the maintenance of nuclear power plants is not a pure scheduling problem: the impact of maintenance decisions is calculated with the expected production cost to fulfill the power demands. Maintenance operations for nuclear power plants imply an outage, a null production phase during several weeks. Planning the outages is crucial in the overall production costs: the part of nuclear production in the French power production is around 60-70%, and the marginal costs of nuclear production are lower than the ones of other thermal power plants [18,19]. Furthermore, maintenance planning must imply the feasibility of low-level technical constraints for productions and fuel levels. The optimization problem was formulated using 2-stage stochastic programming in [14]. Power demands, production capacities, and costs are uncertain and modeled using discrete stochastic scenarios. Refueling and maintenance decisions for nuclear power plants are taken, as a tactical level. The operational level optimizes production plans implied by the previous decisions for each stochastic scenario, to minimize the average production cost. Notations are presented in Table 1 for the sets and indexes and in Table 2. Table 1. Definitions and notations-the sets and indexes.
c ∈ C CT with CT ∈ [ [14,21]] Scheduling constraints denoted CT14 to CT21 in [14].  Proportional cost to the refueled quantity at outage k. C p j,s,t Production costs for T1 unit j at time period t and scenario s. D s,t Power demands at time step t for the scenario s.
Outage duration for maintenance and refueling at cycle k. F t = F Conversion factor between power and fuel in time step t, constant in the data. I c,w Maximal offline T2 power at week w for CT21 constraint c. L i,k,c Latency in the resource requirement in CT19 constraint c for outage i, k. N c,w Maximal number of outages under maintenance at week w for CT20 constraint c. π s Probability of scenario s. Pmin j,s,t Minimal power for T1 unit j at time period t and scenario s. Pmax j,s,t Maximal power for T1 unit j at time period t and scenario s.
Maximal generated power at time step t .
Proportion of fuel that can be kept during reload in cycle k at plant i R 19 c Maximal quantity of resource in CT19 constraint c.

Rmin i,k
Minimal refueling at outage k.

Rmax i,k
Maximal refueling at outage k. S c Spacing/overlapping value defined by CT14-CT18 constraints. Smax i,k Maximal fuel level at production cycle k. t w Index of the first time step of week w. Ta i,k Last possible beginning week for outage k of T2 plant i. To i,k First, possible outage week for cycle k.
Length of the resource usage in CT19 constraint c for outage i, k. w t Week of the production time step t. Xi i Initial fuel stock of T2 plant i.

Production Assets
To generate electricity, two kinds of power plants are distinguished. On one hand, Type-1 (shortly T1, denoted j ∈ J ) power plants can be supplied in fuel continuously without inducing offline periods. On the other hand, Type-2 (shortly T2, denoted i ∈ I) power plants have to be shut down for refueling and maintenance regularly. T2 units correspond to nuclear power plants. The production planning for a T2 unit is organized in a succession of cycles, an offline period (called outage) followed by an online period (called production campaign). Cycles are indexed with k ∈ K = [[0, K]], each T2 unit i having possibly K maintenance checks scheduled in the time horizon. The number of maintenance checks planned in the time horizon may be inferior to K. Cycle k = 0 denotes initial conditions, the current cycle at w = 0. For units on maintenance at w = 0, cycle k = 0 considers the remaining duration of the outage. For units on production at w = 0, notations and constraints are extended considering a fictive cycle k = 0, with a null duration.

Time Step Discretization
The time horizon is discretized with two kinds of homogeneous time steps. On one hand, outage decisions are discretized weekly and indexed with w ∈ W = [[1; W]]. On the other hand, production time steps for T1 and T2 units are discretized with t ∈ T , hourly time periods from 8 h to 24 h. This fine discretization is used to consider fluctuating demands in hourly periods [20,21].

Decision Variables and Objective Function
Maintenance, refueling, and production decisions are optimized conjointly. The variables for outage dates and refueling levels are related to T2 units, whereas the production variables concern T1 and T2 units. There is another heterogeneity of maintenance and production decisions: production decisions are considered for a discrete set of scenarios s ∈ S to model uncertainty on power demands and production capacities, whereas outage and refueling decisions are common for all the scenarios, and must guarantee the feasibility of a production planning for each scenario. The objective function minimizes the average value of the production costs considering all the scenarios. T1 Production costs are proportional to the generated power. T2 production costs are proportional to the refueling quantities over the time horizon, minus the expected stock values at the end of the period to avoid end-of-side effects.

Description of the Constraints
There are 21 sets of constraints in the Challenge, numbered from CT1 to CT21 in the specification [14]. CT6 and CT12 are relaxed in our study; this is possible as we compute dual bounds relaxing these two sets of constraints. Table 3 describes precisely the constraints.
Constraints CT1 to CT12 concern the production and fuel constraints in the operational level whereas CT13 to CT21 are scheduling constraints in the strategic level. Constraints CT1 couple the production of the T1 and T2 plants with power demands. Constraints CT2 and CT3 define the bounds for T1 productions. Constraints CT4 to CT5 describe the production possibilities for T2 units depending on the fuel stock with outages during the maintenance periods. Constraints CT7 to CT11 involve fuel constraints for T2 units: bounds on fuel stocks, on fuel refueling, and relations between remaining fuel stocks and production/refueling operations.
Constraints CT13 impose Time Windows (TW) for the beginning weeks of outages, and it can impose some maintenance dates. With CT13 constraints, the maintenance checks follow the order of set k ∈ K without skipping maintenance checks: if an outage k + 1 is processed in the time horizon, it must follow the production cycle k. Constraints CT14 to CT21 may be unified in a common format with generic resources. These constraints express minimal spacing/ maximal overlapping constraints among outages, and minimum spacing constraint between coupling/decoupling dates, resource constraints for maintenance using specific tools or manpower skills, and limitations of simultaneous outages with a maximal number of simultaneous outages or maximal T2 power offline. CT14: For all constraint c ∈ C 14 , a subset of outages, Z c have to be spaced by at least S c weeks: for S c > 0, it is a minimal spacing from the beginning of a previous outage to the beginning of a next outage, for S c 0, it is a maximal number of weeks where two outages can overlap.
CT15: CT15 are similar to CT14 with sets C 15 , Z c and parameters S c , applying only on a time interval W c .
CT16: for each constraint c ∈ C 16 , a subset of outages Z c have to be spaced by at least S c > 0 weeks from the beginning of the previous outage to the beginning of the next outage.
CT17: for each constraint c ∈ C 17 , a subset of outages Z c have to be spaced by at least S c > 0 weeks from the end of the previous outage to the end of the next outage.
CT18: for each constraint c ∈ C 18 , a subset of outages Z c have to be spaced by at least S c > 0 weeks from the end of the previous outage to the beginning of the next outage.
CT19: for each constraint c ∈ C 19 , a subset of outages (i, k) ∈ Z c shares a common resource (which can be a single specific tool or maintenance team) with forbid simultaneous usage of this resource. L i,k,c and Tu i,k,c indicate respectively the start and the length of the resource usage and R 19 c the maximal quantity of the resource.
CT20: for each constraint c ∈ C 20 , at most N c,w outages among a subset Z c overlap in weeks w ∈ W c . CT21: for each constraint c ∈ C 21 , a time period W c is associated where the offline power of T2 plants I c ⊂ I due to maintenance operations must be lower than I c,w .

Related Work
This section describes the solving approaches for the problem, focusing on the dual bounds possibilities. We refer to [22] for a general survey on maintenance scheduling in the electricity industry, and to [23,24] for specialized reviews on the approaches for the EURO/ROADEF Challenge 2010.

Solving Methods
Three types of approaches were designed for the EURO/ROADEF Challenge 2010 competition. MILP-based exact approaches used MILP models to tackle a simplified and reduced MIP problem before a reparation procedure to build solutions for the original problem, as in [25][26][27]. Frontal local search approaches used neither simplification nor decomposition in local search approaches. These were the most efficient primal heuristics for the Challenge with outstanding results compared to the other approaches. An unpublished simulated annealing approach had most of the best primal solutions for the challenge [28]. An aggressive local search, LocalSolver's forebear [29], provided most of the reactualized best primal solutions [15]. Heuristic decomposition approaches separate the maintenance and refueling planning decisions from the production optimization, and computes productions independently for each scenario s, as in [30][31][32][33]. Such approaches are less efficient. This is explained by [24].

Reductions by Pre-Processing
Facing the large size of the instances, most of the approaches competing in the EURO/ROADEF Challenge 2010 reduced the problem with pre-processing strategies. Fixing implied or necessarily optimal decisions is a valid pre-processing to compute dual bounds. Heuristic pre-processing strategies were also commonly used. Most of the previous approaches reduced instances using aggregation techniques. Firstly, the power productions are discretized with many time steps t ∈ T , and may be aggregated to their weekly average value, so that production and maintenance decisions have the same granularity. Under some hypotheses, it is proven such aggregation induces lower bounds of the original problem [17]. Secondly, many approaches used computations reduced on single scenarios, aggregating the scenarios into one average deterministic scenario as in [27], or using decomposition as in [25,26,[30][31][32][33]. It is proven that lower bounds of the EURO/ROADEF Challenge 2010 can be guaranteed with S single scenario computations [17].
TW constraints CT13 can be tightened exactly, and we refer to [34] for the most detailed presentation of such pre-processing. Minimal lengths for production campaigns are implied by the initial fuel levels defined in CT8 constraints, maximal T2 productions with CT5 constraints, maximal stock before refueling with CT11 constraints. Such minimal lengths for production cycles imply that TW constraints can be tightened by induction. Tightening TW constraints, some outages can be removed when their earliest completion time exceed the time horizon. In the problem specification, there is no maximal length for a production cycle, and it is possible to have a null production phase when the fuel is consumed. However, such null production phases are not economically profitable as T2 production costs are lower than T1 production costs. Adding maximal lengths for production cycles is a heuristic pre-processing that was proven efficient for primal heuristics [26].
Exact pre-processing applies also for T1 production variables. Once the T2 production is given, optimizing T1 production is equivalent to solve independently knapsack problems for each time step. A greedy algorithm computes optimal T1 production by sorting the marginal costs. T1 units can be aggregated within one unit with piecewise linear costs as in the greedy optimization of T1 production. The most costly T1 units (or equivalently the last points of the piecewise linear cost function) may be removed when the production levels that are not reached considering the demands in power.

Dual Bounds for the Euro/Roadef Challenge 2010
The first dual bounds for the EURO/ROADEF Challenge 2010 were published by [16] and designed independently from their primal heuristic. Two methods were investigated, computing optimal solutions of relaxations. A first method relaxes power profile constraint (CT6), as well as all fuel level tracking (CT7 to CT12) and outage scheduling constraints (CT13 to CT21). The remaining computation to optimality can be processed greedily assigning the production levels to the cheapest power plants for all scenarios. A second method uses a flow network relaxation which considers outage restrictions (CT13) as well as fuel consumption in an approximate fashion (CT7 to CT12) to deduce a tighter lower bound of the objective function than the simple greedy approach. The computational effort to compute the bounds differs dramatically: while the greedy bounds are computed in less than 10 s, solving the flow-network for all scenarios takes up to an hour for the biggest instances.
A preliminary work for that paper improved significantly the previous lower bounds [17]. Using an MILP formulation relaxing only constraints CT6 and CT12, single scenario computations and possibly aggregation of production time steps, it is tractable to compute dual bounds for all the instances given for the Challenge allowing one hour per MILP computation. For the most difficult instances, a dual heuristic using an additional surrogate MILP relaxation improves the results [17].
None of the exact methods used to design primal heuristics for the Challenge published dual bounds, and some reasons will be analyzed in the next section. We note that Semi-Definite Programming (SDP) relaxations were investigated to compute dual bounds for the problem [35]. However, the large size of the instances is still a bottleneck for such SDP approaches [35,36].

Milp Formulations
Several competing approaches for the Challenge were based on MILP formulations. Many constraints have indeed a linear formulation in the specification [14]. There were three types of MILP solving: straightforward B&B based on a compact formulation, a Bender's decomposition, and an extended formulation for a Column Generation (CG) approach.
Several works noticed that scheduling constraints from CT14 to CT21 are modeled efficiently with MILP formulations using time indexed formulation of constraints [25,26]. Straightforward MILP solving is efficiently modeling only CT14 to CT21 constraints [25]. We note that a preliminary work designed a compact MILP model with an exact description of all the constraints (including the CT6 and CT12 ones) but with too many binary variables to hope for a practical efficiency [37].
Only one exact method did not aggregate the stochastic scenarios in a Bender's decomposition [26]. Relaxing constraints CT6 and CT12, an MILP formulation is designed with binary variables only for the outage decisions. Bender's master problem concerns the maintenance dates and the refueling quantities, whereas independent Bender's sub-problems are defined for each stochastic scenario with continuous variables for productions and fuel levels. Even with the weekly aggregation of production time steps, limitations in memory usage and computation time do not allow for deploying entirely the Bender's decomposition algorithm. The heuristic of [26] computes the LP relaxation exactly using Bender's decomposition, a cut and branch approach repairs integrity, branching on binary variables without adding new Bender's cuts, and solutions are repaired considering CT6 and CT12 constraints. The resulting matheuristic approach was efficient for the small instances, and difficulties and inefficiencies occur for the real size instances. Knowing that the time step aggregation keeps dual bounds thanks to [17], dual bounds can be computed after the LP relaxation and the cut and branch phase. However, no such dual bounds are reported by [26].
An exact formulation of CT6 constraints was considered in a CG approach dualizing coupling constraints among units and aggregating the stochastic scenarios to the average one [27]. Production time steps were also aggregated weekly. These two reductions are not prohibitive to compute dual bounds, contrary to a third one: the production domains are discretized to solve CG subproblems by dynamic programming. This is a heuristic reduction of the feasible domain, and it does not provide dual bounds for the original problem. A CG approach was also investigated to complexify the model after the Challenge in [19], with a robust and multi-objective extension searching considering delays in the maintenance operations, as in [38]. In the recent work [19], CG computations are tractable restricting computations to horizons of 2-3 years, instead of the 5-year time horizon of the Challenge.

Mathematical Programming Relaxations
This section presents MILP formulations and reduction techniques that allow to compute dual bounds for the EURO/ROADEF Challenge 2010.

Relaxing Only Constraints CT6 and CT12
Relaxing CT6 and CT12 constraints, an MILP formulation can be designed with binary variables only for the outage decisions, similar to [26]. The binaries x i,k,w in [26] are equal to 1 if and only if the beginning week for cycle (i, k) is exactly w. Similar to [20], we define the binaries d i,k,w as "step variables" with d i,k,w = 1 if and only if the outage beginning week for the cycle k of the unit i is before the week w. We extend the notations with d i,k,w = 0 for k > K,d i,0,w = 1 for w < 0, d i,k,w = 0 for w < 0 and k > 0. CT13 constraints, imposing that maintenance operations (i, k) begin between weeks To i,k and Ta i,k , reduce the definition of variables with d i,k,w = 0 for w < Ta i,k and d i,k,w = 1 for w Ta i,k .
The other variables have a continuous domain: refueling quantities r i,k for each outage (i, k), T2 power productions p i,k,t , fuel stocks at the beginning of campaign (i, k) (resp at the end) j,s,t , and fuel stock x end i,s at the end of the optimizing horizon. We note that T2 power productions p T2 i,k,s,t are duplicated for all cycle k to have a linear model, p T2 i,k,s,t = 0 if t is not included in the production cycle k. These variables are gathered in Table 4. Power production of T2 unit i at cycle k at t for scenario s, 0 if t is not in cycle k. p T1 j,s,t 0 Power production of T1 unit j at time step t for scenario s. r i,k 0 Refueling quantity for outage k of T2 unit i.
Fuel stock of T2 unit i at the end of the optimizing horizon at scenario s.
Fuel level at the beginning of production cycle k of T2 unit i .
x f in i,k,s 0 Fuel level before the refueling k + 1 of T2 unit i, after production cycle k.
Relaxing constraints CT6 and CT12, we have following MILP relaxation, denoted LB.
∀i, k, ∀i, k, s, ∀c ∈ C 14 , ∀w ∈ W, Constraints (2) are required with definition of variables d. Constraints (3) and (4) model CT13 time windows constraints: outage (i, k) is operated between weeks To i,k and Ta i,k . Constraints (5) model CT1 demand constraints. Constraints (6) model CT2 bounds on T1 production. Constraints (7) model CT3, CT4, and CT5 bounds on T2 production. Constraints (8) model CT7 refueling bounds, with a null refueling when outage i, k is not operated, i.e., d i,k,W = 0. Constraints (9) write CT8 initial fuel stock. Constraints (10) write CT9 fuel consumption constraints on stock variables of cycles k x init i,k,s , x f in i,k,s . Constraints (11) model CT10 fuel losses at refueling. Constraints (12) write CT11 bounds on fuel stock levels only on variables x init i,k,s which are the maximal stocks level over cycles k, thanks to (10). Constraints (13) model CT11 minimum fuel stock before refueling, these constraints are active for a cycle k only if the cycle is finished at the end of the optimizing horizon, i.e., if d i,k+1,W = 1, which enforces having disjunctive constraints where case d i,k+1,W = 0 implies trivial constraints thanks to (12). Constraints (14) are linearization constraints to enforce x end i,s to be the fuel stock at the end of the time horizon. x end i,s is indeed the x f in i,k,s such that d i,k,W = 1 and d i,k+1,W = 0, for the disjunctive constraints (14) that write trivial constraints in the other cases thanks to (12), defining S i = max k Smax i,k . Constraints (15)-(22) have a common shape as scheduling and resource constraints, as presented in [24]. The common format for constraints (15)-(19) may be seen as clique cuts, as written in [25,26]. Equivalently, using an additional and unitary fictive resource for each CT14-CT18 constraint that is consumed in the specific spacing/overlapping phases among outages, CT14-CT18 are written (15)- (19) as noticed by [37]. We note that the constraints (15) and (16) are not activated if Da i,k + S c < 0 for a given i, k, c, so that the generic equation is also valid using max(0, Da i,k + S c ) operations. CT19-CT21 are also resource constraints, with a non unitary quantity, which is written similarly in Equations (20)- (22).
LB is a dual bound for the original problem, relaxing only constraints CT6 and CT12, requiring less binary variables than the formulation of [37]. Each dual bound proven for LB is thus a dual bound for the EURO/ROADEF Challenge 2010. To face the resolution limits to calculate LB and its continuous relaxations, more relaxations are considered.

Parametric Surrogate Relaxations
To reduce the number of variables in the MILP relaxation LB, a parametric surrogate relaxation considers only the outages with an index k k 0 with the previous MILP formulation only for cycles k k 0 , and an aggregation of cycles k > k 0 in one production cycle without outages. This gives rise to the following MILP formulation where constraints M ordo ∀i, k k 0 , ∀i, s, ∀s, t, ∀i, k k 0 , s, ∀i, s, x end i,s As mentioned and proven in [17], this MILP formulation gives lower bounds for the EURO/ROADEF Challenge 2010: LB. Hence, each dual bound for LB(k 0 ) is a dual bound for the EURO/ROADEF Challenge 2010.

Preprocessing Reductions and Dual Bounds
This section aims to reduce the size of the MILP computation to provide dual bounds. The crucial point is to prove that the pre-processing reductions are exact processing operations, valid for optimal solutions of LB or LB(k) relaxations, or guarantee to provide lower bounds for LB or LB(k). Firstly, we mention the results proven in [17] that, under some conditions, time step aggregations applied to problems LB or LB(k) provide dual bounds for LB and LB(k). Therefore, we denote the following MILP as a general expression for LB or LB(k): Here, y t 0 denotes the time-indexed variables, i.e., the production variables p T2 i,k,s,t and p T1 j,s,t , whereas x denotes the vector of the other variables. A first important point is to notice that W 1 , W 2 do not depend on t in the MILP equations defining LB and LB(k). Two additional hypotheses are crucial.
Firstly, F t is constant over time, and we denote F = F t in that context. Secondly, we have to suppose that T1 production costs are constant over weeks, defining c w quantities: To define the time step aggregated version of (39), we define α = F/128, representing the proportion of a time period in a weekly period, and the aggregations T 2 Defining weekly production variables y w , which represent y w = ∑ t,w t =w αy t , we have the following aggregated version of MILP (39): We mention and prove now two propositions allowing for dealing with less variables in the previous MILP computations giving dual bounds for the EURO/ROADEF Challenge 2010. Propositions 3 and 4 are exact pre-processing results, allowing for deleting some variables in the MILP computations.
To i,k are first tightened by induction with k increasing, and then Ta i,k are computed using To i,k values.
is a first lower bound for the length of production cycle (i, k). Indeed, it considers a minimum refueling and a maximal fuel consumption: the maximal fuel consumption is given by a bounds of maximal power and the minimum fuel level after refueling is at least the minimal refueling Rmin i,k with the positivity constraint CT11. Denoting This minoration is valid for all feasible solutions, taking the lower bound of the LHS induces: is thus a valid pre-processing to tighten values of Ta i,k .
, we have following relations which allow Proof. We notice that Equations (10) and (11) induce a linear system of equalities considering as variables. A recursion formula can be deduced for x init i,k,s where the initial condition is x init i,0,s = Xi i : Lemma 1 (proven by induction in [23]) allows also to compute x init i,k,s as linear expressions of the variables r i,k and p i,k ,s,t for k < k. x f in i,k,s are also linear expressions of the variables r i,k and p i,k ,s,t for k < k. reporting last inequality in (10). Lemma 1. Let (u n ) n∈N be defined by induction with u n+1 = a n u n + b n , a n , b n ∈ R. We have:

Dual Bounds by Scenario Decomposition
In this section, we face another bottleneck for an efficient MILP solving: the number of scenarios. In the preliminary work [17], it was proven that a decomposition scenario by scenario allows for computing a valid dual bound for the original stochastic problem using S independent computations with the size of the deterministic problem. This section extends these results, so that a fixed maximal number of scenario can be chosen to group sub-computations. It is proven that generalized scenario decomposition improves the previous decomposition scheme from [17]. Numerical interest in the way of partitioning scenarios is also discussed in this section.
Here, x denotes the refueling and maintenance planning variables, i.e., d i,k,w , r i,k , whereas y s 0 gather the other continuous variables indexed by scenarios. Note that matrices A, W, and T do not depend from s. Let S 0 ⊂ S. Let v S 0 be the restriction of the MILP v sto to the subset of scenarios Considering several scenarios improves the scenario decomposition from [17].
Proof. ∑ N n=1 v S n v sto is proven using Lemma 2 with S 0 = S, and the partition (S n ), using v sto = v S and ∑ s∈S π s = 1. Using Lemma 2 for all n ∈ [[1; N]] with S 0 = S n and with the partition of S n into singletons S n = s∈S n {s}, we have ∑ s∈S n v {s} v S n . Summing these N inequalities, we have: In [17], only inequality allows for computing better dual bounds considering a restricted number of scenarios. The point is here to maximize the expected dual bound ∑ N n=1 v S n with an appropriate choice of the partition S n . Supposing that MILP computations of v S n are tractable for a given number of scenarios M > 1, we focus now on how to choose a partition S n with at most M scenario per subset, while trying to maximize the dual bound ∑ N n=1 v S n . Considering conjointly M identical scenarios, v S n computes the optimal maintenance and refueling planning for each scenario, it is the same optimal response for each scenario. Grouping the scenario in (48) induces the same optimal response, and the same dual bound, and it is an equality case in Lemma 2. Similarly, there may exist a very good solution x common for all the scenarios when similar scenarios are considered, implying a little gap in the inequalities from Lemma 2. To maximize the quality of dual bounds, it is thus preferable to have the most diversified scenarios. Hence, we partition the scenarios in subsets of size at most M, while maximizing the diversity among each subset of scenario.
Clustering problems maximizes the similarity inside clusters. To use clustering algorithms, with cardinality constrained versions, we should define a "distance" between scenarios s, s such that d s,s is higher than scenarios s, s are similar. Here, we prefer using a distance to measure the similarity of clusters, with d s,s = 0, and thus maximizing the diversity in the selected subsets of scenarios. To define such distance, we use the following measure: We use only the power demands to define the dissimilarity between scenarios, assuming it has the most influence in the costs of maintenance planning. The maximal number of scenario per cluster is denoted M > 0. Then, N = S M is the maximal number of clusters. We define binary variables x n,s ∈ {0, 1}, such that x n,s = 1 if and only if scenario s ∈ S is included in cluster n ∈ [[1; N]]. Maximizing the total distance intra-clusters, we have the following mathematical programming formulation of the maximal diversity in clusters: max ∑ n,s,s d s,s x n,s x n,s ∀n, ∑ s x n,s M ∀s, ∑ n x n,s = 1 (52) Objective (50) maximizes the total distance intra-clusters, using a quadratic function in variables x n,s . Contraints (51) bound the cardinal of each subset, to deal with at most M scenarios. Contraints (52) impose that each scenario is assigned to exactly one subset, to define a partition of the number of scenarios. To solve this last partitioning problem with an MILP solver, we use a standard linearization techniques for quadratic optimization problems with binary variables, introducing new binaries y n,s,s , such that y n,s,s = x n,s x n,s , to have a linear formulation in the space of variables x, y: ∑ n x n,s = 1 ∀n, s, s y n,s,s x n,s ∀n, s, s y n,s,s x n,s ∀n, s, s y n,s,s x n,s + x n,s − 1 Such formulation may be inefficient for a straightforward MILP solving with B&B solvers. Here, we need to have a good solution in a reasonable time. When the optimization problem (54) is not easily solvable with an MILP solver, one can design a matheuristic Hill-Climbing (HC) local search. Starting from an initial solution, one consider the optimization problem (54) in local re-optimizations with at most 2M scenarios and two clusters. One can operate the N(N − 1)/2 = O(N 2 ) optimizations considering all pairs of clusters to define one HC iteration. It is indeed an HC heuristic, each local optimization of (54) furnishing at least the same solution as the current one.
Algorithm 1 describes precisely such HC matheuristic. To define an initial partition of scenarios for the local search, one can consider subsets {1, .., M}, {M + 1, .., 2M}, . . . , {(N − 1)M + 1, .., S}, defining trivially a partition of S in at most N subsets containing at most M scenarios. A better initialization gathers scenarios using a one-dimensional, dispersion problem. Sorting firstly the scenarios of S such that the cumulated demands D s = 1 T ∑ t D s,t is increasing, i.e., s > s =⇒ D s D s , we can define an initial partition with S n = {s aN+n |a ∈ N, aN + n S}.
Algorithm 1 Partitioning HC matheuristic to maximize the diversity of scenarios. Designing and implementing Algorithm 1, some remarks are of interest:

•
Calling local optimization (54) with P n ∪ P n and M = 2, the current assignment is a feasible solution for the MILP (54), which is useful to speed up MILP solving.

•
The local optimization may be stopped before the optimality proof, only the improvement of a feasible solution in a given time limit is on interest. In such case, it is important to implement the warmstart previously mentioned, to ensure having a steepest descent heuristic: the local optimization will improve or keep the current solution given as warmstart.

•
In the loop re-optimizing (54) with P n ∪ P n with M = 2 for all n < n , many local reoptimizations are independent dealing with independent sub-sets, which allows for designing a parallel implementation.

•
One may wonder why the classical k-means algorithm is not used instead of Algorithm 1. A first reason is that cardinality constraints are not provided in the standard k-means algorithm.
One may try to repair the cardinality constraints after standard k-means iterations. Even in such scheme, k-means was not designed to partition the points in many subsets N with few elements, usually, k-means works fine with N M. On the contrary, the numerical property M N is favorable to an Algorithm 1 matheuristic. A second reason to prefer Algorithm 1 is that the objective function is a straightforward indicator of the data. Using a clustering algorithm, one may wish to minimize a dissimilarity. Considering dissimilarity indicators like 1/d − s, s , it is a deformation of the input data, and this is not a mathematical distance with a triangle inequality which also makes the clustering procedure non standard.
Lastly, we note that this section illustrates a reciprocal interest of hybridizing matheuristic and ML techniques. On one hand, ML partitioning techniques are used to improve the expected dual bounds with Proposition 5. On the other hand, a matheuristic methodology is used in Algorithm 1 to solve the specific partitioning problem for a non standard clustering application case.

Computational Results
This section presents the computational results. Before analyzing the quality of the different dual bounds, we report the characteristics of the datasets, the conditions for our computational experiments, and the solving capabilities for the different MILP formulations.

Computational Experiments and Parameters
MILP and LP problems are solved using Cplex version 12.9, using the ILOG CPLEX Optimization Studio (OPL) interface to model linear optimization problems. The OPL script is used to solve iteratively linear optimization problems. MILP computations of dual bounds are obtained with the OPL method getBestObjValue(). Note that, if the B&B algorithm terminates with the default parameter epgap = 10 −4 , the dual bound is given with a gap of 0.01% to the best primal solution found is given. The gap results are presented in our results with a granularity of 0.01%, so that we set parameter epgap = 10 −5 . Another stopping criterion for Cplex computations is to compute the LP relaxation, denoted LP in the results, or to stop the resolution when branching is necessary, denoted LP+cuts. LP+cuts computations show the influence of the cutting planes at the root node, using OPL parameter nodelim = 1 in the MILP solving mode. LP dual bounds are also provided using the MILP solving mode without cutting planes using parameters cutpass = −1 and nodelim = 1, which activates pre-processing for integer variables [39]. Using the LP mode of Cplex implies larger computation times than using the truncated MILP mode, the MILP pre-processing reduces considerably the number of variables for the LP relaxation. Using the integral pre-processing allows also to improve slightly the lower bounds given by the LP mode of Cplex, pre-processing variable fixing being like valid cuts.
Without any specific precision, the time limit for LP and MILP computations is M × 1800 s where M is the number of scenarios, so that the maximal overall computation time with Cplex is at most S × 1800 s. Computation times measures the sequential iterative computation of partial dual bounds. Scenario decomposition induces independent computations of MILP dual bounds and an easy parallelization in a distributed environment like Message Passing Interface (MPI). Dealing with S 0 S independent computations, the computations can be processed in a distributed way with MPI using S 0 machines with a single coordination operation at the end of the MILP computations: a reduction operation to sum the dual bounds given by each sub-problem. The clock time using distributed computations would be similar to the maximal solving time of the S 0 sub-problems.
Two different computers were used for the computational experiments, denoted Machine A and Machine B. Machine A is a quad-core desktop computer with processor Intel(R) Core(TM) i7-4790, 3.20 GHz, running Linux Lubuntu 18.04, with 16 GB of RAM memory. Without specific precision, tests are reported on Machine A. When larger computation power is required, Machine B is used. Machine B is a working station with a bi-processor Intel(R) Xeon(R) CPU E5-2650 v2, 2.60 GHz with 32 CPU cores, 64 GB of RAM memory, and running Linux Ubuntu 18.04.
Lastly, we note that the parameters for Algorithm 1 are M = 1 and M = S. In Algorithm 1, two iterations of the HC matheuristic are processed, ie nbIter = 2. Computing of dual bounds, the maximal number of considered scenarios is M = 10. It implies MILP computations of the optimization problem (54) with at most 20 scenarios. For such MILP computations, warmstarting was activated and the maximal solving time for Cplex was set to 10 s by precaution, optimal computations often converged much faster.

Data Characteristics
Three datasets were provided for the EURO/ROADEF Challenge 2010. These datasets are now non-confidential and available online [36]. The data characteristics are provided in Table 5. Dataset A is composed of five small instances given for the qualification phase. Production time steps are discretized daily for the instances of dataset A, in a horizon of five years, with 10 to 30 T2 units having six production cycles, and 10 to 30 stochastic scenarios. Instances B and X are more representative of real-world size instances. Production time steps are discretized with 8h time steps to analyze the impact of daily variability of power demands, in a horizon of five years, with 20 to 30 T1 units, around 50 T2 units, and 50 to 120 stochastic scenarios. Dataset X was secret for the challenge, it had to be similar to dataset B. This is not the case looking to the number of binaries, i.e., the cumulated amplitude of TW, given in in the column nbVar0 of Table 5. Instances B8 and B9 induce more binary variables, with very few TW constraints for cycles k 3. Actually, instances B8 and B9 are the most realistic for the real life application [24]. Table 5. Characteristics of the instances of the EURO/ROADEF Challenge 2010: values of I,J,K,S,T,W as defined in Table 1, nbVar0 and nbVar2 are respectively the total number of variables without and with pre-processing of Proposition 3, i.e., the total amplitude of the time windows, nbVar1 and nbVar3 are respectively the remaining variable after Cplex pre-processing with and without pre-processing of Proposition 3. Gaps show the improvements due to Proposition 3. To compare different lower bounds on an instance i, we use for each instance i the cost of the best primal solution known (BKS), denoted BKS(i). The BKS are mainly the one reported by [15] allowing 10 h computation time per instance, except for A1, A2, A3 given by [32] and for instance B6 where the BKS obtained for the Challenge had not been improved since, and is available at [28]. We compare lower bounds using the gap indicator:

Data
where v(i) denotes the considered lower bound for instance i.

Branch & Bound Solving Characteristics
We note that Table 5 analyzes also the impact of pre-processing strategies. The specific pre-processing of Proposition 3 is not redundant with the generic MILP pre-processing of Cplex [39], and the gain in MILP reduction is significant, especially for the difficult instances B7, B8, and B9. Pre-processing of Proposition 3 has a positive impact influence in our computations of complex MILP models to reach the most advanced phases of the B&B algorithm with memory limitations. Proposition 4 reduces the number of continuous variables that will be decisive, reaching memory limits. In the following analyzes, pre-processing of Propositions 3 and 4 are always enabled.
With Cplex 12.5, dual bounds can be computed for each instance with a single scenarios and weekly production time steps in less than one hour for the LP and LP+cuts modes, the B&B convergence characteristics are detailed using Cplex 12.5 in [17], and reactualized results with Cplex 12.8 are available in [24]. We note that 1 hour resolution time was not enough to solve the LP relaxations of such problems for instances B8 and B9 with Cplex 12.3 [23]. Computations LB(k) are quicker than LB computations, reaching more advanced phases of the B&B convergence.
The number of continuous variables is decisive for solving capabilities. Indeed, there are the same number of binary variables with disaggregated or aggregated time steps for production variables, and the computation capabilities differ dramatically. If aggregated computations on a single scenario are tractable, disaggregated computations are tractable only for dataset A, and instances from datasets B and X induce a prohibitive memory consumption. Similarly, considering several stochastic scenarios does not change the number of binary variables, but limits are reached in terms of memory resources. This seems paradoxical, exponential complexity for MILP comes with integer variables. In this problem, increasing the number of continuous variables makes polynomial algorithms intractable. We note furthermore that this occurs also for the MILP relaxation with lighter CT6 constraints presented in [24]. Such formulation considers the same number of binary variables, but the MILP solving capability is also not efficient for the instance sizes of the Challenge. This highlights the interest of Proposition 4, in reducing the number of continuous variables to increase the MILP solving capability. With the recent versions of Cplex, dual bounds can be computed with truncated MILP solving considering up to 10 scenarios for LB(k) and LB problems, with aggregated time steps for all the instances, and with disaggregated for the instances from the dataset A.
Lastly, we note that the variable definition has an influence in the efficiency of the branching quality. Variables d i,k,w d i,k,w+1 imply better branching and B&B convergence than using binaries x i,k,w = d i,k,w − d i,k,w−1 as in [26] with constraints ∑ w x i,k,w 1. These results are coherent with [20], and the standard branching on variables implemented in Cplex induces for d i,k,w variables well balanced B&B trees. However, using variables x i,k,w was slightly better for the search of primal solutions. To design a primal heuristic in [26], and using few branching, variables x are relevant for the primal heuristic in [26], whereas, in our application, only the quality of dual bounds matters and choosing step variables d i,k,w is more efficient.

Lower Bounds for Dataset A
For none of the instance of dataset A, the hypothesis (40) holds to guarantee dual bounds with time step aggregation and Proposition 2. However, computations with disaggregated time steps are tractable with few stochastic scenarios, and dual bounds are computable using Proposition 5. Table 6 compares the dual bounds for the EURO/ROADEF Challenge 2010, computing dual bounds of LB with LP and LP+cut relaxations, and MILP computations truncated in 30 min, to the bounds obtained with surrogate relaxations LB(k) with k = 4 and k = 5. Table 6. Dual bounds for the dataset A, comparison of former dual bounds of the literature with the dual bounds LB(k 0 ) with k 0 > 3 and LB without aggregation of production time steps and scenario decomposition. M indicates the maximal number of scenarios in sub-problems. MILP computation of dual bounds are stopped after 2400 × M seconds using Machine B. Bold values denote stopped computations with optimal dual bounds in the 0.001% tolerance for MILP computations. LP and LP + cuts are also the terminating values, no additional gap is implied by the time limit.  In Table 6, it is preferable to compute LB bounds instead of LB(k) bounds for the dataset A. Indeed, dual bounds are computable with LB bounds, avoiding the additional gaps to the primal solutions with approximations LB(k). For instances A1, A2 and A3, LB bounds are computable to MILP optimality with up to 10 scenarios. It induces lower bounds of a very good quality for these instances. For instances A4 and A5, the computation times do not allow for converging to optimality with MILP computations of LB in the defined time limits. Furthermore, the memory space is limiting, and computations with 3600 s per scenario are too much memory-consuming. We note that, for the A5 instance, it is preferable to optimize conjointly by groups of 3 instead of 5. This is due to the stopped MILP computation, and it is preferable to reach more advanced phases of the B&B algorithm with three scenarios rather than using more scenarios and stopping earlier than the B&B algorithm. Otherwise, increasing the number of scenarios increases the quality of the dual bounds in Table 6, and some relations were proven thanks to lemma 2 with optimal computations. This holds also empirically for the stopped computations on instance A4.
Computation times are not reported in Table 6. Using Machine A, 6395 s (i.e., more than 1 h 45 min) are measured to solve to optimality LB with the 10 scenarios of instance A1, but the time printed by Cplex is only 1014 s (less than 20 min). With three scenarios, the total time measuring the iterated MILP computations using Machine A, compared to the sum of the Cplex solving time are: • for instance A1: we measure 1191 s (nearly 20 min) for the total MILP computations, and summing Cplex time we find nearly 135 s (2 min 15 s). • for instance A2: we measure 10,385 s (nearly three hours) for the total MILP computations, and summing Cplex time we find 2100 s (i.e., 35 min). • for instance A3: we measure 11,835 (more than three hours) for the total MILP computations, and summing Cplex time we find 2418 s (i.e., nearly 40 min).
The differences are not due to Algorithm 1, very quick and not counted in these reported time measures. The differences correspond to the loading time of MILP matrix, before the MILP pre-processing. We note that using Machine B instead of Machine A for the results in Table 6 solves this problem. With Machine B, most of the computation time is devoted to the B&B solving, and loading the MILP models was around 10% of the total computational time. To explain the difficulty met with Machine A, huge matrices are stored by Cplex, before removing many variables and constraints with MILP pre-processing. Using OPL and having to define variables with rectangular multi-indexes arrays, all the variables d i,k,w and p T2 i,k,s,t are initialized firstly before pre-processing operations detect that constraints (3), (4) and (7) allow for deleting most of the variables d i,k,w and p T2 i,k,s,t . Production variables p T2 i,k,s,t have the most impact on these difficulties. One way to reduce the MILP loading time would be an implementation using Cplex API to define only the required variables d i,k,w and p T2 i,k,s,t taking into account constraints (3), (4) and (7), so that Cplex will load a smaller matrix before the MILP solving.
In Table 6, the improvements with partial scenario decomposition are significant, but not so important to close the gap to the optimal dual bounds. For instance A1, we have the optimal values of best and the worst lower bounds of LB with scenario decomposition from Proposition 5: any scenario decomposition implies a gap to the BKS between 0.04%, the worst bound with M = 1, and 0.01%, the best bound with M = 10. The little impact of scenario decomposition is explained with Figure 1, showing the demand curves for each scenario in instance A5. The stochastic scenarios are distributed with a tube structure related to the meteorological seasonality, which is highly correlated to the power demands. It explains that few improvements are observed, the scenarios are relatively close, and induce optimal maintenance dates in similar zones of the time horizon. One may find, for any subset of scenarios, a good compromise in the dates of maintenance operations for each considered scenario. It induces also that choosing badly partitions of scenarios in Proposition 5 may induce very few differences with the scenario decomposition with M = 1. The significant gap improvements show the interest of a careful partitioning procedure with Algorithm 1.
Power demands (power unit is unknown and obfuscated in the Challenge datasets) time steps

Lower Bounds for Datasets B and X
The hypothesis (40) is valid for datasets B and Proposition 2 ensures that the aggregation of production time steps provides dual bounds of LB for datasets B and X. In this section, we always aggregate production time steps, comparing the impact of choosing LB or LB(k) formulations and analyzing the impact of the scenario decomposition. The huge loading times required by Cplex mentioned for the instances of the dataset A is not so significant for the datasets B and X. It is due to the aggregation of production time steps which reduces drastically the number of continuous variables.

Dual Bounds Computable in Less Than One Hour
In this section, we focus on dual bounds that can be obtained in less than one hour for datasets B and X, as in [16]. One hour is reasonable for the real-life application, we note that computing separately dual bounds from a primal heuristic is helpful to detect that primal solutions have a prohibitive cost.
Using MILP formulations LB(k 0 ) with k 0 2, one obtain our quickest computations of dual bounds with MILP sub-problems dealing with a single scenario (M = 1). Such bounds and the computation times are reported in Tables 7 and 8 with respectively LP+cuts and MILP computations. With k 0 = 0 and M = 1, LB(k 0 ) problems are solved to optimality as MILPs, and the solving time is, on average, around 200 s for the total measured time. Such lower bounds have a comparable quality to the best bounds from [16]. However, the computation times are not reported precisely in [16]. It is only mentioned that solving the flow-network for all scenarios takes up to an hour for the biggest instances. The quickest lower bounds from [16] are reported with the column [16].1 in Tables 7 and 9 and have a poor quality, but the solving times are computed in less than 10 s. Our approaches cannot provide dual bounds in such restricted times with sequential implementations, using MILP distributed computations with MPI, the dual bounds in Tables 7 and 8 with LB(0), and M = 1 can be computed in less than 20 s.  Table 7 shows that LB(1) lower bounds with M = 1 and MILP computations outclass significantly the dual bounds from [16], requiring less than one hour computations. Tables 7 and 8 illustrate that computation times differ among the different instances. B8 and B9 are also the most difficult instances for the computations of LB(k 0 ) with k 0 3. B6,B7, X10, and X11 instances have the quickest computation time, and this is due to their minimal number of stochastic scenarios, as shown in Table 5. Generally, the LP+cuts computations of root node of the B&B tree, as in Table 9, induce a little degradation of the lower bounds with full MILP computations presented in Table 7, for a significant speed up comparing the computational times in both result tables. In a defined time limit, it is more efficient to increase the value k than reaching more advanced phases in the B&B algorithm.
This induces results of Table 9, considering in our computational experiments only the results having a total computational time inferior to one hour, which corresponds to some LP and LP+cuts computations. We note that we did not optimize the results obtained in one hour. Table 9. Comparison of our best dual bounds for datasets B and X allowing at most 1 h total computation to the ones from [16]. The results of Tables 7-9 show a graduation of lower bounds with an increasing quality when computation times increase. Having the computation times of the lower bounds [16] or other dual bounds, one can compare the quality of lower bounds in equivalent computational times using equivalent hardware. Our dual bounds are very scalable, and, if a total computation time is allowed, one can divide this total solving time by the number of clusters of scenarios, in order to define a time limit for each MILP sub-problem. In particular, the lower bounds for instances B8 and B9 should be significantly improved allowing one hour total computations and LP+cuts truncated computations.

Dual Bounds Decomposing Scenario by Scenario
Computations of dual bounds LB are tractable for each instance from datasets B and X with LP computations or truncated MILP computations for one scenario and time step aggregation. Table 10 compares the dual bounds of [16] with one resulting from LB(4), LB (5), and LB formulations with LP, LP+cuts, and MILP truncated computations. We note that, contrary to computations in the average scenario presented in [17,24], the B&B convergence of LB computations was less advanced for each scenario in the computations of Table 10, and the average scenario seems to be favorable for the efficiency of the B&B search with Cplex. In addition to Tables 7 and 9, Table 10 presents the gaps obtained using formulations LB(k 0 ) with k 0 varying. Table 10 extends some empirical results mentioned for Tables 7 and 9. The LP+cuts computations of root node of the B&B tree induce a few degradation of the lower bounds compared with the full MILP computations of the same model, for a significant speed up in the computational times when computations are stopped at the root node. In a defined time limit, it is more efficient to increase the value k to compute LB(k) bounds. Computation times are not reported in Table 10, and the MILP computations reach the defined time limits, 1800 s by sub-problem. Such quality of dual bounds can thus be obtained in about 1800 s in clock time, distributing the S computations of sub-problems in a cluster of S computers with MPI distributed implementation.
In Table 10, it is interesting to compare the columns related to LB(5) and LB dual bounds. For most of the instances, computing LB dual bounds induces slightly better bounds, LB outclasses LB(5) computations by more than 0.5% only for X11 and X13 instances. Contrary to Table 6, the lower bounds obtained with LB(5) dual heuristics may be better than the exact formulation. This is the case for most difficult instances B7, B8, and B9 as highlighted in Table 5. This situation seems paradoxical, and we know that LB(5) LB with Proposition 1. Figure 2 illustrates these situations.
For the easier instance B6, LB lower bounds are always better than LB(5) lower bounds. This situation will continue until the curves converge to the optimal values of LB and LB (5). For the more difficult instance B7, the LB(5) lower bounds are always better than LB lower bounds, which seems paradoxical as in the convergence state we have LB(5) LB with Proposition 1. Considering reasonable computational times for our application, it is more efficient to use LB(5) formulation to compute partial lower bounds or the most difficult instances. For many instances, LB(5) dual heuristics are better than LB with the LP+cuts stopping criterion, as shown in Table 10. Having smaller MILPs with lB(5), more efficient cuts are generated, and it is decisive for the dual bound quality with LP+cuts computations. Good primal solutions generally plan four or five outages for each T2 unit. Hence, a 6th outage induces more costs as T2 production is cheaper than T1 production. It explains that LB(5) dual bounds are good approximations, and the major approximation is due to the approximation of the penalization costs fore the remaining fuel. The possibility to have a 6th outage induces more binaries and continuous variables, which is highly penalizing for the cutting planes generation and MILP solving capacities of LB dual bounds. For instances B10, X12, and X14, the LP+cuts dual bounds are better with LB(5) and allowing branching allows for having better dual bounds with LB relaxation.
This illustrates the concept of dual heuristics: relaxations are parametrized to lead to few degradation quality of the dual bounds, but improving significantly the solving capacities. Our first motivation was to have dual bounds for larger sizes of instances than using only exact formulations, and to have a better scalability in computation times with the parametric surrogate relaxations. The notable point here is that dual heuristics LB(5) may also improve the dual bounds LB, with approximation helping the MILP solvers to focus on a specific crucial part of the problem and relaxing a part of constraints that has little impact on the value of the objective function. This approach is similar to primal heuristics when some constraints have little influence on design solutions of good quality. In such applications, the expert knowledge of the problem and numerical properties are crucial to identify the most important constraints, variables, and difficulties of the problem.  Table 11 compares the dual bounds for datasets B and X using LB(5) and LB formulations with LP+cuts computations varying M the maximal number of scenarios in each sub-problem.
Like in Table 6, the impact of partial scenario decomposition is significant in Table 11, but not so important to close the gap to the optimal dual bounds. The distribution of stochastic scenarios illustrated in Figure 1 holds also for datasets B and X, and this is related to the industrial application. Badly partitioning the scenarios in Proposition 5 also induces few differences from the scenario decomposition with M = 1. The gap improvements with M increasing are significant, and it shows the interest of Algorithm 1. This section provides our best results, using truncated MILP computations instead of LP+cuts stopping criterion for the most promising schemes. Table 12 presents our best lower bounds and compares it to the former ones [16,17]. The parameters leading to the reported best dual bounds are presented in Table 12: the MILP relaxation LB or LB(5), the maximal number of scenario in an MILP sub-problem and the maximal time for an MILP computation, time/S denotes the CPU time in seconds per scenario, i.e., an MILP computation considering m scenarios will be limited to time/S ×m seconds. In Table 12, OPT in the time/S column denotes that the dual bounds are computed optimally in the 0.001% tolerance gap set for MILP solving with Cplex.
The best results for instances A4 and A5 were already given in Table 6. For these instances, the memory space is a limiting factor to reach better dual bounds. For instances A1, A2, and A3, we can compute optimally the LB bounds with 10 scenarios with an accuracy of 0.001%. For instance A1, 10 scenarios correspond to the total number of scenario so that Table 6 gives the best lower bound that we can compute with the LB relaxation, which is also the dual bound that would give the Bender's decomposition from [26]. Our new dual bounds for the datasets B and X improve by more than 0.50% the older ones from the preliminary work [17], which highlights the interest of the partial scenario decomposition and the ML-guided partitioning procedure in Algorithm 1.
The most difficult MILP computations are obtained for instances B8 and B9, increasing for instances B8 and B9 the time/S parameter to 3600 s, which was the solving time for the results [17], Each sub-problem reaches the branching phase which was not the case with time/S = 1800 s. Finishing the cutting planes passes has a significant impact, as already observed in Table 10, but also the first branching operations improve significantly the lower bounds from the root node of the B&B tree in general. This explains the improvement of the lower bounds for instances B8 and B9, from 11.95% and 11.25% in the truncated LP+cuts mode in Table 11, to the final gaps 11.69% and 10.87%. For instances B6, B7, X11, and X12, setting time/S = 3600 s implies an improvement about 0.1%, reaching more advanced phases of the B&B algorithm, whereas the parametrizing time/S = 1800 s allows to reach more than 10,000 nodes in the B&B tree. For the instances of datasets B and X, the memory space is not a limiting factor, and better dual bounds can be obtained allowing more time for the computations of MILP sub-problems.  [16,17], and presenting the configuration leading to our best lower bounds. Gaps are calculated for the best primal solution known from [15,28,32]. Machine A is used for instances of datasets B and X, whereas Machine B was used for the instances of datasets A. of dual heuristics for an industrial MILP: relaxations are parametrized to induce little degradation of the dual bound's quality, but it significantly improves the MILP solving capabilities, so that in defined time limits dual heuristics may provide better dual bounds than exact approaches. The scenario decomposition computes dual bounds for the original problem with a fixed number of scenarios, with a numerical interest to choose carefully the partitions of scenarios, using ML techniques. Combining these techniques leads to tractable computations of dual bounds for the EURO/ROADEF Challenge 2010, outclassing significantly the former best dual bounds of the literature. The quality of our dual bounds justifies the quality of the best primal solutions known, mainly in [15]. This work offers new perspectives. Using more powerful computers and more computation time may improve some reported results. Our approach is scalable to compute dual bounds in a given time limit and hardware configuration, which is useful for a real life application. We note special perspectives for the approach deployed by [26]: parametric MILP relaxations may be useful to accelerate Bender's decomposition, to compute quicker valid cuts and also to stabilize the Bender's decomposition. Considering bi-objective extension, implied by robust optimization considerations [38] or stability costs [24], we can compute dual bound-sets using computations of dual bounds from this paper using scalarization or epsilon-constraint methods. Lastly, we mention as a perspective that the results of Section 5 are valid for a general class of 2-stage stochastic problems, optimizing strategical problems like maintenance planning conjointly with an operational level as in [3].

Conflicts of Interest:
The authors declare no conflict of interest.