Matheuristics to optimize refueling and maintenance planning of nuclear power plants

TL;DR

This paper develops and tests matheuristics combining MILP techniques and local search to optimize maintenance and refueling schedules for nuclear power plants over five-year periods, improving solution quality for large instances.

Contribution

It introduces effective matheuristics for large-scale nuclear maintenance planning, demonstrating their ability to handle complex constraints and reduce over-costs associated with time window restrictions.

Findings

Matheuristics are effective for medium and large instances.

Time window restrictions cause significant over-costs.

Matheuristics enable extension to bi-objective optimization.

Abstract

Planning the maintenance of nuclear power plants is a complex optimization problem, involving a joint optimization of maintenance dates, fuel constraints and power production decisions. This paper investigates Mixed Integer Linear Programming (MILP) matheuristics for this problem, to tackle large size instances used in operations with a time scope of five years, and few restrictions with time window constraints for the latest maintenance operations. Several constructive matheuristics and a Variable Neighborhood Descent local search are designed. The matheuristics are shown to be accurately effective for medium and large size instances. The matheuristics give also results on the design of MILP formulations and neighborhoods for the problem. Contributions for the operational applications are also discussed. It is shown that the restriction of time windows, which was used to ease computations, induces large over-costs and that this restriction is not required anymore with the capabilities of matheuristics or local search to solve such size of instances. Our matheuristics can be extended to a bi-objective optimization extension with stability costs, for the monthly re-optimization of the maintenance planning in the real-life application.

Figures

Tables

Table 1. Table 1: Definitions and notations for the input parameters

$w\in 𝒲=[[1,W]]$	Weekly time steps.
$j\in 𝒥=[[1,J]]$	Flexible (Type 1, T1) power plants.
$i\in ℐ=[[1,I]]$	Nuclear power plants (Type 2, T2).
$k\in 𝒦=[[0,K]]$	Cycles related to T2 units, $k=0$ for initial conditions.
$c\in 𝒞$	Scheduling constraints, type CT14 to CT21 in the specification [38].

Table 2. Table 2: Description of the constraints considered in this paper, with the reference to the corresponding constraints from the EURO/ROADEF Challenge

CT1, demand covering: for all time step $w\in 𝒲$, the total production of T1 and T2 power plants must equalize the demands ${𝐃}_w$.
CT2, T1 production bounds for all $w\in 𝒲$, the production domain of T1 plant $j\in 𝒥$ describes exactly the continuous domain $[{\underset{¯}{𝐏}}_{j,w},{\overline{𝐏}}_{j,w}]$
CT3-5, T2 production bounds The productions of offline T2 power plants are null during outages. Otherwise, the production domain of T1 plant $j\in 𝒥$ describes exactly the continuous domain $[0,{\overline{𝐏}}_{i,w}]$ for all period $w\in 𝒲$.
CT7, refueling quantities: The refueling possibilities for outage $k$ of T2 plant $i$ describes exactly the continuous domain $[{\underset{¯}{𝐑}}_{i,k},{\overline{𝐑}}_{i,k}]$
CT8, initial fuel stock: The initial fuel stock for T2 unit $i$ is ${\mathrm{𝐗𝐢}}_{𝐢}$.
CT9, fuel consumption: The fuel consumption for T2 unit $i$ at period $w$ is proportional to the generated power, with a proportional factor $-{𝐅}_w$.
CT10, losses in refueling After an outage, the new fuel stock is the sum of the refueling quantities (cf CT7) and the amount of unspent fuel after the last production cycle. A proportional loss factor is applied to the residual fuel before refueling.
CT11, bounds on fuel stock The fuel level is bounded in $[0,{𝐒}_{i,k}]$ for cycle $k$ of T2 unit $i$. The fuel level must be lower than ${𝐀}_{i,k+1}$ to process outage $k+1$.
CT13, time windows The maintenance cycles follow the set $k\in 𝒦$ without skipping cycles: if outage $k+1$ is processed, it must follow the cycle $k$. The beginning dates of outage $k$ of T2 unit $i$ must be included in $[{\mathrm{𝐓𝐨}}_{i,k},{\mathrm{𝐓𝐚}}_{i,k}]$.
CT14-21, resource scheduling constraints : For each constraint $c\in 𝒞$, a subset of outages $A_c$ and a subset of time periods $W_c$ is considered. For all outage $a\in A_c$ and for all time period $w\in W_c$, a resource consumption $R_{c,a,w}$ is defined for the outage $a$ being scheduled at $w$, and the global resource consumption at $w$ must be lower than ${\overline{R}}_{c,w}$

Table 3. Table 3: Characteristics of the ROADEF instances : 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

Instances	I	J	K	S	T	W	varBin
A1	10	11	6	10	1750	250	3892
A2	18	21	6	20	1750	250	7889
A3	18	21	6	20	1750	250	8162
A4	30	31	6	30	1750	250	17465
A5	28	31	6	30	1750	250	15357
B6	50	25	6	50	5817	277	24563
B7	48	27	6	50	5565	265	35768
B8	56	19	6	121	5817	277	69653
B9	56	19	6	121	5817	277	69306
B10	56	19	6	121	5565	265	29948
X11	50	25	6	50	5817	277	20081
X12	48	27	6	50	5523	263	27111
X13	56	19	6	121	5817	277	30154
X14	56	19	6	121	5817	277	30691
X15	56	19	6	121	5523	263	27233

Table 4. Table 4: Set of variables

$d_{i,k,w}\in \{0,1\}$	Beginning dates of outages of unit $i$ at cycle $k$ at $w$.
$r_{i,k}⩾0$	Refueling levels of unit $i$ at cycle $k$ .
$p_{i,k,w}⩾0$	Production of T2 unit $i$ at cycle $k$ at $w$.
$p_{j,w}⩾0$	T1 production of unit $j$ at week $w$ .
$r_{i,k}⩾0$	Refueling levels for unit $i$ at cycle $k$.
$x_i^f⩾0$	Fuel stock of T2 unit $i$ at the end of the optimizing horizon.
$x_{i,k}^i⩾0$	Fuel levels of T2 unit $i$ at the beginning of cycle $k$.
$x_{i,k}^f⩾0$	Fuel levels of T2 unit $i$ at the end of production cycle $k$.

Table 5. Table 5: Comparison of the number of binaries for the original instances of the ROADEF Challenge, and the reduced and extended instances. The numbers of binaries are shown without pre-processing and using the exact pre-processing from [ 16 ]

	ROADEF	ROADEF	_3_120	_3_120	Ext3	Ext3	Ext2	Ext2
Instance		+PP		+PP		+PP		+PP
A1	3878	497	171	171	-	-	-	-
A2	7853	1139	522	345	-	-	-	-
A3	8124	1017	398	298	-	-	-	-
A4	17384	2115	956	554	-	-	-	-
A5	15293	3008	1314	1163	-	-	-	-
B6	24485	3697	801	801	43069	10747	56179	20488
B7	35586	7500	3373	1232	43390	12736	52639	21519
B8	69378	26014	10825	3178	70876	28683	72367	30362
B9	69026	29208	10808	3663	70549	31915	72042	33942
B10	29843	4533	934	934	46332	11571	60302	23010
X11	20015	3610	810	810	43015	11657	56149	22390
X12	27000	4837	1249	1014	40362	11237	51672	20837
X13	30019	5169	1994	1118	49822	13236	63406	23241
X14	30559	5824	1975	1165	49831	14630	63404	26079
X15	27140	4144	983	983	45966	12272	59840	23785

Table 6. Table 6: Cost values for the BKS of the medium and large instances, the over-cost gap to the BKS of the ROADEF instances are computed for the Ext2 and Ext3 extensions to show the impact of a heuristic reduction of time windows.

	ROADEF	Ext3	Ext3	Ext2	Ext2
	BKS	BKS	gap	BKS	gap
B6	76966,1M	76734,8M	0,30 %	76407,1M	0,73 %
B7	74232,0M	73739,0M	0,66 %	73544,9M	0,93 %
B10	69501,5M	69143,1M	0,52 %	68834,6M	0,96 %
X11	73017,8M	72023,1M	1,36 %	71801,2M	1,67 %
X12	70599,5M	70240,4M	0,51 %	70255,8M	0,49 %
X13	69227,1M	68352,4M	1,26 %	68063,2M	1,68 %
X14	68395,3M	67799,8M	0,87 %	67496,5M	1,31 %
X15	66028,5M	65638,8M	0,59 %	65271,8M	1,15 %
Total	567967,8M	563671,4M	0,76 %	561675,1M	1,11 %

Table 7. Table 8: Comparison to the gaps to the BKS for different LB and UB at different stage of the B&B solving: LB of the LP relaxation, LB and UB at the root note after cutting plane passes and before branching, and after 1h and 3h computations.

	LP	LP	nod1	nod1	nod1	1h	1h	3h	3h
	t	LB	t	LB	UB	LB	UB	LB	UB
B6	4,4	2,59 %	25,6	0,92 %	0,45 %	0,20 %	0,00 %	0,06 %	0,00 %
B7	17,8	4,74 %	220,4	3,37 %	-	2,49 %	0,03 %	2,02 %	0,02 %
B8	147	12,23 %	1,1k	11,82 %	-	10,81 %	-	10,81 %	-
B9	167	11,62 %	1,5k	10,92 %	-	10,39 %	-	10,39 %	-
B10	8,3	3,97 %	36	1,76 %	0,12 %	0,44 %	0,00 %	0,23 %	0,00 %
X11	3,5	1,82 %	30	0,73 %	0,85 %	0,17 %	0,01 %	0,10 %	0,00 %
X12	7,7	3,25 %	45	1,20 %	0,39 %	0,37 %	0,01 %	0,12 %	0,00 %
X13	5,5	3,11 %	131	1,96 %	-	1,38 %	0,16 %	1,14 %	0,01 %
X14	6,3	3,56 %	64	2,08 %	1,75 %	1,34 %	0,07 %	1,12 %	0,05 %
X15	6,3	3,95 %	25	0,86 %	0,15 %	0,09 %	0,00 %	0,05 %	0,00 %

Table 8. Table 9: Comparison to the gaps to the BKS for the solutions of matheuristics, to Cplex solving and the dual bounds from [ 18 ] .

	$\overline{t}$	$\overline{g}$	$\sigma g$	Q1	Q3	$N_{0.01}$	$N_{0.05}$	$N_1$
Cplex 1h	217	0,00 %	0,00 %	0,00 %	0,00 %	18	18	18
RRF	238	0,07 %	0,23 %	0,01 % %	0,01 %	11	16	18
RRF+1VND	270	0,01 %	0,05 %	0,00 %	0,00 %	17	17	18
RRF+VND	313	0,01 %	0,05 %	0,00 %	0,00 %	17	17	18
Simpl	59	0,08 %	0,24 %	0,00 %	0,02 %	10	15	18
Simpl+1VND	93	0,01 %	0,05 %	0,00 %	0,00 %	17	17	18
Simpl+VND	137	0,01 %	0,05 %	0,00 %	0,00 %	17	17	18
CMSA	63	0,03 %	0,07 %	0,01 %	0,03 %	6	16	18
CMSA+1VND	96	0,01 %	0,05 %	0,00 %	0,00 %	16	17	18
CMSA+VND	133	0,01 %	0,05 %	0,00 %	0,00 %	16	17	18
Best VND		0,01 %	0,05 %	0,00 %	0,00 %	17	17	18

Table 9. Table 10: Comparison of termination time of B&B to optimality: without and with the exact pre-processing of [ 16 ] , pre-processing and B&B warmstart with the optimal (or best known) solution, and with penalisation costs to the baseline solution given by 5.1. These results were obtained using Cplex 12.5.

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
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

Table 10. Table 11: Comparison of the gaps to the BKS for the POPMUSIC-VND and for the VND with single types of large neighborhoods

Instances	Init	${𝒩}_{(1,3)}^{TW}$	${𝒩}_{k,k+1}^{cycles}$	${𝒩}_I^{units}$	POPMUSIC-VND
				$\left|ℐ\right|=8$
B6	3,22%	0,13%	0,13%	0,04%	0,00%
B7	3,68%	0,35%	0,4%	0,02%	0,00%
B8	52,54%	16,05%	27,63%	0,36%	0,00%
B9	27,46%	7,54%	12,12%	0,05%	0,00%
B10	0,21%	0,02%	0,03%	0,03%	0,00%
X11	0,45%	0,22%	0,26%	0,09%	0,00%
X12	0,25%	0,07%	0,1%	0,05%	0,00%
X13	7,37%	0,93%	0,28%	0,07%	0,00%
X14	6,02%	0,45%	0,6%	0,13%	0,00%
X15	0,12%	0,06%	0,08%	0,05%	0,00%
Total	10,31%	2,63%	4,26%	0,09%	0,00%

Table 11. Table 12: Result comparison for constructive primal math-heuristics

	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 12. Table 13: Quality of local minimums considering one single type of neighborhoods in the VND of Algorithm 3

Instances	Init	${𝒩}_{(0,1)}^{TW}$	${𝒩}_{(1,2)}^{TW}$	${𝒩}_{(0,3)}^{TW}$	${𝒩}_{(1,3)}^{TW}$	${𝒩}_{(0,5)}^{TW}$
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%
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	10,31%	6,17%	3,05%	4,78%	2,63%	3,87%

Table 13. Table 14: Over-Costs of the outage solutions with the MILP of section 4.2 in the model adding CT6 constraints, and local improvements with two iterations of VND with neighborhoods 𝒩 ( 0 , 2 ) T ​ W superscript subscript 𝒩 0 2 𝑇 𝑊 {\mathcal{N}}_{(0,2)}^{TW} .

	BKS	Projected 4.2	+LS
A1_3_120	7612M	0,01%	0,00%
A2_3_120	7213M	0,01%	0,00%
A3_3_120	7638M	0,02%	0,00%
A4_3_120	5 122M	0,05%	0,00%
A5_3_120	5 852M	0,08%	0,03%
B6_3_120	31 950M	0,10%	0,00%
B7_3_120	30 364M	0,18%	0,06%
B8_3_120	29 438M	0,14%	0,13%
B9_3_120	29 239M	0,18%	0,00%
B10_3_120	27 957M	0,34%	0,06%
X11_3_120	30 199M	0,20%	0,00%
X12_3_120	29 669M	0,25%	0,00%
X13_3_120	29 148M	0,36%	0,00%
X14_3_120	28 438M	0,23%	0,00%
X15_3_120	27 722M	0,37%	0,00%

Equations

$$v_0=\underset{d\in \{0,1{\}}^N,r,p,x⩾0}{\min }$$

$$\forall i,k,w,$$

$$\forall i,k,$$

$$\forall w,$$

$$\forall j,w,$$

$$\forall i,k,w,$$

$$\forall i,k,$$

$$\forall i,$$

$$\forall i,k,$$

$$\forall c\in \mathcal{C},w\in W_c,$$

$${\mathbf{C}}^{roadef}=\sum _{i,k}{\mathbf{C}}_{i,k}^r{r}_{i,k}+\sum _{j,w}{\mathbf{C}}_{j,w}^p{\mathbf{F}}_w{p}_{j,w}-\sum _i{\mathbf{C}}_i^f{x}_i^f$$

$${\mathbf{C}}^{desorg}=\sum _{i,k,w}{\mathbf{C}}_{i,k,w}^{pen}(d_{i,k,w}-d_{i,k,w-1})$$

$$\sum _{i,k}(1+d_{i,k,W^0(i,k)-1}-d_{i,k,W^0(i,k)})⩽\mathbf{Nmax}$$

$$v^{simpl}=\min$$

$$\forall i,k,$$

$$\forall j,w,$$

$$\forall i,w,$$

$$\forall w,$$

$$v_{RRF}(k^0)=\min$$

$$\forall i,k⩽k_0+1,w,$$

$$\forall i,k⩽k_0+1,$$

$$\forall w,$$

$$\forall j,w,$$

$$\forall i,k⩽k_0+1,w,$$

$$\forall i,k⩽k_0+1,$$

$$\forall i,k⩽k_0,$$

$$\forall i,k⩽k_0+1,$$

$$\forall c\in \mathcal{C},w\in W_c,$$

$$\sum _{i,k,w:w⩾W_{i,k}^o}{d}_{i,k,w}+\sum _{i,k,w:w<W_{i,k}^o}(1-d_{i,k,w})⩽k$$

$$\forall i,k,w,x_{i,w}⩽x_{i,k}^i-\sum _{w^{\prime }⩽w}{\mathbf{F}}_{w^{\prime }}{p}_{i,k,w^{\prime }}+M_i(1-d_{i,k,w}+d_{i,k-1,w})$$

$$\forall i,k,w,m>0,\frac{p_{i,k,w}}{{\overline{\mathbf{P}}}_i^t}⩽\frac{{\mathbf{c}}_{i,k,m-1}-{\mathbf{c}}_{i,k,m}}{{\mathbf{f}}_{i,k,m-1}-{\mathbf{f}}_{i,k,m}}(x_{i,w}-{\mathbf{f}}_{i,k,m})+{\mathbf{c}}_{i,k,m}$$

$$\forall i,w,m>0\sum _k\frac{p_{i,k,w}}{{\overline{\mathbf{P}}}_{i,w}}⩽\frac{{\mathbf{c}}_{i,m-1}-{\mathbf{c}}_{i,m}}{{\mathbf{f}}_{i,m-1}-{\mathbf{f}}_{i,m}}(x_{i,w}-{\mathbf{f}}_{i,m})+{\mathbf{c}}_{i,m}$$

Full text

This article should be cited as Dupin, N., Talbi, E. Matheuristics to optimize refueling and maintenance planning of nuclear power plants. J Heuristics (2020). https://doi.org/10.1007/s10732-020-09450-0

Keywords : Hybrid heuristics; Matheuristics ; Mixed Integer Programming ; Maintenance planning; Nuclear Power Plants ; Optimization in Energy

1 Introduction

Planning maintenance of nuclear power plants is a crucial issue in energy management, with huge financial stakes in France [40]. Indeed, 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 [30]. At a first glance, the nuclear power plants should be in maintenance when the demands in electricity are low, especially in summer. Technical and operational constraints must be satisfied for the maintenance planning, it makes difficult to plan the maintenance for human experts. Operations research approaches were investigated since the 70’s [11]. With the continuous progress of hardware and algorithms, especially in Mixed Integer Linear Programming (MILP), it became possible to tackle more complex and realistic problems in the decision support [21]. The sizes of the real-life instances has always been a challenging issue. Heuristics and decomposition techniques were used for the short term Unit Commitment (UC) problem, and also for the maintenance planning optimization of nuclear power plants [12, 30]. Recent progress in MILP allows an efficient resolution without decomposition for UC problems with additional real life constraints ([14]). Hybridizing heuristics using mathematical programming, known as matheuristics [34], allows to solve realistic size of instances [19, 17].

A problem of refueling and maintenance planning of nuclear power plants was specified by the French Utility company (EDF) in the 2010 ROADEF/EURO Challenge, a biennial competition of operations research [38]. The best results of the challenge were mainly obtained in [36, 22] with local search approaches. MILP formulations were also investigated, many constraints being linear in the specification [38]. Matheuristics were also used to tackle the large size instances of the problem, as investigated by [29, 33, 41]. The large size of the instances was a bottleneck for many approaches, especially on two instances, denoted B8 and B9, which can be seen as the most realistic instances for the real-life operations. In this paper, we focus on a deterministic version of the problem given for the Challenge. The difficulty of the instances B8 and B9 is explained and reproduced generating new instances from all the instances of the Challenge. The first motivation is to solve accurately the large instances which are relevant for real life applications. The second motivation is to give results related to the design of MILP formulations and heuristics for the problem.

This paper is organized as follows. In Section 2, we describe precisely the problem with its industrial context. In Section 3, we mention related works to appreciate our contributions. In Section 4, we present a MILP compact formulation. In Section 5, constructive matheuristics are designed using the previous MILP formulation. In Section 6, a local search matheuristics with MILP neighborhoods provides local improvements once feasible solutions are known. In Section 7, the computational results are discussed. In Section 8, our contributions are summarized, discussing also future directions of research. Appendix A discusses some industrial and numerical justifications related to the model definition, and Appendix B gathers tables of intermediate results.

2 Problem statement

This section describes the problem of refueling and maintenance planning, its constraints and objectives and some insights of the data characteristics and implementation constraints which induced our scientific problematic. The differences with the 2010 ROADEF/EURO Challenge are also highlighted.

2.1 Simplifications from the 2010 ROADEF/EURO Challenge

We consider in this paper the following simplifications of the constraint of the 2010 ROADEF/EURO Challenge [38]:

• •

weekly time steps: in the Challenge, there were two types of time steps, a weekly discretization for the outage dates and a hourly or daily discretization for power production and demands. In this study, we consider only weekly time steps, production and demand levels are aggregated weekly.

• •

determistic model: In the Challenge, discrete scenarios modeled the uncertainty of power demands, production capacities and costs. In this study, we consider only one deterministic scenario in the optimization.

• •

relaxation of CT6 constraint: CT6 constraints express deterministic production curves for nuclear power plants when the fuel is lower than a known threshold, as illustrated Figure 4. The relaxation of these constraints is justified numerically in Appendix A.

• •

relaxation of CT12 constraint: CT12 constraints express an aggregated constraint, limiting the modulation possibilities of nuclear power plants, to enforce the nuclear power plants to produce mostly at their maximal level. The relaxation of these constraints is justified in Appendix A, mentioning realistic constraints related to modulation of nuclear power plants.

2.2 Set and indexation

Table 1 gathers the notation for the sets and indexation. Two kinds of power plants are modeled. On one hand, Type-2 (or T2) power plants indexed with $i\in{\mathcal{I}}$, 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 $j\in{\mathcal{J}}$, model other power plants. Outages and production campaigns are indexed with the cycles $k\in{\mathcal{K}}$ for all T2 plants $i$. By convention, a cycle $k$ begins with the outage period for maintenance and refueling, before the production campaign $k$. The outage and production decisions are discretized weekly and indexed with $w\in{\mathcal{W}}=[\![1;W]\!]$. Lastly, $c\in{\mathcal{C}}$ gather scheduling constraints in outages. Such constraints were detailed in the sets from CT14 to CT21 in the Challenge specification [38]. Actually, one can model such constraints in one common framework, as mentioned in Table 2.

2.3 Constraint description

Table 2 defines the constraints with their nomenclature from CT1 to CT21 in the specification [38]. 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 process maintenance. 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 maintenance 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 Objective function

The objective function minimizes the power generation cost while satisfying the customer load for all time steps. Production costs of T1 units $j$ are proportional to the production levels. Production costs of T2 units are calculated proportionally to the fuel consumption: refueling costs are considered, and reduced with the costs of the remaining fuel with a proportional factor $\mathbf{C}_{i}^{f}$ to avoid end-of-side effects. This makes an aggregated cost function to optimize the global financial cost.

In operations, the maintenance planning is re-optimized monthly, the baseline planning impacts the re-optimized one. The maintenance operations are outsourced to several external companies fixing the periods of operations. Time window constraints CT13 can model such decisions with hard constraints. Actually, the maintenance planning can be slightly modified if the financial stakes worth the reorganization efforts. Hence, we consider an objective of stability, defining with $\mathbf{C}^{pen}_{i,k,w}$ the stability cost to reschedule the maintenance $(i,k)$ at week $w$. Denoting $W_{i,k}^{0}$ the initial beginning week of maintenance $(i,k)$, $\mathbf{C}^{pen}_{i,k,W_{i,k}^{0}}=0$. Any penalization with $\mathbf{C}^{pen}_{i,k,W_{i,k}^{0}}=0$ can be considered, like a linear penalization $\mathbf{C}^{pen}_{i,k,w}=|w-W_{i,k}^{0}|$, a quadratic one $\mathbf{C}^{pen}_{i,k,w}=(w-W_{i,k}^{0})^{2}$, or constant costs for modifications with $\mathbf{C}^{pen}_{i,k,w}=1$ for all $w\neq W_{i,k}^{0}$.

2.5 Data characteristics

The datasets provided for the Challenge are now public. The instance characteristics are provided in Table 3. Dataset A contains five small instances given for the qualification phase, in a horizon of five years, with 10 to 30 nuclear power plants having $6$ production cycles. Instances B and X are more representative of real-world size instances, with 20 to 30 T1 units and around 50 T2 units. These datasets are now public, dataset X was secret for the challenge, with instance characteristics that had to be symmetric to dataset B. This symmetry holds for the cardinal of the different sets , but this does not hold for the number of binaries (i.e. the total amplitude of time windows). Table 3 shows in the column varBin the number of binary variables. Instances B8 and B9 are more difficult. X13 and X14 are not symmetric with B8 and B9 in the number of binaries.

There is no time window constraints for cycles $k\geqslant 3$ for both instances B8 and B9, which is not the case for the others. Actually, B8 and B9 are the most relevant instances for the operational application. Maintenance operations are outsourced with several companies which makes difficult to reorganize maintenance dates in the first years of the planning horizon, whereas few such constraints apply for the last cycles of the planning horizon. A special interest is given in our study to solve efficiently and accurately instances like B7, B8 and B9.

2.6 Hardware and software limitations

For the Challenge, the hardware and software for the evaluation were specified. The computer was a bi-processor Intel Xeon 5420 2.5Ghz quad core (8GB of memory under Linux, 3GB under windows, 12MB of cache). Cplex could be used for MILP solving in the version 12.1.2.0 for the evaluation. Such hardware limitations are relevant for the industrial application, the software limitations may be updated to the current versions. The time limit was set to one hour per instance for the evaluation of the Challenge. This constraint was justified to be able to evaluate the approaches of many teams in the defined period (one month) between the deadline for submitting the final program and the announcement of the challenge results. For the real-life application, the planning is re-optimized monthly, several hours for computing time are reasonable, and possibly computations by night [13, 26].

3 Related work

This section describes the state of the art of the problem. The approaches of the 2010 EURO/ROADEF Challenge are of special interest, focusing on the related works both in the heuristic and mathematical programming approaches.

3.1 General facts of the Challenge ROADEF

The approaches in competition led to several publications of full papers ([2, 23, 24, 22, 29, 33, 41]). Slides (mainly presented at the EURO 2011 conference) are also available at http://www.roadef.org/challenge/2010/fr/methodes.php. We distinguish three main types of approaches: matheuristics based on MILP exact approaches, decomposition heuristics following the 2-stage structure and straightforward local search approaches. These three types of approaches are analyzed precisely in the next section.

A major difficulty came with the large size of instances. Many approaches used pre-processing reductions: tightening the time windows with exact or heuristic pre-processing; aggregating production time steps to week or aggregating the stochastic scenarios to the average one. The heuristic pre-processing requires to repair the solutions obtained in a post-processing considering the whole problem. Reparation strategies were in general effective to recover feasible solutions. The instances B8-B9 raised the most difficulties among the competing approaches , with many approaches leading to large over-costs for these instances whereas other instances have better results (we refer to the official results: http://www.roadef.org/challenge/2010/en/results.php).

3.2 Matheuristics based on MILP exact approaches

Two matheuristic approaches were designed for the Challenge. These matheurstics solve a simplified MILP problem before repairing the solution to ensure the feasibility for all the constraints and computing the cost in the full model. In both cases, the MILP model is solved heuristically using truncated exact methods, to solve the instances in one hour.

Lusby et al. (2013) relaxed the constraints CT6 and CT12, leading to a MILP formulation with binaries only for the maintenance decisions [33]. The production time steps were aggregated weekly. The stochastic scenarios were not aggregated, leading to 2-stage stochastic programming structure solved by Bender’s decomposition. The master problem concerns the maintenance decisions and the refueling quantities, independent sub-problems are defined for each stochastic scenario with continuous variables for power productions and fuel levels. The heuristic of [33] computes first the Linear Programming (LP) relaxation exactly with the Bender’s decomposition algorithm. Then, a Branch&Bound (B&B) phase repairs integrity, branching on binaries without adding new Bender’s cuts. Actually, a pool of solution is generated to have several solutions to repair with respect to the constraints CT6-CT12. Using a pool of solution improves feasibility issues and provides better solutions than using only the best solution. The resulting heuristic approach was efficient for the small dataset of the qualification, difficulties occur for the final instances of the competition, especially B8 and B9 instances.

Rozenkopf et al. (2013) ([41]) considered an exact formulation of CT6 constraints, in a Column Generation (CG) approach dualizing coupling constraints among units, i.e. CT1 demands and CT14 to CT21 scheduling constraints. The MILP considered a unique scenario, the average one, with production time steps aggregated weekly. The CG approach computes a LP relaxation relaxing scheduling constraints CT14-CT21, with sub-problems solved by dynamic programming. The CG heuristic adds these scheduling constraints in the ILP with the columns generated. The ILP solution furnishes the final maintenance dates and week where the CT6 constraints are activated, the 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.

We note that a compact MILP formulation for the full problem was provided in [28], introducing binaries for all time steps, cycle and production mode $m\in{\mathcal{M}}$. 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 CG work similar with [41].

3.3 Heuristics based on the 2-stage decomposition

A very common idea was to decompose the problem in a 2-stage structure, distinguishing the optimization of maintenance and refueling of T2 units from the power production problems, as in [2, 7, 8, 23, 24, 29]. The first-level problem fixes the maintenance dates and the refueling levels modifying the current solution while respecting constraints CT7-CT11 and CT13-CT21. The current planning is modified in a neighborhood search, where the feasibility of constraints CT13-CT21 can be guaranteed using MILP of Constraint Programming models. The second-level sub-problems compute independently for each 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. Such decomposition approaches in competition were not among the most efficient. We note that the operational approach in EDF was in this scope before the Challenge, we refer to [30].

3.4 Local search approaches

Local search approaches, tackling the problem without decomposition, were used by the two best approaches of the Challenge, having similar results, on one hand a methodology similar with LocalSolver [22, 4] and on the other hand an unpublished Simulated Annealing local search [36]. We note that excellent results are reported in [22] with a time limit of few minutes, which is interesting for an industrial application. The choices of the neighborhoods is crucial for such approaches to have efficient moves and cost computation of the locally modified solution, to be able to explore aggressively a lot of solutions. The neighborhoods chosen by [22] are reported in their paper:

• •

k-MoveOutagesRandom: select k outages among T2 plants randomly and move them randomly in “feasible” time intervals, 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 scheduling constraints 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.

One can give explanations to the better quality solutions of local search approaches. Many exact optimization and decomposition heuristics restricted the model with heuristic pre-processing, it induce over-costs, which also explains the better results of [22, 36] that avoided such over-costs. Furthermore, heuristic decomposition and exact approaches are slowed down by the use of MILP, LP or CP computations. Local search approaches explore aggressively more solutions. Furthermore, the neighborhoods used for decomposition heuristics were usually small (e.g allowing to move only one an outage from only one week). Such neighborhoods were used by [36] in the qualification phase, and led to bad solutions. The effort to have diversified and large neighborhoods, as mentioned in [22], is crucial.

3.5 Other works after the Challenge

After the Challenge, other issues were investigated. One open question was to compute efficiently dual bounds. The first dual bounds were furnished by[8] using dual heuristics, and were improved in [18, 16]. Semi-definite relaxations were investigated, to improve bounds obtained with LP relaxations, but the size of the instances is a bottleneck for such approaches [25]. Using the optimal control theory, one can consider the fuel and production constraints as a discretized Partial Derivative Equations system. It raised new perspectives for local search approaches: a sensitivity analysis can indicate moving directions, as in [3]. However, the method do not allow moves that swap maintenance orders, which is interesting among the nuclear reactors from the same site with constraints CT14, as used by [22]. Mathematical Programming methodologies were investigated to complexify the model, searching for a robust scheduling considering delays in the maintenance operations. Such extension were studied with CG approaches [37, 26], and also with multi-objective matheuristics [15, 26]. For these last studies, computations are tractable restricting computations to horizons of 2-3 years.

4 MILP formulation

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

4.1 Definition of variables

The definition of variables is highlighted in Table 4. We define variables $d_{i,k,w}$ for maintenance dates for each cycle $i,k$ with $d_{i,k,w}=1$ if and only if the outage $i,k$ begins before week $w$. These are the only integer variables thanks to the relaxation of CT6 and CT12, as in [33].

Other continuous variables are introduced: refueling quantities $r_{i,k}$ for each outage $(i,k)$, T2 power productions $p_{i,k,w}$ at cycle $k$, fuel stocks at the beginning of campaign $(i,k)$ (resp at the end) $x_{i,k}^{i},x_{i,k}^{f}$, T1 power productions $p_{j,w}$, and fuel stock $x_{i}^{f}$ at the end of the optimizing horizon. We note that T2 power productions $p_{i,k,w}$ are duplicated for all cycle $k$ to have a linear model, $p_{i,k,w}=0$ if week $w$ is not included in the production cycle $k$.

4.2 MILP formulation

The variable definition induces the MILP formulation below.

[TABLE]

(2) is required with definition of variables $d$. (3) and (4) model CT13 time windows constraints: outage $(i,k)$ is operated between weeks $\mathbf{To}_{i,k}$ and $\mathbf{Ta}_{i,k}$. (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 $i,k$ is not operated, ie $d_{i,k,W}=0$. (9) writes CT8 initial fuel stock. (10) writes CT9 fuel consumption constraints on stock variables of cycles $k$ $x_{i,k}^{i},x_{i,k}^{f}$. (11) models CT10 fuel losses at refueling. (12) writes CT11 bounds on fuel stock levels on variables $x_{i,k}^{i}$ thanks to the decreasing law (10). (13) models CT11 min 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, ie if $d_{i,k+1,W}=1$, which enforces to have disjonctive constraints where case $d_{i,k+1,W}=0$ implies a trivial constraints thanks to (12). (14) linearizes the constraints to enforce $x_{i}^{f}$ to be the fuel stock at the end of the time horizon. $x_{i}^{f}$ is indeed the $x_{i,k}^{f}$ such that $d_{i,k,W}=1$ and $d_{i,k+1,W}=0$, for the disjonctive constraints (14) that write a trivial constraints in the other cases thanks to (12), we define $\overline{S}_{i}=\max_{k}\mathbf{S}_{i,k}$. (15) is a common framework for scheduling constraints from CT14 to CT21, which was noticed independently in [33, 29].

4.3 Adding stability cost for dynamic re-optimization

The objective function (1) considers only the financial cost similarly with [38]. With the extension of stability cost defined in Section 2.2, we consider two objectives, the financial cost $\mathbf{C}^{roadef}$ and the desorganization cost $\mathbf{C}^{desorg}$:

[TABLE]

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 $\mathbf{C}^{pen}_{i,k,w}$ is defined by $\mathbf{C}^{pen}_{i,k,w}=1$ for all $w\neq W_{i,k}^{0}$. and $\mathbf{C}^{pen}_{i,k,W_{i,k}^{0}}=0$, i.e. $\mathbf{C}^{desorg}$ counts the number of desorganization. This makes the objective $\mathbf{C}^{desorg}$ integer with a granularity of $1$, which is interesting considering an $\varepsilon$-constraint method with $\varepsilon=1$ ([32]). In this case, the constraint to bound the number of modifications of the original planning by at most $\mathbf{Nmax}$ is given by:

[TABLE]

5 Constructive matheuristics

5.1 Exact and heuristic pre-processing

As noticed during the Challenge, heuristic and exact pre-processing may reduce the size of the problem to improve the solving capabilities. Exact pre-processing strategies fix necessarily optimal decisions, whereas heuristic strategies reduce the problem with approximations inducing possibly over-costs.

Time window constraints CT13 can be tightened in a exact pre-processing. initial stocks defined in CT8, maximal T2 productions CT5, max stock before refueling with CT11. This tightens by inference the time window constraints, as described in [10, 16, 33]. Some outages can also be removed when their earliest completion time exceed the time horizon.

There is no maximal length on a production cycle in our model and in [38]: decreasing profile phases can be followed by a nil production phase without any maximal duration in this state. However, such situations are not economically profitable as the T2 production cost is lower than the T1 production cost. Thus, [33] avoided such situations, defining heuristically maximal duration for production cycles and tightening by inference the latest dates to begin outages $\mathbf{{Ta}}_{i,k}$. This is a heuristic pre-processing, the reduction of the solution space can be useful for exact methods. Similarly with [17], another heuristic pre-processing fixes variables guided by the LP relaxation. A LP-based heuristic reduction can be applied tightening the time window constraints with the integral values in the LP relaxation of variables $d_{i,k,w}$.

5.2 Post-processing to repair partial solutions

Another application of the stability extension is to repair an initial solution $W^{1}(i,k)$ of maintenance planning, that may be infeasible, looking similarly with Feasibility Pump (FP, [5]) for the closest feasible solution using $\mathbf{C}^{pen}_{i,k,w}=|w-W_{i,k}^{1}|$ or $\mathbf{C}^{pen}_{i,k,w}=(w-W_{i,k}^{1})^{2}$ and minimizing $\mathbf{C}^{desorg}$ or $\mathbf{C}^{desorg}+\epsilon\mathbf{C}^{roadef}$ with $\epsilon>0$ small. Furthermore, one can also restrict the time window possibilities around the baseline solution to avoid having a too large model

5.3 A simplified MILP formulation

One can approximate the previous MILP formulation, assuming that the T2 production is either null during outages, either at its maximal level, ie $p_{i,k,w}=\mathbf{\overline{P}}_{i,w}(d_{i,k,w-{\mathbf{Da}_{i,k}}}-d_{i,k+1,w})$. Imposing such constraints can lead to infeasibilities, it requires some other relaxations. T2 power can be upper than some low demands $\mathbf{D}^{w}$, leading infeasibilities with 5. Relaxing 5 with minimal production inequalities fix such problems. Furthermore, a maximal T2 production can be infeasible with the fuel constraints. More specifically, a maximal T2 production (and thus a maximal fuel consumption), can be in contradiction with CT13 earliest dates of outages and positivity of fuel levels To avoid such infeasibilities, we introduce continuous variables $\Delta_{i,k}$ to authorize negative fuel stocks in soft constraints with high values $C_{i,k}^{\Delta}$ in the objective function. It leads to following simplified MILP formulation:

[TABLE]

This MILP formulation computes a first maintenance planning. Then, the induced production costs are computed as an LP, ensuring feasibility of all the constraints. If this high level planning cannot lead to feasible production and fuel levels, a local reparation is required, as described in 5.2.

5.4 Restrict- Relax-and-Fix algorithm

Relax-and-Fix (RRF) algorithms ([27]) construct iteratively a solution following a partition of the variables, fixing values from the previous computations, and relaxing continuously others. A RRF strategy can partition the variables following the cycle index $k$. It computes successively the solutions for cycles $k$, the dates of cycles $k^{\prime}<k$ are fixed by the previous optimizations, and the variables related to the cycles $k^{\prime}>k$ are relaxed continuously or simplified using a parametric aggregation [18].

Feasibility difficulties occur when a cycle optimization fixes outages and no feasible solution exist to place the next outages because of constraints CT13-CT21 and/or fuel constraints CT11. This induces in Algorithm 1 to consider binaries for cycle $k$ and $k+1$ for the iteration placing outages $k$, and to relax continuously/simplify formulations only for cycles $k^{\prime}>k+1$. The remaining fuel costs are relaxed also to focus on the local optimization around the cycles to optimize, aggregating the cycles $k>k_{0}$ in the cycle $k_{0}+1$ relaxing the outage and the fuel constraints, as in [18]. The MILP formulation considered for each iteration is thus:

[TABLE]

5.5 Construct, Merge, Solve & Adapt

Another decomposition iterates unit by unit, fixing in the current planning of the computed sub-problems, declaring variables as integer for the current unit to optimize, and relaxing continuously the other units. Such heuristic has several weaknesses. Firstly, having one infeasibible sub-problem stops the construction. Secondly, the choice of the iterated unit order has some influence in the final solution.

In Algorithm 2, a Construct, Merge, Solve & Adapt (CMSA, [6]) strategy computes independently (for a possible parallel computation) planning for all units, as if the other units were T1 power plants. Merging these planning is likely unfeasible for the scheduling constraints CT14-CT21, a second phase repair the feasibility around this ”ideal” 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.

6.1 General algorithm

We present here the local search algorithm described in the Algorithm 3. The neighborhoods are defined using MILP definitions, similarly to [1, 19]. Many MILP neighborhoods can be designed, and large neighborhoods may be explored. Algorithm 3 alternates the choice of the neighborhood , similarly with multi neighborhood search approaches. It induces that a local optimum for Algorithm 3 is a local optimum for all the considered neighborhoods. The stopping criterion could be a maximal time limit or a maximal number of iterations, or being in a local extremum for all neighborhoods.

MILP neighborhoods may be explored for small sub-problems with optimal B&B solving, as in [31]. In this work, MILP neighborhoods are defined with three characteristics for an efficient B& B search in defined time limits. Parametrization is empiric, for a good trade-off between solution quality and time spent in MILP solving. The MILP neighborhoods are 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 an absolute or relative threshold between the best primal and dual bounds can avoid to waste time without improving significantly the best primal solutions.

• •

a **specific parametrization for MILP solving **: for an efficient B&B search in the defined time limit. Small neighborhoods need few specific MILP parametrization. For a better efficiency in short resolution time for difficult sub MILPs, one can emphasize the heuristic search , disable or limit cutting plane passes and MILP preprocessing (with parameters mipemphasis, cutpass, presol 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 MILP resolution is at least as good as the current solution at each iteration. This algorithm is thus a steepest descent algorithm, similarly to MILP-VND as in [31].

6.2 MILP neighborhoods

Multiple types of large and variable neighborhoods can be defined. We denote with ${W}_{i,k}^{o}$ the beginning week of outage $(i,k)$ in the current incumbent solution. Neighborhoods define strategies to fix variable $d_{i,k,w}$ around the solution defined by ${W}_{i,k}^{o}$. Generic neighborhoods are first derived from [9, 20]:

• •

${\mathcal{N}}^{rins}$ : Similarly to RINS heuristic [9], variables are fixed if they are a common integer value in the LP relaxation and in the current solution.

• •

${\mathcal{N}}^{LB}_{k}$: for k an integer, k-Local Branching neighborhood allows only $k$ modifications to the incumbent, adding constraints as in [20]:

[TABLE]

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

• •

${\mathcal{N}}_{I}^{units}$ unit selection: only T2 units $i\in I\subset{\mathcal{I}}$ are re-optimized.

• •

${\mathcal{N}}_{k,k^{\prime}}^{cycles}$ : variables relative cycles $k^{\prime\prime}$ with $k\leqslant k^{\prime\prime}\leqslant k^{\prime}$ are re-optimized.

• •

${\mathcal{N}}_{(a,b)}^{TW}$: outages are re-optimized in time windows $[{W}_{i,k}^{o}-ak-b,{W}_{i,k}^{o}+ak+b]$.

${\mathcal{N}}_{(0,b)}^{TW}$ is a re-optimization with a fixed radius $b$ around the incumbent, ${\mathcal{N}}_{(0,1)}^{TW}$ is larger than neighborhood widely used for the Challenge, moving one outage of one week. The linear dependence with $a>0$ 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. ${\mathcal{N}}_{k,k^{\prime}}^{cycles}$ is a first answer to have complementary neighborhoods to swap the order of outages. ${\mathcal{N}}_{I}^{units}$ can be implemented for all single units, in this case, it authorizes large neighborhoods, and the 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 when the substituted production cost is higher.

6.3 Sequences of neighborhoods

A key point in the Algorithm 3 is the sequence of neighborhoods. Applying Algorithm 3 with a single type of neighborhoods, it allows to analyze the quality and number of local extrema. In this case, different neighborhoods are compared starting with the same initial solutions. Usually, VND local searches increase the size of the neighborhoods when a local minimum is reached. It can be implemented for Local Branching neighborhoods ${\mathcal{N}}^{LB}_{k}$ with $k$ increasing , as developed in [35] or with neighborhoods ${\mathcal{N}}_{(0,k)}^{TW}$. In our case study, long computation time are allowed in the operational process for the monthly optimization of the maintenance planning. We do not seek to find the best solutions in short resolution time. 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 Algorithm 3 when no neighborhood improves the current best solution.

Using neighborhoods ${\mathcal{N}}_{I}^{units}$ or ${\mathcal{N}}_{k,k^{\prime}}^{cycles}$, one can design different partitions of the binary variables. A partition of units $P\subset{\mathcal{P}}({\mathcal{I}})$ give rise to neighborhoods ${\mathcal{N}}_{p}^{units}$ for $p\in P$ partitioning the variables of the problem. $P=\{\{i\}\}_{i\in{\mathcal{I}}}$ is a first partition units by units, reoptimizing all the planning separately. One can consider partitions with subsets of at least $k$ elements where $k$ is given. In the applications, partition occurs in sites with several nuclear reactors, in subsets of the constraints CT14-CT15. VND iterations can successively re-optimize along different partitions, similarly with POPMUSIC (Partial optimization meta-heuristic under special intensification conditions, see [42]). ${\mathcal{N}}_{k,k}^{cycles}$ neighborhoods define also a partition, ”orthogonal” to the unit re-optimization.

7 Computational results and analyses

The computational experiments were computed with a personal computer having a processor Intel(R) Core(TM) i7-4790S CPU @ 3.20GHz using up to 8 threads, running Linux Ubuntu 18.04 with 8 GB of RAM memory. MILP and LP were solved with Cplex version 12.8.0. We note that, contrary to the straightforward MILP solving with Cplex, the final matheuristic is even scalable on lighter hardware, with also older versions of Cplex.

7.1 New instances

We used the dataset from the EURO/ROADEF 2010 challenge presented in Section 2.5 and Table 3, aggregating production time steps to weeks and the stochastic scenarios to a single scenario of demands and costs as in [41]. Small instances were also generated to have comparisons with known optimal values. Instances suffixed with _3_120 consider at least $3$ cycles in the $120$ first weeks of the ROADEF instances. All these instances _3_120 are solvable in less than one hour, except B9_3_120, as shown in Table 7.

Instances are generated from the dataset B and X from the Challenge, respectively suffixed with Ext2 and Ext3, removing the time window for cycles $k>2$ and $k>3$ respectively. These instances were generated to reproduce instances with a similar difficulty from B8 and B9. Table 5 compares the number of binaries of these new instances to the original ones. Table 5 shows that the number of binaries is increasing with Ext2 and Ext3 datasets, Ext2 reduces considerably the gaps between B8 - B9 and the other instances regarding the number of binaries. This gap closing is also observable after Cplex solving, as shown in Table 8.

Table 6 compares the Best Known primal Solutions (BKS) of the extensions Ext2, Ext3 to the one of the original dataset for the instance sets B and X, except for B8 and B9. Table 6 shows that the heuristic reduction of time windows has an impact of $1\%$, which is very significant regarding the financial value, and also comparing to the over-cost induces by different optimization methods for the Challenge. The numerical results aim to provide an accurate resolution of the difficult instances, which are the most relevant for the industrial application, with a validation on $18$ small instances where the optimal values are known.

7.2 Table results

To compare different approaches on an instance $i$, we computed for each instance $i$ the cost of the best primal solution known, denoted $\emph{BKS}(i)$, without restriction in the time limit. Lower bounds (LB) or upper bounds (UB), denoted $v(i)$ on instance $i$, are compared with the gap indicator: ${\emph{gap}_{i}=\dfrac{\mid v(i)-\emph{BKS(i)}\mid}{\emph{BKS}(i)}}$. Over-cost gaps over $100\%$ are considered as non significant solutions, and considered as a failure.

To show results on a group of instances, we provide statistics using the following indicators:

• •

$\bar{t}$, the average time solving in seconds;

• •

$\bar{g}$, the average gap;

• •

$\sigma_{g}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\emph{gap}_{i}-\bar{g})^{2}}$ to measure dispersion among $n$ instances;

• •

Q1, Q2, Q3 denote the quartiles of the gap sequence, Q2 is the median;

• •

$N_{x}$ represent the number of instances with a gap below $x\%$;

• •

$N_{F}$ represent the number of failures, i.e. the number of instances with a gap upper $100\%$;

7.3 MILP solving characteristics

The characteristics of MILP solving is analyzed with Tables 7, 8 and 10.

Tables 7 and 8 compare the gaps to the BKS for different LB and UB at different stage of the B&B solving: LB of the LP relaxation, LB and UB at the root note after cutting plane passes and before branching, and B&B computation in time limits of 1h and 3h. Table 7 highlights that the B&B algorithm is very efficient on the 18 instances where the optimal value is known, with high quality LB and UB of the B&B solving at the root node.

On the most difficult instances, the UB and the MILP primal heuristics may be of a poor quality, as shown by Table 8. Actually, the LP relaxations are not useful to compute primal solutions for instances like 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 induces difficulties for MILP primal heuristics, relying on the LP relaxation such as Feasibility Pump, RINS or Local Branching, and tailored variable fixing heuristics. This explains two bottlenecks for straightforward MILP solving in the instances with large time windows shown in Table 8: the LP relation is significantly of a worse quality and the search of primal solutions is inefficient, compared to the easier instances with smaller time windows. On one hand, Such difficulties have to be overcome with the matheuristics. On the other hand, matheuristics can rely on the efficiency of B&B search for restricted sizes of instances.

We note that the LB obtained with the straightforward B&B algorithm and with the dual heuristics of [18], are of a good quality, these dual bounds are reported and denoted ”Best LB” in Table 9 for the instances where the optimal value is not proven by the straightforward B&B computations. These dual bounds improve significantly the ones from the straightforward MILP solving, which allows relevant comparisons with the primal heuristics.

Table 10 compares the termination time of B&B to optimality on the smallest instances: without and with the exact pre-processing of [16], the influence of warmstarting the B&B algorithm with the optimal (or best known) solution. The exact pre-processing of [16], tightening time windows thanks to constraints CT5 and CT8-11, has an strong influence in the computation times to optimality. Warmstarting has less influence for the small instances. Warmstarting is relevant when the search of primal solutions is inefficient, which is the bottleneck of MILP solving noticed for the most difficult instances. For the MILP computations inside the MIP-VND matheuristic, warmstarting is also relevant, especially for the largest neighborhoods.

7.4 Optimizing around an initial solution

In the instances from the ROADEF Challenge, no initial solution was given. To analyze the impact of having an initial solution and minimizing with stability costs, we considered as initial solutions the ones given by the matheuristic of Section 5.1 without the repairing post-processing. Such initial solutions are not always feasible, and are of a medium quality.

Adding penalization to a baseline solution in the objective function changes the MILP convergence characteristics. Guiding the solution search around the baseline solution improves very significantly the solving time, as shown in Table 10. The higher are the stability coefficients $\mathbf{C}^{pen}_{i,k,w}$, the more the B&B algorithm is accelerated. Many solutions with similar ROADEF costs exist which was already noticed for the Challenge, especially in [22]. Similarly with symmetries in the B&B tree search, it induces that few cut-off are processed, and the optimal computations explore many nodes where primal solutions with similar costs are found. On the contrary, the tree search is guided around the LP relaxation with high stability costs, with fewer nodes in the B&B tree to explore. As is, it can be used to repair efficiently infeasible solutions.

These characteristics allow to consider an $\varepsilon$-constraint method to draw the Pareto curve of non dominated solutions in the bi-objective optimization cost/stability ([32]). The granularity for the $\varepsilon$-constraint method can be set to $\varepsilon=1$ without granularity errors, using the integer objective functions for stability penalization. Figure 1 presents such Pareto frontiers for instances B7-3-120 and X11-3-120. 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. The price of stability is higher than the price of robustness shown in [15].

7.5 Constructive matheuristics

With the results of Table 8, the first challenging issue is to compute accurately primal solutions for the most difficult instances. The numerical results of constructive heuristics are highlighted in Table 9. Simplified MILP computations of section 5.1, the CMSA and the R-R&F heuristics succeeded to provide feasible solutions for all types of instances, without failures, contrary to straightforward Cplex solving in one or three hours. This was the first motivation to design these primal heuristics.

The MILP formulation of section 5.1 is easier for the MILP primal heuristics, it accelerates significantly the search of primal solutions. For the difficult instances like B8 and B9, feasible solutions were found in 15 minutes including the repairing post-processing, as shown in Table 9. The CMSA algorithm requires a similar computation time than the simplified MILP, with a similar (slightly better) quality of solutions. R-R&F approach gives significantly better solutions, but requiring longer computation times to converge to a feasible solution.

We note that the computation times from 15 minutes to two hours for the constructive heuristics were not optimized, the objective was to solve accurately the instances without over-fitting. Contrary to the R-R&F algorithm where the time limit in MIP solving is a bottleneck, the computation times can be decreased for the CMSA and the simplified MILP approaches. Such computation times are reasonable for the industrial application as a first phase before local search improvements. At this stage, the goal was to be able to compute accurately feasible solutions for all the instances, it is reached. The VND matheuristic can be initialized with feasible solutions from the three types of constructive matheuristics.

7.6 POPMUSIC-VND matheuristics

Table 13 compares the quality of the local minimum considering each type of single neighborhoods separately in the MILP-VND local search.

Neighborhoods ${\mathcal{N}}_{(0,1)}^{TW}$ include the neighborhood to move one single outage from one week, considered by many approaches for the ROADEF Challenge. ${\mathcal{N}}_{(0,1)}^{TW}$ includes also 1-translation neighborhoods, translating from one week the maintenance planning. Table 13 shows that it induces a poor quality of local minimums. Many local minimums exist with neighborhoods ${\mathcal{N}}_{(0,1)}^{TW}$ and also for the smaller neighborhoods. The quality of local minimums is improved increasing the radius with ${\mathcal{N}}_{(0,k)}^{TW}$ neighborhoods. However, significant over-costs remain. It justifies to consider other types of neighborhoods in a local search method for the problem. We note that “funnel” neighborhoods ${\mathcal{N}}_{(a,b)}^{TW}$ with $a>0$ are interesting instead of using only ${\mathcal{N}}_{(0,b)}^{TW}$ neighborhoods.

Cycle neighborhoods ${\mathcal{N}}_{k,k^{\prime}}^{cycles}$ allow to swap the order of outages, which is unfeasible with small ${\mathcal{N}}_{(a,b)}^{TW}$ neighborhoods with spacing constraints CT14-CT15. Significant improvements are observed in Table 13, but it requires more computation time than the other neighborhoods. Swap moves are also possible with ${\mathcal{N}}_{I}^{units}$ neighborhoods, when subset $I$ contains units in a site defined by constraints CT14-CT15. Local minimums are of an excellent quality considering unit re-optimization with $|{\mathcal{I}}|=1$, with MILP sub-problems that converges very quickly (order of one second).

Table 9 shows the final results considering in the POPMUSIC-VND neighborhoods ${\mathcal{N}}_{(1,3)}^{TW}$, ${\mathcal{N}}_{k,k}^{cycles}$ and ${\mathcal{N}}_{I}^{units}$ with $|{\mathcal{I}}|=5$. The VND local search starting with the CMSA solutions is very efficient, improving quickly primal solutions and requiring few VND iterations to compute a local minimum for all neighborhoods. The resulting matheuristic is more efficient in the ratio improvement of solution/computation times than the R-R&F heuristic. Furthermore, comparing the local extrema given by the VND to the optimal solutions of small instances, we had only one example where the VND induces over costs more than the tolerance gap of Cplex parametrized to $0,01\%$. The tolerance gap induces sometimes that the VND solutions were better than the “optimal” solutions with the terminantion of the B&B algorithm, within the $0,01\%$ tolerance gap.

This illustrates that VND-MILP allows to define easily neighborhoods induced by operational expertise. The genericity of MILP modeling and resolution avoid specific implementations. An application of MILP-VND is to design and test neighborhoods with few implementation efforts, to select neighborhoods that should be implemented properly in a local search as in [36, 22]. These results give some insights to explain some results of local search heuristics used for the Challenge.

7.7 Explanations of some results of the Challenge

This section give some insights provided by our intermediate results to analyze the design of previous approaches, in mathematical modelling, and also for the operators of local search heuristics.

7.7.1 MILP formulations and efficiency of B&B solving

Our formulation is similar to deterministic version of the MILP formulation in [33]. The major difference is in the definition of binary variables $d_{i,k,w}$: binaries $x_{i,k,w}$ defined in [33] are equal to $1$ if and only if outage beginning week for cycle $(i,k)$ is exactly $w$. Actually, these definitions can be seen as isomorphic (as in the analysis of [14]) with $x_{i,k,w}=d_{i,k,w}-d_{i,k,w-1}$ and $d_{i,k,w}=\sum_{w^{\prime}\leqslant w}x_{i,k,w^{\prime}}$. We note that this transformation applies also for the matheuristic approach [29]. If the MILP formulations are polyhedrally equivalent, it can have an influence in the efficiency of MILP solvers, and different activation of cuts with Cplex, as in[14]. This is not crucial in the master problem of the Bender’s decomposition in [33] or in the MILP model of scheduling constraint in [29] as this problems are quickly solved by Cplex. This had more impact for the exploration of our neighborhoods, we observed than our formulation provides better branching (well-balanced branching similarly with SOS1 branching), whereas the formulation of [33] was slightly better for the efficiency of the MILP primal heuristics.

Adding continuous variables with the production level makes the MILP solving with Cplex difficult, as shown in Table 8 and in the Appendix A. This justifies the Bender’s decomposition approach to tackle the problem with many stochastic scenarios, as in [33], this would be also a perspective to generate only the CT6 constraints that are useful. The difficulty met with an approximate description of CT6 constraints using only continuous variables makes the exact MILP formulation of [28] intractable even for small sizes of instances. Lastly, we mention that solving the LP with fixed maintenance dates, which corresponds to ${\mathcal{N}}_{(0,0)}^{TW}$ computations, are solved easily by Cplex (in less than one second). There is no need anymore for decomposition approaches like price decomposition, as used in ancient work with less efficient LP solvers and hardware like in [12].

7.7.2 Extended formulations for CT6 constraints?

Several approaches used extended formulation and column generation, as in [41, 26], to deal with an exact description of the CT6 ”stretch” constraints. The Appendix A give results that show that CT6 constraints have few impact in the quality of the final solutions: a little over-cost is given repairing the stretch constraints from an optimal solution relaxing the CT6 constraints, and very local modifications allow to decrease significantly the over-costs to negligible over-costs. Actually, this can be analyzed and explained numerically: stretch constraints are not often active, it is more costly to have a non maximal T2 production as it is cheaper. These situations can be unavoidable for feasibility and time windows constraints. In such cases, it concerns few weeks in the planning, and the local optimization is efficient to incorporate the CT6 constraints and improve the cost. We note that the column generation approaches are more limited than our matheuristics: [41] was not able to implement the Branch&Price resolution because of time limit constraints, with 40 minutes to solve the LP relaxation of the extended model, and [26] restricted the time horizon to 3 years instead of 5 in the Challenge. Willing to face a 5 year time horizon, with the large time windows like datasets Ext2 and Ext3, the capacity of Branch&Price resolution would be still limited nowadays.

7.7.3 2-stage decomposition approaches

The 2-stage decomposition is a natural approach, taking separately maintenance planning and production decisions having fixed the other decisions. For the production decisions, sub-problems are independent for all the stochastic scenarios like in the Bender’s decomposition approach [33]. Many heuristic approaches used this 2-stage decomposition: [2, 7, 8, 23, 24, 29]. Having the maintenance and refueling decisions fixed, it makes easier to compute the productions. For a maintenance planning, it computes the production, indicates if the planning is feasible and what is the cost, using a LP or a greedy algorithm. In our simplified model, it corresponds to ${\mathcal{N}}_{(0,0)}^{TW}$ computations, which are very fast. Modifying the maintenance decisions were in these approaches dealt by very local transformations: translating the maintenance planning of one weeks for consecutive outages of a T2 power plant, or modifying one outage by one week. Such neighborhoods are smaller than our ${\mathcal{N}}_{(0,1)}^{TW}$ neighborhoods, with many local optimums and a poor quality of local minima as shown in Table 13. This is an explanation of the over-costs given by approaches [2, 7, 8, 23, 24, 29].

Another difficulty of such approaches are that the CP or MILP computations are more time consuming than straightforward local search approaches, it induces to explore less solutions in the solving time. In our VND, the MILP explorations cover larger neighborhoods than just computing the costs. These MILP iterations are more costly than ${\mathcal{N}}_{(0,0)}^{TW}$ computations. The VND has the advantage of furnishing directions of optimizations, which helps the coordination phase to modify the maintenance planning. On the contrary, 2-stage decompositions modify the maintenance planning problem blindly with random moves from the current planning. This is another explanation of the limitations of such approaches, while the optimal exploration with ${\mathcal{N}}_{(0,1)}^{TW}$ is more efficient to accomplish a similar task. We note that motivations of the work [3] was also to give directions of re-optimization for a 2-stage decomposition approach.

7.7.4 Local search without decomposition

The analyses of neighborhoods with VND shows that the “Winter Jump” neighborhoods are crucial to avoid being struggled in local minimum, as shown in Table 13. The neighborhoods ”k-MoveOutagesConsecutive” in [22], as mentioned in section 3.4, incorporate these winter jumps. We note that the ”k-MoveOutagesRandom” and ”k-MoveOutagesConstrained” in [22] contains similar neighborhoods than our ${\mathcal{N}}_{(a,b)}^{TW}$ and ${\mathcal{N}}_{k,k^{\prime}}^{cycles}$ neighborhoods. The benefit of combining these three types of neighborhoods may have a similar influence than our Table 13 to have local minimums of an excellent quality. Such choice of neighborhoods is a key component to design a local search heuristic for this problem. This result was analyzed with matheuristics, the implementation of [22] allows to have quicker moves and thus a more scalable approach. The approach of [22] is able to explore aggressively a lot of solutions, an is efficient also within 15 minutes of resolution time, contrary to our matheuritics which require heavier computations.

We note that [36] did not publish elements on the final design of neighborhoods. They mentioned having used only neighborhoods which are smaller than ${\mathcal{N}}_{(1,0)}^{TW}$ in the first phase, like the 2-stage decomposition approaches. The improvements in the quality results of [36] between the qualification phase and the final result are due to the choice of other neighborhoods. The over-costs was in average of $2,8\%$ for the qualification phase, whereas they had 6 best solutions and very few over-cost in general for the final results. This is another proof of the crucial impact of the design of neighborhoods for this problem. It explained also the worst quality of solutions obtained by 2-stage decomposition approaches.

7.8 Contributions for the industrial application

Our results give also contributions for the industrial application. A first result is the validation of the design of the heuristic of [22], which explains its quality. This was not obvious after the Challenge why such approaches were the most efficient, and practitioners may mistrust optimization approaches without guarantee of optimality or dual bounds. We note that dual bounds are given now in [16], which validates the quality of the solutions furnished for the Challenge. With this work, one may trust the design of such heuristics to give good solutions for instances that are met in the operational application. We note that experts could have the intuition that ”winter jumps” may be en efficient neighborhood, the results of [22] and Table 13 justify it is a crucial point. The following of this section highlights the impact of intermediate results for the industrial application.

7.8.1 Relevancy of simplification hypotheses

Before the Challenge, the operational resolution of the problem in EDF contains the approximations of our model mentioned in section 2.1, we refer to [30, 26]. This paper contains even more constraints and less approximations, relaxation of constraints CT6,CT12 but also CT14-CT21 and also no fuel optimization were common assumptions. Another heuristic reduction was to tighten the constraints of time windows, to ease the computations (or make it possible in reasonable time). Looking to the instances of the Challenge in Table 3, the heuristic reduction apply to all the instances except B8-B9 (and thus the reproduced instances of this study). One reason motivating the Challenge, was to give the problems without the approximations and simplifications used in their homemade resolution, and to analyze if the approximations are relevant.

The aggregation of production time steps and stochastic scenarios were also used in [41]. It induced an overall over-cost of $0,5\%$ compared to the best solutions, mainly the one from [22] without such aggregation. This gap main include also over-costs due to the non termination of the Branch&Price algorithm, with hierarchic computations in [41]. In Table 6, the impact of the heuristic tightening of time windows in a pre-processing phase, can be evaluated to $1\%$, which is huge. Our matheuristics and the local search of [22] are still efficient with large time windows, the reduction is not required anymore for solving capabilities, and this hypothesis seems to have the highest impact in terms of over-costs. Furthermore, the impact of relaxing CT6 constraints is negligible, as analyzed in Appendix A and mentioned in section 7.7. This work gave the impact of simplification hypotheses, the assumptions taken in this paper are a very good compromise in minimizing the over-costs and minimizing the computational difficulties, which is useful for the operational application.

7.8.2 Towards dynamic optimization

The impact of the optimization with large time windows, as previously highlighted, has also incidences in the dynamic re-optimization process considering stability costs. The Figure 1 illustrated that the initial baseline solution can have a huge impact in the over-costs while searching for a solution close to the baseline solution. In these cases, the over-costs may be of $5\%$, whereas starting from a very good solution has less impact. By comparison, the over-costs of having more robust solution was in order of $0,1-0,5\%$ in [15]. For this applications, it is interesting to have a maximal flexibility to optimize the latest outages with $k\geqslant 3$, so that the stability considerations would start with a good solution for these outages when the considered outage have to be fixed for stability while term is approaching.

8 Conclusions and perspectives

This paper investigated matheuristic approaches to optimize refueling and maintenance scheduling of nuclear power plants. Instances with time window constraints only for the first maintenance outages induce difficulty for state-of-the art MILP solvers. Several constructive matheuristics use decompositions to have more accurately feasible solutions solving with smaller problems in heuristic decomposition schemes. Once feasible solutions are built, a POPMUSIC-VND matheuristic allows to have accurately solutions of an excellent quality, with standard computers in reasonable time limits. Using matheuristics allows also a genericity in the objective functions, to derive a bi-objective optimization to find the best compromises solutions between the planning stability and the financial costs.

Furthermore, intermediate results help to understand better the problem characteristics and some ranking after the EURO/ROADEF 2010 Challenge. This article justified some approximations made by [33, 41, 29] to have tractable MILP formulations to face the large size of real size instances: aggregation of production time steps to week, at least partial aggregation of stochastic scenarios. The relaxation of the constraints CT6 are also justified thanks to the matheuristics. The matheuristic methodology, from MILP formulation work to the design of constructive and local search approaches, allows to analyze the impact of neighborhoods in the quality of local minima. This analysis helps to understand why many decomposition approaches failed for the challenge, because of too short neighborhoods. Our matheuristic study explains why approaches developed in [22, 36] for the Challenge were more efficient: using no time consuming LP, MILP or CP computations, more time could be devoted to explore solutions, and the use of variable and larger neighborhoods than ${\mathcal{N}}_{(0,1)}^{TW}$ for a better exploration. Especially, the ”Winter Jump” neighborhoods are crucial.

These results offer several perspectives for the industrial application. If the design of approaches [22, 36] are validated by this study, perspectives to improve other approaches are discussed in details in Section 7.7. The impact of the reduced time windows is shown in this study to have an enormous impact,compared to the other simplification hypothesis. We note that the heuristic reduction is an hypothesis to ease computations, the solving capabilities of our matheuristics and the frontal local searches [22, 36] do not need such simplification to run in reasonable time. Furthermore, we note than fixing local optimization around the known solution in the dynamic optimization application with the monthly re-optimization has similar effects, our bi-objective optimization work recommends to care about it in the real-life application. Another perspective is to consider other objective functions: maximizing robustness to lessen financial risks, or sustainable development issues like minimizing nuclear wastes or minimizing carbon emissions in regards to financial over-costs. This works illustrates that many solutions with similar costs exist, which was known already after [22]. One may choose a maintenance planning among the best ones with the cost optimization given for the Challenge and looking for a better trade-off using several objective functions.

Appendix A: relaxation of constraints

In this appendix, we justify the relaxation of modulation constraints CT12, and decreasing profile constraints CT6. Such constraints were often relaxed in the approaches for the EURO/ROADEF 2010 Challenge, inducing more complexity and difficulty in the optimization modeling and solving. This appendix provides justifications, regarding the industrial application and some numerical analyses.

A.1 Modulation constraints CT12

In this paper, modulation constraints, denoted CT6 in the Challenge description, are relaxed. In the Challenge, it is formulated as a maximal volume of energy generated when a unit $i$ is not producing at the maximal power $\mathbf{\overline{P}}_{i,w}$ at time step $w$, the cumulated energy being counted for all cycle $k$ when the fuel level is superior to $\mathbf{Bo}_{i,k}$. When the fuel level is inferior to $\mathbf{Bo}_{i,k}$, the production is at the maximal power defined by the CT6 decreasing profile.

Actually, such modulation constraints are approximated and aggregated from technical constraints of nuclear power plants. The fine grain technical constraints are as following, similarly to [39, 14], and illustrated Figure 2 :

• •

Min-up/min-down constraints: every unit $u$ has a minimum up time $\Delta^{on}_{u}$ online and a minimum down time $\Delta^{off}_{u}$ offline.

• •

Max duration with a non maximal power constraints: if a nuclear power plant is online, there is a maximal duration to have modulated power, i.e. non maximal power.

• •

Min duration between consecutive modulations: between consecutive modulations, there is a minimal duration to have maximal power.

• •

Max number of modulation constraints: for each day, there is a maximal number of periods where the production is not maximal.

If a T2 unit is online, the three last constraints induce a maximal volume of modulated power. Actually, the CT12 modulation constraints consider that the T2 units are always on-line in a production cycle. This last assumption has its justification as T2 production is less expensive than the T1 production. However, modulation may be necessary with low demands (lower than the available T2 production). In such cases, T2 units can be set offline with the min-up/min-down constraints, defining a volume of modulation that is not bounded. In these cases, there will be minimal durations offline with Min-up/min-down constraints. These extreme modulations happen in Figure 5 using the model of section 4.

A.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 $i\in{\mathcal{I}}$ is inferior to the level $\mathbf{Bo}_{i,k}$ in the cycle $k\in K_{i}$. The production of $i$ is deterministically defined, following a decreasing profile, piece-wise linear function of the stock level, as in Figure 4. More precisely:

CT6, ”stretch” constraints: During every time step $w\in{\mathcal{W}}$ of the production campaign of cycle $k\in K_{i}$, if the current fuel stock of plant $i\in{\mathcal{I}}$ is inferior to the level $\mathbf{Bo}_{i,k}$ , the production of $i$ is deterministic, following a decreasing profile, piece-wise linear function of the stock level, as in Figure 4. $m\in{\mathcal{M}}_{i,k}=[\![1,\mathbf{Np}_{i,k}]\!]$ define the points of the CT6 profile ending the cycle $(i,k)$. The points are denoted $(\mathbf{f}_{i,k,m},\mathbf{c}_{i,k,m})$ for each cycle $(i,k)$, $\mathbf{f}_{i,k,m}$ denote the fuel levels, $\mathbf{c}_{i,k,m}\in[0,1]$ is defined as the loss coefficient in power compared to the maximal power.

We will consider a lighter constraint , imposing only the UB constraints as illustrated in Figure 4. This MILP formulation was already introduced in [16]. New fuel variables $x_{i,w}$ are introduced, representing the residual fuel at week $w$ for T2 unit $i$. These variables require the following constraints to enforce $x_{i,w}$ to be the residual fuel:

[TABLE]

where $M_{i}=\max_{k}\mathbf{\overline{S}}_{i,k}$, such that $M_{i}$ verifies $x_{i,t}\leqslant M_{i}$. Indeed, if $d_{i,k,w}-d_{i,k-1,w}=1$, week $w$ happens in cycle $k$, the active constraint is $x_{i,t}\leqslant x_{i,k,s}^{i}-\sum_{w^{\prime}\leqslant w}\mathbf{F}_{w^{\prime}}\>p_{i,k,w^{\prime}}$, otherwise we have $x_{i,t}\leqslant M_{i}$, which is always true. The UB of the CT6 constraints, as illustrated Figure 4, are given with:

[TABLE]

If ( $\mathbf{c}_{i,k,m},\mathbf{f}_{i,k,m-1}$) do not depend on indexes $k$, we can have an equivalent MILP formulation using aggregated constraints, writing the constraints for the global production power $\sum_{k}{p_{i,k,w}}$:

[TABLE]

Adding the lighter CT6 constraints, it makes MILP solving difficult. Computations were not possible for the datasets B and X, inducing memory errors. We note that the aggregated formulation (35) provided improvements in the computation times compared to (34), but the formulation is still untractable for the real size instances. It can be explained in terms of number of variables, the extension adds $I\times W$ continuous variables $x_{i,w}$, whereas $I\times W\times K$ continuous variables were already in the model with T2 productions $p_{i,k,w}$. The main difference appears in the number of constraints, there were mainly $I\times W\times K$ constraints in the MILP of Section 4 with constraints (7), (33-35) require to add $I\times W\times(K+N_{p})$ constraints, whereas (34-33) require to add $I\times W\times(N_{p}+1)\times K$ constraints. This slows down very significantly the LP computations, when computable.

Adding CT6 constraints (33-34) or (33-35) has furthermore few impact in the quality of the LP relaxation, which was already noticed in [16]. We analyze here whether it is insightful in searching primal solutions. Actually, decreasing profile phases are rarely activated in our model. An explanation is that the upper bounds restrict the T2 production capacities. The T2 units having the lowest marginal costs of production, the optimization avoids the over-costs of production with mainly maximal T2 productions. Maximizing the T2 production available tends to avoid the decreasing profile phases, preferring to have earlier outages to produce at the maximum power in the production cycles. Optimal uses of decreasing profile occur when it is imposed by time windows constraints to have longer cycles than the duration to reach $\mathbf{Bo}_{i,k}$ fuel levels, which is the situation met in Figure 5. In such cases, the optimization tends to produce at the upper bound of production thanks to the difference in marginal costs between T1 and T2 units, which justifies the lighter formulation of CT6 constraints.

The over-costs to project the optimal (or best known) solutions computed thanks to the MILP of section 4.2 into a MILP with stretch constraints is given for the small instances in Table 14. Little over-costs are observed, around $0.3\%$ in average. The full relaxation of CT6 constraints allows to generate solutions of very good quality after repairing of CT6. Figure 5 suggests that some local modifications around the CT6 situations can improve the previous projected solution. Table 14 gievs also the improvements of the projected solutions with only two iterations of local search with MILP neighborhoods ${\mathcal{N}}_{(0,2)}^{TW}$ using (34), which are tractable computations. These very local modifications around the CT6 situations improve significantly the projected costs.

These conclusions justify the relaxation of CT6 constraints in a matheuristic approach, as developed in this paper, but also in the approaches of [33].

Appendix B: Intermediate results