Matheuristics to optimize maintenance scheduling and refueling of nuclear power plants

by   Nicolas Dupin, et al.

Scheduling the maintenances of nuclear power plants is a complex optimization problem, formulated in 2-stage stochastic programming for the EURO/ROADEF 2010 challenge. The first level optimizes the maintenance dates and refueling decisions. The second level optimizes the production to fulfill the power demands and to ensure feasibility and costs of the first stage decisions. This paper solves a deterministic version of the problem, studying Mixed Integer Programming (MIP) formulations and matheuristics. Relaxing only two sets of constraints of the ROADEF challenge, a MIP formulation can be written using only binary variables for the maintenance dates. The MIP formulations are used to design constructive matheuristics and a Variable Neighborhood Descent (VND) local search. These matheuristics produce very high quality solutions. Some intermediate results explains results of the Challenge: the relaxation of constraints CT6 are justified and neighborhood analyses with MIP-VND justifies the choice of neighborhoods to implement for the problem. Lastly, an extension with stability costs for monthly reoptimization is considered, with efficient bi-objective matheuristics.



There are no comments yet.


page 1

page 2

page 3

page 4


Mixed integer formulations using natural variables for single machine scheduling around a common due date

While almost all existing works which optimally solve just-in-time sched...

Column generation for the discrete Unit Commitment problem with min-stop ramping constraints

The discrete unit commitment problem with min-stop ramping constraints o...

Dominance inequalities for scheduling around an unrestrictive common due date

The problem considered in this work consists in scheduling a set of task...

Exploring Mixed Integer Programming Reformulations for Virtual Machine Placement with Disk Anti-Colocation Constraints

One of the important problems for datacenter resource management is to p...

Maintenance scheduling of manufacturing systems based on optimal price of the network

Goods can exhibit positive externalities impacting decisions of customer...

The Maintenance Location Choice Problem for Railway Rolling Stock

Due to increasing railway use, the capacity at railway yards and mainten...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The 2010 ROADEF/EURO Challenge was specified by the French utility company (EDF) in [25] to address the complex industrial optimization problem to schedule nuclear power plant outages to process maintenances and refueling. Scheduling the maintenances of nuclear power plants is a complex industrial problem: the impact of the outages is calculated with the expected production costs to fulfill the power demands using the available power plants. Furthermore, the feasibility of a high-level maintenance planing must be ensured for low-level production constraints. The optimization problem was formulated using 2-stage stochastic programming for the challenge. Power demands, production capacities and costs were stochastic, discrete scenarios model this uncertainty. The first stage optimization problem concerns the maintenance decisions for nuclear power plants and the refueling quantities. The second stage optimization problem computes production plans implied by the first stage decisions for all stochastic scenarios, to minimize the average production cost.

This ROADEF Challenge is a complex industrial optimization problem with varied constraints and types of variables. The best results were mainly obtained in [13] with an aggressive local search approach. Several approaches in competition were based on exact methods like ([18, 22, 27]

), this required heuristic reductions, constraint relaxations and repairing procedure to compute final solutions. The large size of the instances was a bottleneck for the exact methods. With the short development times of the Challenge, the competition focused on having efficient primal heuristics for the whole problem. Intermediate questions were opened for a better understanding of the problem structure. In this paper, we focus on a deterministic version of the problem with few relaxations/simplifications. This work is both motivated to design an efficient algorithm for an industrial application, and to understand better some structures of the problem and some results of the ROADEF Challenge.

This paper is organized as follows. In section 2, we give an overview of the problem. In section 3, we discuss related state-of-the-art elements to appreciate our contributions. In section 4, we present a new MIP compact formulation. In section 5, constructive matheuristics are derived from the previous MIP formulation. In section 6, a Variable Neighborhood Descent with MIP neighborhoods provides local improvements once feasible solutions are known. In section 7, the computational results are discussed on real world instances. In section 8, our contributions are summarized, discussing also future directions of research.

Sets and indexes

Weekly time steps.
Flexible (Type 1, T1) power plants.
Nuclear power plants (Type 2, T2).
Cycles related to T2 units, for initial conditions.
Points of the CT6 profile for the cycle .

Temporal notations

Power demands at time step .
Conversion factor between power and fuel.

Notations for T1 units

Production Costs proportional to the generated power at .
Minimal power to generate at time step .
Maximal generated power at time step .

Notations for T2 units

Maximal fuel level remaining to process outage .
Fuel level ”Bore O” of cycle
Proportional cost to the final remaining fuel levels at .
Loss coefficient of generated power at mode and cycle .
Additional (stability) cost to schedule outage at week .
Initial week where outage was scheduled
Outage duration for maintenance and refueling at cycle .
First possible outage week for cycle of T2 plant .
Last possible beginning week for outage of T2 plant .
Fuel levels for the stretch decreasing profile at mode and cycle .
Maximal generated power at time step .
Proportion of fuel that can be kept during reload in cycle at plant
Minimal refueling at outage .
Maximal refueling at outage .
Maximal fuel level of T2 plant at production cycle .
Initial fuel stock of T2 plant .
Table 1:

Definitions and notations for the imput parameters

2 Problem statement

We consider in this paper some simplifications of the constraint of the 2010 ROADEF/EURO Challenge, as specified in [25], and an extension to optimize stability of monthly reoptimization. In the Challenge, there were two types of time steps, weekly discretization for the outage dates and a hourly discretization for production time steps. In this study, we consider only weekly time steps, production and demand levels are aggregated with their average value in each week. In the Challenge, discrete scenarios modeled the uncertainty of power demands, production capacities and costs. In this study, we consider only one scenario, aggregating the uncertain parameters with their average value over all the scenarios of the Challenge. Lastly, we relax in this study constraints CT6 (”stretch” decreasing profile) and CT12 (”modulation” constraints).

2.1 Set and index

Two kinds of power plants are modeled. On one hand, Type-2 (or T2) power plants indexed with , correspond to nuclear power plants. T2 power plants have to be shut down for refueling and maintenance regularly. On the other hand, Type-1 (or T1) power plants are indexed with , model other power plants with more flexibility in the production. Outages and production campaigns are indexed with the cycles for all T2 plants . By convention, a cycle begins with the outage period for maintenance and refueling, before the production campaign. The outage and production decisions are discretized weekly and indexed with .

2.2 Objective function

The objective function of the Challenge minimizes the power generation cost while satisfying the customer load for all time steps. Production costs of T1 units are proportional to the production levels at a given time step . Production costs of T2 units are calculated proportionally to the fuel consumption. This makes an aggregated cost function with weights defined to optimize the global financial cost.

In this study, we extend the objective function to consider stability costs for the operational application. Indeed, the French Utility Company reoptimizes monthly the planning of maintenances. The baseline maintenance planning impacts the reoptimized one. The maintenance operations are outsourced and formalized/organized by contracts to external companies fixing the periods of operations. In the model, time window constraints CT13 can model such decisions with hard constraints. In the reality, maintenances can be slightly modified if the financial stakes worth the reorganization efforts, from a contractual and organizational point of view. To model this fact, we consider an objective of stability, defining with the stability cost to reschedule the maintenance at week . Denoting the initial beginning week of maintenance , . Any penalization with can be considered A linear penalization can be written with , or a quadratic function are natural choices. Another penalization function can consider constant costs for modified maintenances, with for all . We note that the case of the Challenge is extended with .

CT1, demand covering: for all time step , the total production of T1 and T2 power plants must equalize the demands .
CT2, T1 production bounds for all , the production domain of T1 plant describes exactly the continuous domain
CT3-5, T2 production bounds The productions of offline T2 power plants are null during outages. Otherwise, the production domain of T1 plant describes exactly the continuous domain for all period .
CT6, ”stretch” constraints (relaxed fully or partially in this study) During every time step of the production campaign of cycle , if the current fuel stock of plant is inferior to the level , the production of is deterministic, following a decreasing profile, piecewise linear function of the stock level, as in Figure 2.
CT7, refueling quantities: The refueling possibilities for outage of T2 plant describes exactly the continuous domain
CT8, initial fuel stock: The initial fuel stock for T2 unit is , known and common for all scenario .
CT9, fuel stock variation during a production campaign: The fuel stock variation during a production campaign of a cycle between and , is proportional to the power produced by at time step , with a proportional factor .
CT10, fuel stock variation during an outage During an outage, the fuel stock variation is the sum of the decisional refueling bounded with CT7 and a certain amount of unspent fuel, calculated with a proportional loss with factor to the residual fuel before refueling.
CT11, bounds on fuel stock The fuel level is in for cycle of T2 unit . The fuel level must be lower than to process outage .
CT13, time windows for maintenances The beginning dates of outage of T2 unit must be included in . Implicit constraints: The maintenances follow the order of set without skipping maintenances: if outage is processed, it must follow the cycle .
CT14-21, resource constraints in maintenance scheduling: For all constraint , a subset of outages (possibly ) and a subset of time periods (possibly ) is considered. For all outage and for all time period , a resource consumption is defined for the outage being scheduled at , and the global resource comsumption at must be lower than
Table 2: Constraints coupling productions, fuel stock levels and outage decisions

2.3 Constraints description

Table 2 defines the constraints with their nomenclature from CT1 to CT21 in the challenge specification [25]. CT1 are the demand constraints to equalize productions and demands for each time step and each scenario. Constraints CT2 to CT6 and CT12 model production constraints: production bounds for T1 and T2 power plants and specific technical constraints of nuclear power plants. Constraints CT7 to CT11 model stock level constraints: bounds on stock levels and refueling and CT9 is the equation linking T2 production and fuel consumption. The remaining constraints (CT13 to CT21 ) are specific to T2 plants outage scheduling. CT13 defines time window constraints to begin the maintenances. Constraints CT14 to CT21 can be unified in a common format of resource constraints applying for the maintenance decisions, as defined in Table 2. These constraints express minimal spacing/ maximal overlapping constraints among outages, and Minimum spacing constraint between coupling/decoupling dates, resource constraints for maintenances using specific tools or manpower skills and limitations of simultaneous outages with maximal number of simultaneous outages or maximal T2 power off-line.

2.4 Justification of the relaxation of constraints

In this section, we justify the reasons to relax fully modulation constraints CT12, and the stakes to consider partial decreasing profile constraints CT6.

2.4.1 Modulation constraints CT12

In this paper, modulation constraints, denoted CT6 in the Challenge description, are fully relaxed. In the Challenge, it is formulated as a maximal cumulated volume of energy generated when a unit is not producing at the maximal power at time step , the cumulated energy beeing counted for all cycle when the fuel level is superior to . When the fuel level is inferior to , the production is at the maximal power defined by the CT6 decreasing profile.

A first reason to relax constraints CT12 (especially in [22]) is that modeling CT12 induces more difficulties, requiring additional binary variables for a MIP formulation. A second reason is the evaluation of the impact of such constraints. T2 production is indeed less expensive than T1 production, the cost optimization induces that T2 units are likely to produce at the maximal powers when they are online, so that the modulation capability is a marginally used. The case where modulations cannot be avoided is when the total power of online T2 units exceed the global demands. This phenomenon is local and predictible knowing the demand curves, it suggests to relax firstly constraints modulation CT12 before to repair locally the violated modulation constraints to tackle the problem of the Challenge.

A structural reason to relax constraints CT12 comes from the real technical constraints of nuclear power plants. The real constraints are that the duration where nuclear units do not produce at the maximal power is bounded, with a maximal duration in a intermediate power level similarly with [9]. This induced for tha large time periods of the Challenge to limit the global volume of modulations in cycles. Actually, it is possible to shut down completely a nuclear unit in a production cycle, without maximal duration on the offline period, with minimal durations for the offline and online periods, similarly with [26].

2.4.2 Decreasing profile for T2 units CT6

In the Challenge, the “stretch decreasing profile” CT6 are defined when the fuel stock level of a T2 unit is inferior to the level in the cycle . The production of is deterministically defined, following a decreasing profile, piecewise linear function of the stock level, as in Figure 2.

Figure 1: Illustration of the production domain for T2 power plants in the Challenge ROADEF
Figure 2: Illustration of the production domain for T2 power plants with light CT6 constraints
Figure 1: Illustration of the production domain for T2 power plants in the Challenge ROADEF

In this paper, we will fully relax constraints CT6 or we will consider only the upper bounds constraints as illustrated in Figure 2. A first reason to relax constraints CT6 is is that modeling CT6 induces more difficulties in the MIP modeling and in the MIP solving, requiring additional binary variables for a MIP formulation. A second reason is due to the impact of considering such constraints: T2 production has lower marginal costs than the T1 production. Hence, the optimization tends to have T2 units producing at their maximal power. It justifies the first relaxation with production domain as illustrated in Figure 2. The full relaxation of CT6 seems a priori interesting, to study if complete solutions can be repaired starting from solution of the model relaxing fully CT6 constraints. One of the applicative issue of this work is to analyze the impact for primal heuristics of such modeling variants in the real life solutions.

3 Related work

This section describes the state of the art of the solving approaches for the 2010 EURO/ROADEF Challenge, focusing on the related works both in the heuristic and mathematical programming approaches.

3.1 General facts

Although the teams of the challenge worked independently, similar ideas occurred. It seems strongly related to the specificity of the problem and the size of instances. A major difficulty comes with the large size of instances. Many approaches reduced the size of instances with different preprocessing strategies: reducing the time windows with exact or heuristic preprocessing, or aggregating production time steps to week or aggregating the stochastic scenarios to the average one. The heuristic preprocessing requires a repairing postprocessing, which made no difficulty to recover feasible solutions for the whole problem.

We distinguish three main types of approaches for the ROADEF challenge, matheuristics based on MIP exact approaches, decomposition heuristics following the 2-stage structure and straightforward local search approaches. We note only one attempt of population heuristics, with an ant colony optimization approach, with reported results using more memory than specified. An explanation comes from the memory limitations in [25] and the size of instances.

An open question after the Challenge was to prove dual bounds for the large size instances of the competition. Brandt et al. (2013) furnished the first dual bounds [7] using dual heuristics, the bounds proven in [11, 10] improved these last bounds.

3.2 Matheuristics based on MIP exact approaches

Two matheuristic approaches were designed for the ROADEF Challenge. A simplified MIP problem is solved, and the solution is repaired to ensure the feasibility for all the constraints and computing the cost in the full model. In both cases, the MIP model is solved heuristically using truncated exact methods, the simplified MIP being too difficult to be solved in one hour. We note that Joncour et al (2010) provided an exact and compact MIP formulation for the full problem, introducing binaries for all time steps, cycle and production mode [17]. This size of problem is not reasonable for a B&B search, even in a truncated mode, this work was a preliminary work before a Column Generation (CG) work similar with [27].

Rozenkopf et al. (2013) considered an exact formulation of CT6 constraints, in a CG approach dualizing coupling constraints among units [27]

, i.e. CT1 demands and CT14 to CT21 scheduling constraints. The MIP considered a unique scenario, the average one, with production time steps aggregated weekly. The CG approach is deployed to compute a LP relaxation relaxing scheduling constraints CT14-CT21, with subproblems solved by dynamic programming. The further CG heuristic incorporate these scheduling constraints in the integer resolution with the columns generated. The solution of this MIP gives the outage dates and week where the CT6 constraints are activated, the final cost and feasibility issues for the whole problem are given computing the production by Linear Programming (LP). This approach was one of the most effective for the challenge, up to

of the best solution known.

Lusby et al. (2013) relaxed fully the constraints CT6 and CT12, leading to a MIP formulation with binaries only for the outage decisions [22]. The production time steps were aggregated to weeks for size reasons. The stochastic scenarios were not aggregated, leading to 2-stage stochastic programming structure solved by Bender’s decomposition. The master problem concerns the dates of outages and the refueling quantities. Independent sub-problems are defined for each stochastic scenarios with continuous variables for productions and fuel levels. The heuristic of [22] computes first the LP relaxation exactly with the Bender’s decomposition algorithm. Then, a cut&branch approach repairs integrity, branching on binary variables without adding new Bender’s cuts. The resulting heuristic approach was efficient for the small dataset of the qualification, difficulties and inefficiencies occur for the final instances of the competition.

3.3 Heuristics based on the 2-stage decomposition

A natural idea to solve the problem by decomposition is to follow the 2-stage structure of stochastic programming, distinguishing the high-level maintenance and refueling problem of T2 units from the lower level production problems, as in [2, 6, 7, 14, 15, 18].

The high-level problem fixes the maintenance dates and the refueling levels modifying slightly the current solution and respecting scheduling constraints CT7-CT11 and CT13-CT21. The previous maintenance planning of the last iteration (or the initial planning) is slightly modified in a defined neighborhood, where the feasibility of constraints CT13-CT21 can be ensured using MIP of Constraint Programming models.

The low-level subproblems compute independently for all scenario the production plans optimizing the production costs and fulfilling the constraints CT1-CT6 and CT12 having fuel levels and maintenance dates fixed. Such production problems can be solved using greedy strategies following increasing production costs or using LP problems.

We note that the operational approach in the French Utility Company was in this scope before the Challenge, we refer to [19]. The approaches in competition were not among the most efficient, which can be analyzed comparing with the solving characteristics of frontal local search approaches.

3.4 Local search approaches

The last category gathers local search approaches, tackling the problem straightforward without decomposition. In this category, we find the two best approaches of the Challenge in terms of solution quality, having both very similar results, a methodology similar with LocalSolver [13, 3] and an unpublished Simulated Annealing local search [24].

The crucial points are in this case the choices of the neighborhoods to have efficient moves with the possibility to calculate quickly the cost of a locally modified solution, to be able to explore aggressively a lot of solutions. The final neighborhoods chosen by Gardi et al (2011) are reported in [13]:

  • k-MoveOutagesRandom: select k outages among T2 plants randomly and move them randomly in “feasible” time intervals, that is, ensuring the respect of earliest and latest starting dates (CT13 constraints) and maximum stocks before and after refueling (CT11 constraints);

  • k-MoveOutagesConstrained: select T2 plants which are involved in combinatorial constraints related to outage scheduling (CT14-CT21) randomly, select k outages of these plants randomly, move these outages randomly in feasible time intervals;

  • k-MoveOutagesConsecutive: select a T2 plant randomly, select k consecutive outages of these plant randomly, move these outages randomly in feasible time intervals.

To explain the better quality of primal solutions of frontal local search compared to heuristic decompositions, heuristic decomposition are slowed down by the use of MIP,LP or CP computations. On the contrary, frontal local search explored aggressively more solutions. An other explanation is due to the neighborhoods chosen for decomposition approaches, that were mostly very small (e.g allowing to move only one an outage from only one week). Such neighborhoods were used by [24] in the qualification phase, and lead to bad solutions. The effort to diversify and get larger neighborhoods is mentioned in [13] to have a neighborhood structure compatible to their aggressive exploration. We note that excellent results are reported in [13] with a time limit of few minutes, which is interesting for an industrial application.

Exact approaches are also slowed down by the use of MIP, LP or CP computations, with also memory space problematic, as in [22] to tackle large size instances with few approximations. Tractable exact optimization based approaches cannot avoid over-costs due to their approximation model, which explains the superiority of approaches [13, 24] in this case.

4 MIP formulation

In this section, we provide a MIP formulation for the problem, relaxing only constraints CT6 and CT12 as in [22]. It leads to a MIP formulation where the only binary variables are the decisions of outage beginning weeks .

4.1 Definition of the variables

A major modeling difference with [22] is in the binary variable definitions of : we define if and only if the outage beginning week for unit ’s cycle is before week . Binary variables in [22] are equal to if and only if outage beginning week for cycle is exactly . This choice allows to have efficient branching following results of [9]. We extend the notations with for , for , for and .

Other continuous variables must be introduced: refueling quantities for each outage , T2 power productions at cycle , fuel stocks at the beginning of campaign (resp at the end) , T1 power productions , and fuel stock at the end of the optimizing horizon. We note that T2 power productions are duplicated for all cycle to have a linear model, if week is not included in the production cycle .

Beginning dates of outages of unit at cycle at .
Refueling levels of unit at cycle .
Production of T2 unit at cycle at .
T1 production of unit at week .
Refueling levels for unit at cycle .
Fuel stock of T2 unit at the end of the optimizing horizon.
Fuel levels of T2 unit at the beginning of cycle .
Fuel levels of T2 unit at the end of production cycle .
Table 3: Set of variables

4.2 MIP formulation relaxing CT6 and CT12

Relaxing CT6 and CT12, it gives rise to the MIP formulation below. (2) is required with definition of variables . (3) and (4) model CT13 time windows constraints: outage is operated between weeks and . (5) models CT1 demand constraints. (6) models CT2 bounds on T1 production. (7) models CT3, CT4 and CT5 bounds on T2 production. (8) models CT7 refueling bounds, with a null refueling when outage is not operated, ie . (9) writes CT8 initial fuel stock. (10) writes CT9 fuel consumption constraints on stock variables of cycles . (11) models CT10 fuel losses at refueling. (12) writes CT11 bounds on fuel stock levels only on variables which are the maximal stocks level over cycles . thanks to (10). (13) models CT11 min fuel stock before refueling, these constraints are active for a cycle only if the cycle is finished at the end of the optimizing horizon, ie if , which enforces to have disjonctive constraints where case implies a trivial constraints thanks to (12). (14) linearizes the constraints to enforce to be the fuel stock at the end of the time horizon. is indeed the such that and , for the disjonctive constraints (14) that write a trivial constraints in the other cases thanks to (12), we define . (15) is a common framework for scheduling constraints from CT14 to CT21, which was noticed independently in [22, 18].


4.3 Adding stability cost for dynamic reoptimization

Previously, the objective function (1) consider only the financial cost. With the extension of stability cost as defined in Section 2.2, we shall consider two objectives, the financial cost and the stability cost , defined with:


The two objectives can be adversarial, the solution of the bi-objective optimization is thus a set of Pareto non-dominated solutions. We note that in the case is defined by for all . and , i.e. counts the number of desorganization, this makes the objective integer with a granularity of , which is interesting considering an -constraint method with ([21]). In this case, the constraint to bound the number of modifications of the original planning by at most is given by:


Another application of this extension is to repair an initial solution of maintenance planning, that may be infeasible, looking for the closest feasible solution as in Feasibility Pump (FP, [4]) using or and minimizing or with small.

4.4 Adding lighter CT6 constraints

This section provides a lighter formulation for stretch constraints enforcing only upper bounds on the production as illustrated Figure 2. As illustrated in Figure 2, the equation linking modes et is:

where is the ratio maximal power in stretch / maximal power , and denotes the residual fuel stock.

To write mathematically constraints to have upper bounds for the T2 production, we use that stretch decreasing profile is concave. Thus, the production domains are defined with the intersection of the semi spaces defined with the equations of the different modes. Adding in the model variables denoting the residual fuel stock of the unit at time period , the intersection of the semi spaces gives rise to following stretch constraints


The definition of new variables requires following linking constraints to enforce to be the residual fuel:


where , such verifies . Indeed, if , week happens in cycle , the active constraint is , otherwise we have , which is trivial thanks to the definition of .

A special case is interesting: if ( ) do not depend on indexes , we can have a MIP formulation using less constraints than previously, writing the constraints for the global production power :


5 Constructive matheuristics

This section aims to design constructive matheuristics. The modeling and solving facilities with MIP ensures constraint feasibility, with heuristic reductions to have tractable MIP resolutions for large scale instances.

5.1 A very simplified MIP formulation

A first heuristic approximates the previous MIP formulation, noticing that nuclear energy is cheaper that flexible energy from T1 power plants. A MIP formulation is considered where the nuclear production is either null (during outages), either at its maximal level, ie . This can lead to infeasibilities, it requires some relaxations. Nuclear power can be upper than some low demands , leading infeasibilities with 5 written with an equality. Relaxing 5 with minimal production inequalities fixes such problems. Furthermore, if a maximal nuclear production can be seen as economically optimal, it can be infeasible with the fuel constraints. Indeed, a maximal nuclear production (and thus a maximal fuel consumption), can be in contradiction with CT13 earliest dates of outages and positivity of fuel levels Hence, we introduce continuous variables to authorize negative fuel stocks in soft constraints with high values of in the objective function, in a following simplified MIP formulation:

This simplified MIP formulation is used to build a first maintenance planning. Such planning requires only a LP computation to compute the optimal production cost related, ensuring feasibility of all the constraints. If this high level planning cannot lead to feasible production and fuel levels, a local reparation is required.

5.2 Restrict- Relax-and-Fix algorithm

To decrease the size of the MIP computation, Restrict- Relax-and-Fix (RRF, [16]) strategies partition variables and iterate with partial MIP computations with fixing and relaxation reductions. A RRF strategy can be implemented partitioning the variables following the cycle index . It computes successively the solutions for cycles , the dates of cycles are fixed by the previous optimizations, and the variables related to the cycles are relaxed continuously or even simplified thanks to the aggregation used in [11].

As is, some feasibility difficulties can occur when a cycle optimization fixes outages in a way that no solution exist thereafter to place the remaining outage because of scheduling constraints CT13-CT21 and/or fuel constraints CT11. A crucial point is also to have enough flexibility in the iterated fixing decisions to have feasible solutions for the next iterations. This induces in Algorithm 1 to consider binaries for cycle and for the iteration placing outages , and to relax continuously/simplify formulations only for cycles . The remaining fuel costs are relaxed also to focus on the local optimization around the cycles to optimize, aggregating the cycles in the cycle relaxing the outage and the fuel constraints. The MIP formulation considered for each iteration is thus:

Algorithm 1: Restrict- Relax-and-Fix algorithm
Initialisation: define matrix for all .
for to
solve the MIP with for
store the solutions
end for
solve the MIP with for
return the solution of the last MIP

5.3 Construct, Merge, Solve & Adapt algorithm

Another RRF decomposition iterating unit by unit could be considered. Such heuristic has several weaknesses a priori. Firstly, having an infeasibility at an iteration does not allow to finish the resolution with the sequentiality of the algorithm. Secondly, the choice of the iterated unit order can be crucial.

To face these two points, a Construct, Merge, Solve & Adapt (CMSA, [5]) is implemented, computing independently (for a possible parallel computation) planning for all units as if they were alone, seeing the other units as T1 power plants. Merging these planning has significant chances to be unfeasible for the scheduling constraints CT14-CT21, a second phase repairs the feasibility around this ”ideal” solution, which makes no asymmetry in the choice of the units.

Algorithm 2: Construct, Merge, Solve & Adapt
Input: a repairing operator
Initialisation: define matrix for all .
for each T2 unit :
define new demands
solve the MIP with and new demands
store the maintenance planning
end for
Merge the maintenance plannings .
Compute as an LP the related production and fuel decisions.
if the current planning is unfeasible:
repair the current solution with
end if
return the last current solution

6 POPMUSIC-VND matheuristic

Once a feasible solution is built with previous constructive matheuristics, this section improves the current solution in a local search procedure.

Algorithm 3: multi neighborhood descent with MIP neighborhoods
Input: an initial solutions, a set and order of neighborhoods to explore
Initialisation: currentSol = initSolution, initial neighborhood.
while the stopping criterion is not met
define the MIP with incumbent currentSol and the neighborhood )
define currentSol as warmstart
currentSol = solveMIP(MIP,timeLimit( ))
end while
return CurrentSolution

6.1 General algorithm

The general local search algorithm is described in the Algorithm 3. The neighborhoods are defined using MIP definitions, similarly to [1]. MIP neighborhoods allow to explore large neighborhoods using recent progresses of MIP solvers, to have fewer and better local optimums than small neighborhoods approaches. Many MIP neighborhoods can be defined, Algorithm 3 changes systematically the choice of the neighborhood within the local search, similarly with multi neighborhood search approaches. It induces that a local optimum for the whole local search is a local optimum for all the neighborhoods considered in Algorithm 3. This a second property which ensures to have less and better local optimums than classical local search approaches. The stopping criterion could be a maximal time limit or a maximal number of iterations, or being in a local extremum for all neighborhoods.

Usually, MIP neighborhoods are defined for small sub-problems where the exact B&B computation to optimality is quick, as in [20]. In this work, MIP neighborhoods are defined with three characteristics for an efficient MIP neighborhood searc of primal solutions in small time limits. Parametrization is empiric, for a good trade-off between solution quality and time spent in MIP solving. The MIP neighborhoods are thus defined as following:

  • The restriction of search space: Variable fixations or other extra constraints that the current solution satisfies to have an easier B&B solving.

  • a B&B stopping criterion: it must defined so that the B&B search is efficient in a short solving time. Parametrizing short time limits is a first step, additional parameters like stopping a threshold in the gap between the best primal and dual bounds can be defined in order to not waste time in proving optimality without improving significantly the best primal solutions.

  • a specific MIP solving parametrization: for an efficient B&B search in the defined time limit. Easy neighborhoods need no specific MIP parametrization, for a better efficiency in short resolution time for difficult sub MIPs, emphasizing the heuristic search with mipemphasis parameter, disabling or limiting cutting plane passes with cutpass parameter (in the terminology of Cplex).

The current solution is the primal solution given by the last B&B resolution and it is also defined as warmstart for the next B&B resolution to improve the efficiency of B&B primal heuristics, enabling RINS or Local Branching heuristics from the beginning. This ensures that the solution given by the MIP resolution is at least as good as the current solution at each iteration. This algorithm is thus a steepest descent algorithm, similarly to MIP-VND as in [20].

6.2 MIP neighborhoods

Multiple types of large and variable neighborhoods can be defined. We denote with the beginning week of outage in the current incumbent solution. Neighborhoods define strategies to fix variable around the solution defined by . Generic neighborhoods are first derived from [8, 12]:

  • : Similarly to RINS heuristic [8], variables are fixed if they are a common integer value in the LP relaxation and in the current solution.

  • : for k an integer, k-Local Branching neighborhood allows only modifications to the incumbent, adding constraints as in [12]:

We note that neighborhoods have the same number of variables than the original problem, the acceleration opf the MIp solving is induced by the less nuimber of branching to explore the space thanks to the additional constraints.

The multi-index structure allows to define neighborhoods fixing variables along a type of indexation:

  • unit selection: only T2 units are reoptimized.

  • : all outages are reoptimized in the time window .

  • : variables relative cycles with are reoptimized.

would be a reoptimization with a fixed radius around the incumbent, is thus a natural neighborhood widely used for the Challenge. The linear dependence with induces a ”funnel” structure, to model that a move on the first outages can imply larger moves for the succeeding outages.

The weakness of such neighborhoods is that it may not be possible to swap the order of outages constrained with scheduling constraints like CT14 with short neighborhoods. is a first answer to have complementary neighborhoods to swap the order of outages. can be implemented for all single units, in this case, it authorizes large neighborhoods, and the operational “winter jump” neighborhoods: it may be interesting to change radically an outage scheduling, from an Autumn period to a Spring period, to avoid Winter periods where the substituted production cost is higher.

This illustrated that VND-MIP allows to define easily neighborhoods induced by operational expertise. In that case, the genericity of modeling and resolution is useful to avoid specific implementations imposed by traditional local search. An application of MIP-VND scheme is also to design and test neighborhoods with few implementation efforts, to select neighborhoods that should be investigated for a proper implementation.

6.3 Sequences of neighborhoods

A key point in the Algorithm 3 is the sequence of neighborhoods. Applying Algorithm 3 with a single neighborhoods, it allows to analyze the impact of neighborhoods in the quality and number of local extrema. In this case, different neighborhoods are compared starting with the same initial solutions, and analyzing the different qualities of the local minimum where the steepest descent local search converges. Here, there is a specific interest to an analyze the quality of neighborhoods , that contains widely used neighborhoods with decomposition approaches for the Challenge. It is interesting also to analyze the impact to add swap possibilities with or “winter spring” with .

A traditional idea with VND is to increase the size of the neighborhoods when a local minimum is reached, for a better compromise between the computation time and the local minimum reached. It can be implemented for Local Branching neighborhoods with increasing , as developed in [23]. It is similar with neighborhoods . In our case study, long computation time are allowed in the operational process for the monthly reoptimization of the maintenance planning (for instance several hours of computation by night). We do not seek to find the best solutions in short resolution time, where the previous sequencing strategy of neighborhoods is relevant. Having many types of neighborhoods, the stopping criterion is firstly to reach a local minimum for all the types of neighborhoods. The choice of neighborhoods is nested, alternating deterministically the order of neighborhoods and stopping the resolution when no neighborhood can improve the current best solution.

Using neighborhoods or , we can design different partitions of the binary variables. A partition of units give rise to neighborhoods for partitioning the variables of the problem. Many such partition can be considered, the partition units by units with is a first natural choice, reoptimizing all the planning separately. To control the size of the MIP problems, one can consider partitions with subsets of at least elements where is given. In the applications, natural partition occurs in sites with several nuclear reactors, in subsets of the constraints CT14-CT15. VND iterations can also try successively to reoptimize along different partitions, similarly with the partitioning local search POPMUSIC (Partial optimization meta-heuristic under special intensification conditions, see [28]). neighborhoods for all define also a partition, ”orthogonal” to the unit reoptimization. Considering all the does not define a partition of the variables anymore, but this is not a fundamental property for the VND.

7 Computational results

The computational analyses have two goals. On one hand, we are interested in methodological results, examining the limits of MIP solving and matheuristics and trying to explain some results/ranking after the ROADEF Challenge. On the other hand, industrial results analyze the impact of different modeling, e.g. the impact of CT6 constraints and of the extension with stability costs, for a better operational application. Tests were computed with a laptop running Linux Ubuntu 14.04 with an Intel Core2 Duo processor, 2.80GHz. Our implementation used the modeling language OPL to solve MIP with Cplex 12.5 and OPL script to iterate MIP computations. If this implementation is far from being optimal, this step is useful to understand better an industrial problem, to design efficient meta-heuristic operators with a low implementation effort, for a further proper implementation of selected operators.

Instances I J K S T W varBin Troncat. varBin
A1 10 11 6 10 1750 250 3892 A1_3_120 245
A2 18 21 6 20 1750 250 7889 A2_3_120 663
A3 18 21 6 20 1750 250 8162 A3_3_120 568
A4 30 31 6 30 1750 250 17465 A4_3_120 1305
A5 28 31 6 30 1750 250 15357 A5_3_120 1868
B6 50 25 6 50 5817 277 24563 B6_3_120 1519
B7 48 27 6 50 5565 265 35768 B7_3_120 4658
B8 56 19 6 121 5817 277 69653 B8_3_120 11057
B9 56 19 6 121 5817 277 69306 B9_3_120 11146
B10 56 19 6 121 5565 265 29948 B10_3_120 1816
X11 50 25 6 50 5817 277 20081 X11_3_120 1470
X12 48 27 6 50 5523 263 27111 X12_3_120 1927
X13 56 19 6 121 5817 277 30154 X13_3_120 2838
X14 56 19 6 121 5817 277 30691 X14_3_120 2844
X15 56 19 6 121 5523 263 27233 X15_3_120 1784
Table 4: Characteristics of the instances ROADEF and their truncated ones to 3 cycles and 120 weeks: I,J number of T2 and T1 units, K number of cycles for all T2 unit, S number of scenarios, T,W, number production and weekly time steps, varBin denotes the number of binaries

7.1 Instances characteristics

We used the datasets provided by EDF for the Challenge, after a preprocessing to aggregate production time steps to week and the scenarios using weighted sums with respective weights the length of the time stem and the probability of occurrence of a scenario. The full data characteristics are provided in Table

4. Dataset A contains five small instances in a horizon of five years, with 10 to 30 nuclear units having production cycles. Instances B and X are more representative of real-world size instances, with 20 to 30 T1 units, around 50 nuclear units. These datasets are now public, dataset X was secret for the challenge, with instance characteristics that had to be similar to dataset B. This symmetry holds for the number of units, cardinal of sets of Table 1 and number of constraints of each type, but this does not hold for the number of binaries; Table 4 shows in the column nbVarBin the number of binary variables, the cumulated amplitude of time windows. Instances B8 and B9 are more combinatorial, with no time window constraints for cycles , which are likely the most representative instances from the real-world application. A special interest is given to be able to solve efficiently instances B8 and B9, B7 is also difficult considering the nbVarBin indicator. Lastly, we used truncated instances planning at least maintenances in the first weeks of the previous instances to have easier instances and possible comparisons of heuristics to optimally proven solutions. Such instances are suffixed with _3_120. Table 4 shows also the number of binaries of these smaller instances.

Instances no PP + PP +warmstart penal
A1_3_120 0,12 0,12 0,04 0,13
A2_3_120 0,31 0,29 0,27 0,58
A3_3_120 0,23 0,21 0,21 0,14
A4_3_120 1,55 0,89 0,45 0,39
A5_3_120 127 82 70,86 15,3
Total A 129,21 83,51 71,83 16,54
B6_3_120 32,3 17 13,8 4,98
B7_3_120 916,9 512 437,9 20,1
B8_3_120 3600 1483 1461 45
B9_3_120 3600 3600 3600 264
B10_3_120 60,8 52 43,8 5,4
Total B 8210 5664 5556,5 339,48
X11_3_120 60,5 40,7 21,75 5,7
X12_3_120 111 69,4 43,6 17,2
X13_3_120 27,6 22 16,2 5,98
X14_3_120 98,8 56,1 53,2 16,3
X15_3_120 79,3 14 8,7 15,4
Total X 377,2 202,2 143,45 60,58
Total 8716,41 5949,71 5771,78 416,6
Table 5: Comparison of termination time of B&B to optimality: without and with the exact pre-processing of [10], preprocessing and B&B warmstart with the optimal (or best known) solution, and with penalization costs to the baseline solution given by 5.1

7.2 MIP solving characteristics

Table 5 compares the termination time of B&B to optimality on the smallest instances: without and with the exact pre-processing of [10], the influence of warmstarting the B&B algorithm with the optimal (or best known) solution. On short instances, the B&B algorithm converges in less than one hour proving the optimality for all instances except for B8_3_120 and B9_3_120. The exact pre-processing of [10], tightening time windows thanks to constraints CT5 and CT8-11, has an strong influence in the computation times to optimality. It allows toprove the optimality in one hour for B8_3_120. Warmstarting has less influence, it is however relevant for the application in Algorithm 3.

Table 6 indicates the quality of the solutions found in one hour computation with B&B for the largest instances comparing the gaps th the Best Known Solutions (BKS) and provides the dual bounds in terms of gaps to the BKS using the results of [10] with 1h computations. These results on the deterministic average scenario justify the commonly used simplification to aggregate production time steps to weekly time steps, and to aggregate scenarios. Using Cplex 12.3, the frontal resolution is completely inefficient on instances B8 and B9, 1h resolution time is not enough to compute the LP relaxation. The sizes of these instances are a strong limiting factor. For the other instances, the frontal MIP resolution is efficient, with low gaps between the best primal and dual bounds in a 1h resolution time. Using Cplex 12.5, LP relaxations can be computed for all instances in 1h.

However, the solutions of the LP relaxation are not useful to compute primal solutions for B8 and B9. In these cases, several production cycles can overlap and the LP relaxation gives very few integer variables, specially for last cycles. Hence, this is a big handicap for MIP primal heuristics, relying on the LP relaxation such as Feasibility Pump, RINS or Local Branching, and tailored variable fixing heuristics. The MIP frontal resolution is always inefficient on B8 and B9, if LP relaxation can be computed in one hour, MIP primal heuristics are inefficient. These resolution facts explains and justify the approximations made for the Challenge by [22, 18, 27].

Figure 3: Instance B7-3-120
Figure 4: Instance X11-3-120
Figure 5: Pareto frontiers for the best compromises cost/stability.

7.3 Optimizing around an initial solution

In the instances from the ROADEF Challenge, there were no initial solution given. To analyze the impact of having an initial solution and minimizing with stability costs, we considered the initial solutions given by the matheuristic of Section 5.1 without the repairing post-processing. Such initial solutions are not always feasible, and are of medium quality, which seemed a good choice to represent monthly reoptimization, the current solution computed the previous month can become infeasible or improvable after the realization of real-life uncertainties.

Adding penalization to a baseline solution in the objective function changes fundamentally the MIP convergence characteristics, guiding the solution search around the baseline solution improves very significantly the resolution time, as shown in Table 5. The higher are the stability coefficients , the more the B&B algorithm is accelerated. An explication is that a lot of solutions with similar ROADEF costs exist which was already noticed for the Challenge. It induces “pseudo-symetries” in the B&B tree search, which is known to be a bottleneck for the B&B search. It explains the difficulty of the B&B convergence previously outlined.

An application of these properties is to repair partial infeasible solutions built by a constructive heuristic, with a high penalization of the distance to the partial solution, similarly to Feasibility Pump ([4]).

These good B&B characteristics allow easily to considering an -constraint method to draw the Pareto curve of non dominated solutions in the bi-objective optimization cost/stability ([21]). We note that the granularity for the -constraint method can be set to along integer objectives measuring the stability penalization, which induces the whole Pareto frontier without granularity errors. Figure 5 presents such Pareto frontiers for instances B7-3-120 and X11-3-120 with the baseline solution given by the constructive heuristic of section 5.1. The quality of the baseline solution has a consequent impact for the Pareto frontiers, with a poor quality solution like the ones given by the constructive heuristic of section 5.1, an important trade-off exist between cost and stability, whereas good initial solutions produce a flat Pareto frontier and in the case where the baseline solution is optimal the Pareto frontier would be reduced to a single point.

BKS Dual Frontal Simp CMSA R-R&Fix VND
B6 76966 0,13% 0,00% 3,22% 3,8% 0,34% 0,00%
B7 74234 0,52% 0,02% 3,68% 0,75% 0,25% 0,00%
B8 73240 7,95% NS 52,54% 52,21% 2,02% 0,00%
B9 72812 6,69% NS 27,46% 28,95% 1,15% 0,00%
B10 69501 0,07% 0,02% 0,21% 0,2% 0,21% 0,00%
X11 73018 0,37% 0,00% 0,45% 0,53% 0,38 0,00%
X12 70604 0,20% 0,00% 0,25% 0,24% 0,31% 0,00%
X13 69231 1,09% 0,04% 7,37% 9% 0,71% 0,00%
X14 68395 0,89% 0,02% 6,02% 1,21% 0,54% 0,00%
X15 66029 0,10% 0,06% 0,12% 0,12% 0,38% 0,00%
Total 714031 0,00% NS 10,31% 9,89% 0,63% 0,00%
Table 6: Result comparison for constructive primal math-heuristics
Instances Init
B6 3,22% 0,54% 0,13% 0,26% 0,13% 0,13%
B7 3,68% 1,14% 0,41% 0,47% 0,35% 0,19%
B8 52,54% 38,34% 19,70% 31,12% 16,05% 26,73%
B9 27,46% 15,91% 7,54% 12,35% 7,54% 8,89%
B10 0,21% 0,1% 0,05% 0,06% 0,02% 0,05%
Total B 17,40% 11,18% 5,55% 8,83% 4,80% 7,18%
X11 0,45% 0,31% 0,10% 0,26% 0,22% 0,27%
X12 0,25% 0,13% 0,10% 0,09% 0,07% 0,06%
X13 7,37% 2,23% 1,30% 1,47% 0,93% 1,01%
X14 6,02% 1,64% 0,48% 0,61% 0,45% 0,54%
X15 0,12% 0,08% 0,06% 0,07% 0,06% 0,06%
Total X 2,82% 0,87% 0,41% 0,50% 0,34% 0,39%

B6 3,22% 0,26% 0,03% 0,04% 0,31% 0,13%
B7 3,68% 0,34% 0,06% 0,02% 0,55% 0,40%
B8 52,54% 1,18% 0,36% 0,36% 35,79% 27,63%
B9 27,46% 1,31% 0,05% 0,05% 14,48% 12,12%
B10 0,21% 0,05% 0,03% 0,03% 0,08% 0,03%
Total B 17,40% 0,63% 0,11% 0,10% 10,21% 8,04%
X11 0,45% 0,16% 0,09% 0,09% 0,30% 0,26%
X12 0,25% 0,08% 0,05% 0,05% 0,15% 0,10%
X13 7,37% 0,51% 0,07% 0,07% 2,43% 0,28%
X14 6,02% 0,28% 0,13% 0,13% 1,40% 0,60%
X15 0,12% 0,09% 0,05% 0,05% 0,10% 0,08%
Total X 2,82% 0,22% 0,08% 0,08% 0,87% 0,26%
Table 7: Quality of local minimums considering one single type of neighborhoods in the VND of Algorithm 3

7.4 Constructive matheuristics

With the characteristics of B&B search previously outlined, the matheuristic effort to design constructive approaches focused on the ability to get feasible solutions for the most difficult instances B7, B8 and B9. The numerical results of constructive heuristics are highlighted in Table 6. Simplified MIP computations of section 5.1, the CMSA and the R-R&F heuristics succeeded to provide feasible solutions for all types of instances.

Focusing the binaries of the problem, the MIP formulation of section 5.1 accelerated significantly the search of primal solutions, easing the MIP primal heuristics. For the instances B8 and B9, feasible solutions were found in 15 minutes including the repairing post-processing. R-R&F approach gives very good solutions, but require sometimes long calculus times to converge to a feasible solution in the MIP iterations (around 20 or 30 minutes to compute the most difficult computations in intermediary cycles, for a total solving time in more than one hour). The CMSA strategy using the previous repairing strategy to build a feasible solution from partial solutions built by single unit optimization. It allows to build quickly primal solutions for all instances, but with a lower quality of primal solutions than the relax-and-fix strategies.

7.5 POPMUSIC-VND matheuristics

Table 7 compares the quality of the local minimum considering each type of single neighborhoods separately int the MIP-VND local search.

Neighborhoods include the natural neighborhood to move one single outage from one week. These neighborhoods were considered by many approaches for the ROADEF Challenge. includes also 1-translation neighborho