%\documentclass[a4paper,11pt]{article} % \documentclass{svjour3} % onecolumn (standard format) %\documentclass[smallcondensed]{svjour3} % onecolumn (ditto) %\documentclass[smallextended]{svjour3} % onecolumn (second format) %\documentclass[twocolumn]{svjour3} % twocolumn % \smartqed % flush right qed marks, e.g. at end of proof \usepackage{amssymb} \usepackage{amsmath} \usepackage{mathrsfs} \usepackage{subeqnarray} \usepackage{graphics} \usepackage{tikz} \usepackage{algorithmic,algorithm} \usepackage{pstricks} \usepackage{a4wide} \usepackage{verbatim} \usepackage{graphicx} \usepackage{url} \usepackage{mathptmx} \RequirePackage{fix-cm} %\renewcommand{\baselinestretch}{1.1} \newtheorem{thm}{Theorem} \newtheorem{cor}[thm]{Corollary} \newtheorem{lem}{Lemma} \newtheorem{prop}[thm]{Proposition} \newtheorem{defn}[thm]{Definition} \newtheorem{rem}[thm]{Remark} % ----------------------------------------------------------------------------- \newcommand{\norm}[1]{\left\Vert#1\right\Vert} \newcommand{\abs}[1]{\left\vert#1\right\vert} \newcommand{\set}[1]{\left\{#1\right\}} \newcommand{\Real}{\mathbb R} \newcommand{\eps}{\varepsilon} \newcommand{\To}{\longrightarrow} \newcommand{\mb}[1]{\mathbf{#1}} \newcommand{\gb}[1]{\text{\boldmath$#1$}} \newcommand{\ov}[1]{\overline{#1}} \newcommand{\mml}{\{} \newcommand{\mmr}{\}} \newcommand{\lh}[1]{\widehat{#1}} \newcommand{\pname}[1]{\textsc{#1}} \newcommand{\lc}[1]{\mathscr{#1}} \newcommand{\lcs}[1]{\ov{\mathscr{#1}}} \newcommand{\lch}[1]{\widehat{\mathscr{#1}}} % ----------------------------------------------------------------------------- \newcommand{\zub}{z(UB)} %%%%% Obj function \newcommand{\zf}[1]{z(#1)} \newcommand{\zfp}[1]{z'(#1)} \newcommand{\zfh}[1]{\hat{z}(#1)} \newcommand{\zft}[1]{\tilde{z}(#1)} \newcommand{\zfo}[1]{\overline{z}(#1)} % ----------------------------------------------------------------------------- %%%%% Nuclear Power Plant \newcommand{\npp}{J} %%%%% Thermal Power Plant \newcommand{\tpp}{I} %%%%% Scenario \newcommand{\sce}{S} %%%%% Fuel cost \newcommand{\fc}{r} %%%%% Thermal variable cost \newcommand{\vc}{c} %%%%% Set-partitioning production planning index \newcommand{\spi}{\ensuremath{\ell}} %%%%% Mathematical formulation \newcommand{\mf}{F} %%%%% Dual of the mathematical formulation \newcommand{\dmf}{D} %%%%% Set-partitioning formulation variables \newcommand{\spv}{x} %%%%% Set-partitioning planning \newcommand{\spr}{R} %%%%% Problem name \newcommand{\edf}{EDF} %%%%% Problem name \newcommand{\mip}{$MIP$} %%%%% Problem name \newcommand{\plu}{LP_R} %%%%% Problem name \newcommand{\pld}{LP_s} %%%%% Problem name \newcommand{\pla}{\ensuremath{LP_a}} %%%%% Time horizon \newcommand{\tmax}{{\ensuremath{|\timeline|}}} %%%%% Time steps \newcommand{\timeline}{{{T}}} \newcommand{\timesteps}[1]{{\ensuremath{\timeline_{#1}}}} %%%%% Week \newcommand{\week}{w} %%%%% Week of timestep \newcommand{\weekof}[1]{{\ensuremath{\week(#1)}}} %%%%% Weeks \newcommand{\weeks}{{\ensuremath{W}}} %%%%% Weeks of an arc \newcommand{\weeksof}[1]{{\ensuremath{\weeks(#1)}}} %%%%% Week max \newcommand{\wmax}{{\ensuremath{|\weeks|}}} %%%%% Arc Weight \newcommand{\weight}{{\ensuremath{c}}} %%%%% Path \newcommand{\arcs}[1]{{\ensuremath{\mathcal{A}_{#1}}}} % ----------------------------------------------------------------------------- \DeclareMathOperator*{\argmin}{\text{argmin}} \def \union{\mathop{\mathrm{\cup}}} \def \equi{\Longleftrightarrow} % Insert the name of "your journal" with \journalname{journal of scheduling} % \begin{document} \title{Solving electricity production planning by column generation} %\subtitle{Do you have a subtitle?\\ If so, write it here} %\titlerunning{Short form of title} % if too long for running head \author{Antoine Rozenknop \and Roberto Wolfler Calvo \and Laurent Alfandari \and Daniel Chemla \and Lucas L\'etocart } %\authorrunning{Short form of author list} % if too long for running head \institute{ Antoine Rozenknop \and Roberto Wolfler Calvo \and Daniel Chemla \and Lucas L\'etocart\at LIPN (UMR 7030) - Universit\'e Paris 13, 99 av. J-B Clement 93430 Villetaneuse\\ Tel.: +33-1-49403894\\ Fax: +33-1-48260712\\ \email{name.surname@lipn.fr} \and L. Alfandari \at ESSEC Business School, av B. Hirsch 95021 Cergy-Pontoise \\ \email{alfandari@essec.edu} \and Daniel Chemla \at CERMICS \'Ecole des Ponts ParisTech, 6-8 av. Blaise Pascal, Cit\'e Descartes - Champs sur Marne, 77455 Marne la VallĂ©e Cedex 2 \\ \email{daniel.chemla@enpc.fr} } \date{Received: date / Accepted: date} % The correct dates will be entered by the editor \maketitle \begin{abstract} This paper presents a heuristic method based on column generation for the EDF (Electricit\'e De France) electricity production planning problem proposed as subject of the ROADEF/EURO 2010 Challenge. This is to our knowledge the first-ranked method among those methods based on mathematical programming, and was ranked $4^{th}$ overall. The problem consists in determining a production plan over the whole time horizon for each thermal power plant of the French electricity company, and for nuclear plants, a schedule of plant outages which are necessary for refueling and maintenance operations. The average cost of the overall outage and production planning, computed over a set of demand scenarios, is to be minimized. The method proceeds in two stages. In the first stage, dates for outages are fixed once for all for each nuclear plant. Data are aggregated with a single average scenario and reduced time steps, and a set-partitioning reformulation of this aggregated problem is solved for fixing outage dates with a heuristic based on column generation. The pricing problem associated with each nuclear plant is a shortest path problem in an appropriately constructed graph. In the second stage, the reload level is determined at each date of an outage, considering now all scenarios. Finally, the production quantities between two outages are optimized for each plant and each scenario by solving independent linear programming problems. \end{abstract} {\bf Keywords}: Energy, electricity production, column generation, mathematical programming, heuristics. \allowdisplaybreaks \section{Introduction} Optimization in the energy sector is a large subject. With the raise of raw materials prices and the growing awareness for environment concerns, energy consumption must be considered in a less greedy way in order to save resources. For years, several researchers have worked on optimizing the energy industry. In \cite{NRJbook}, the authors gather such type of problems, one of which is the Electricity production planning. Electricity can be produced through different ways: nuclear and other thermal energies are the major part of production, but hydraulic energy or other renewable energies have sometimes a significant role (in Norway almost all the production is done using hydraulic energy). Producers can mix these different production resources in order to satisfy the energy demand while decreasing their costs and/or their impact on the environment. Electricity companies look for effective optimization tools since even small percentages of gain through optimization translate into large savings. The problem of optimizing energy production planning has several difficulties, one of which is that the exact amount of the future demand and the exact price of raw materials are unknown. This difficulty naturally increases with the length of the time horizon, which can be in some cases very large. A typical energy production planning problem with a short term horizon is the Unit Commitment Problem (UCP). In the UCP, several generators are available for producing energy. Costs and constraints are associated with generators depending on their type: nuclear and other thermal units, hydraulic energy, wind turbines. For example, in a short term horizon, the level of water in dams has to stay above some threshold. A forecast of the demand over a period of time is given. The objective is to find the best production assignment in order to satisfy the demand while minimizing production costs. Kazarlis {\em et al.} \cite{GenUCP} propose a genetic algorithm to solve the UCP. Claudia Sagastiz\`abal {\em et al.} \cite{FINARDI2005} present a Lagrangian relaxation method. Frangioni {\em et al.} in \cite{Frangioni2008316} extend Lagrangian relaxation methods to the UCP with ramping constraints, i.e., constraints that bound the difference of production of thermal units from one time step to the next. Claude Lemar\'echal {\em et al.} \cite{A-DubGonLem05} also studied UCP in the French context. In long-term planning, typical constraints on energy production relate to maintenance operations that have to be executed periodically on nuclear plants. No production can be planned during these maintenance periods. The problem of scheduling the shut down of nuclear power plants for refueling and maintenance (see for example \cite{EDF1}) has been proposed as subject of the EURO/ROADEF $2010$ challenge. EDF power generation facilities in France stand for a total of $98.8$ GW of installed capacity and mix thermal energy ($90$ \% in $2008$, among which $86$ \% is produced by nuclear power and the rest by coal, fuel oil and gas), hydraulic and other renewable energies. The research project consists in modeling the electricity production assets, finding an optimal schedule of outages for each nuclear plant over a given time horizon, and determining an optimal production plan to satisfy demand. Moreover, in long-term electricity production planning, numerous uncertainties have to be taken into account, such as demand variation, generation units availability, spot market prices, quantities that can be bought or sold. This uncertainty leads to consider multiple scenarios and asks for a robust outage schedule that has to be feasible for every scenario. The scheduling of nuclear plant outages has to comply with various safety constraints and limitations on resources which are necessary to perform refueling and maintenance operations. The production planning fixes the quantity of energy to be produced by each plant at each time step for each scenario, and is also subject to technical constraints. The objective of the problem is to minimize the expected production cost over all scenarios. The proposed approach is based on mathematical programming and problem decomposition. An other approach used for solving the challenge ROADEF, based on a different decomposition method has been proposed by Lubsy {\em et al.} \cite{roadef10Danois}. In this paper a column generation method (see \cite{CGbook} and \cite{CGsurvey} for surveys on this approach) is proposed to solve it. %This subject is focused on thermal power plants and especially on those that have to be periodically %shut down for refueling and maintenance, e.g. nuclear ones. %is latter problem that has been studied for years now for instance) was the starting idea %The method developed in this paper is based on mathematical programming and column generation (see \cite{CGsurvey} for a survey). Column generation is often used to solve large-scale optimization problems when the problem can be formulated as a mixed integer linear program with a huge number of variables, or \textit{columns} -- called the \textit{master problem} -- as it is the case in the EDF problem. Creating and keeping all variables is not possible and useless as most of them would not enter the optimal solution. A \textit{restricted master problem} is then defined, where instead of solving the problem on all its solution space, only a subset of variables (or columns) is considered. The column generation method iterates between two different problems: the restricted master problem and a \textit{pricing subproblem}. The pricing subproblem aims at looking for new columns with negative reduced cost to add to the restricted master problem to improve its optimal value. % If the decision variables are integer, these constraints are % relaxed to have continuous decision variables. Each restricted master problem is then solved to optimality and iteratively enriched with % new columns produced by the \textit{pricing subproblem} which uses as input data the optimal dual variables of the % restricted master problem. Even if the set of columns of the original master problem has a huge or exponential size, it is often possible to find the minimum reduced cost column in short or polynomial time, depending on the structure of the pricing problem. The method stops when no more negative reduced cost column is produced by the pricing subproblem, which means that the linear relaxation of the original master problem is solved to optimality. Nevertheless, in the case of integer variables the optimal solution of the linear relaxation could be far away from the optimal solution of the original master problem. In this case it is necessary either to start a branch and price algorithm (see for example \cite{BrAdPr} for a survey on this topic) or to have a complete pricing of the whole set of columns candidate to be in the optimal solution (see for example \cite{ITCVRP} for a reference on this topic). % The column generation method iterates between two different problems: a \textit{restricted master problem} and % a \textit{pricing subproblem}. A \textit{restricted master problem} is the master problem where instead of solving the problem on all % its solution space, only a subset of variables (or columns) is considered. If the decision variables are integer, these constraints are % relaxed to have continuous decision variables. Each restricted master problem is then solved to optimality and iteratively enriched with % new columns produced by the \textit{pricing subproblem} which uses as input data the optimal dual variables of the % restricted master problem. % The \textit{pricing subproblem} aims at looking for new columns with negative reduced cost (associated with the former dual variables), % that if added to the master problem would improve its optimal value (except in case of degeneracy). Even if the set of columns of the % original master problem has a huge or exponential size, it is often possible to find the minimum reduced cost column in pseudo polynomial time. % The method stops when no more negative reduced cost column is produced by the pricing subproblem, which means that the linear relaxation % of the original master problem is solved to optimality. Nevertheless, in the case of integer variables the optimal solution of the % linear relaxation could be far away from the optimal solution of the original master problem. In this case it is necessary either to start % a branch and price algorithm (see for example \cite{BrAdPr} for a survey on this topic) or to have a complete pricing of the whole set of columns % candidate to be in the optimal solution (see for example \cite{ITCVRP} for a reference on this topic). Our approach is a heuristic method since only a subset of the whole set of possible columns are implicitly considered, due to the time limit imposed by the competition. This method is indeed heuristic for several reasons. The first reason is that production levels are discretized over time periods. This implies that the pricing subproblem is approximately solved using a shortest path algorithm in a graph with a finite number of nodes associated with production levels, which is much faster than solving it with continuous variables. The second reason is that scenarios are aggregated into a single mean scenario. Finally, the column generation method leads to find the optimal solution of the relaxed problem on the subset of columns that could be produced with respect to the discretization considered. A branch and price algorithm was not envisaged due to the time limits. At the end of the column generation process, once the outages have been scheduled, two linear programming problems are solved to optimality for deciding the quantity of fuel and the optimal production planning over all scenarios for this set of outages. %The discriminant option chosen for our approach is to spend time to generate a high-quality schedule of outages, as guided by the cost objective of column generation. The main contributions of this paper are the following: \begin{itemize} \item an original set partitioning-like reformulation of the problem; %associated to a Benders decomposition%, and its decompostion into two subproblems \item a pricing subproblem solved by shortest path computations on a particular state graph; \item a diversification method that changes the objective function, adding columns that may have a positive reduced cost; \item two linear programming problems solved to optimality for deciding the quantity of fuel and the production planning for the selected set of outages. \end{itemize} The discriminant strategy chosen for our approach is to spend time to generate a single high-quality schedule of outages, guided by the cost objective of column generation, rather than enumerating several schedules of outages in a shorter time. This paper is organized as follow. The electricity production planning problem is presented in Section \ref{sec:problemstatement}. Section \ref{sec:model} introduces the mathematical model that formulates the problem. Section \ref{sec:decomposition} explains the decomposition applied and describes the overall column generation method used to solve it. Section \ref{sec:pricing} describes the formulation of the pricing subproblem and its resolution method. Section \ref{sec:solving_the_master_problem} gives the model of the final master problem that is solved for fixing dates of outages and critical stock-level periods, ending the first stage of the method. Section \ref{sec:diversification} shows how to improve the pricing subproblems in order to generate a greater number of columns with higher diversity. Section \ref{sec:ReloadProduction} presents the second-stage process used for fixing reload and production quantities for each plant. Section \ref{sec:ComputationalResults} gathers the numerical results obtained on the test instances while section \ref{sec:conclusion} concludes the paper. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Problem statement} \label{sec:problemstatement} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The notations used in this paper mix the notations used in the ROADEF/EURO 2010 Challenge document \cite{sujetOnline} and notations issued from the paper by Lubsy {\em et al.} \cite{roadef10Danois}. The various indexes used are scenarios $s\in S$, time steps $t\in \timesteps = \{1,\ldots,\tmax\}$, weeks $w \in \weeks = \{1,\ldots,\wmax\}$, non-nuclear plants $j\in J$, nuclear plants $i\in I$ and production cycles $k\in K_i=\{1,\ldots,|K_i|\}$ for each nuclear power plant $i$, $|K_i|$ being the maximum number of cycles that can start during the considered period. The decision variables introduced in the Challenge document are: \begin{itemize} \item $ha_{ik}$ : week of cycle $k$ during which nuclear plant $i$ goes offline (i.e. week of decoupling) \item $p_{jts}$ : production of non-nuclear plant $j$ during time step $t$ over scenario $s$ \item $p_{its}$ : production of nuclear plant $i$ during time step $t$ over scenario $s$ \item $r_{ik}$ : reload performed during the outage of cycle $k$ of nuclear plant $i$ \item $x_{its}$ : stock of fuel of nuclear power plant $i$ at the beginning of time step $t$ over scenario $s$ \item $ec_{ik}$ : set of time steps composing the production campaign of cycle $k$ of nuclear plant $i$ \item $ea_{ik}$ : set of weeks composing the outage of cycle $k$ of plant $i$ \end{itemize} Variables $ha$ are integer, all other variables but $ec$ and $ea$ are continuous variables. The objective function is: \[ \min \sum_{i \in I} \sum_{k \in K} C_{ik} r_{ik} + \frac{1}{S} \sum_{s \in S} \left( \sum_{t \in T} \sum_{j \in J} C_{jts} p_{jts} D^t -\sum_{i \in I} C^f_i x^f_{is}\right) \] \noindent where $C_{ik}$ is the proportional cost of fuel during cycle $k$ for nuclear plant $i$, $C_{jts}$ is the proportional production cost for time $t$ and scenario $s$ for non-nuclear plant $j$, $C^f_i$ is the selling price (or the gain) for the remaining fuel $x^f_{is}$ and $D^t$ is the length of time step $t$. The set of constraints is the following, taking the same numbering as in the Challenge document. Some notations used throughout the paper are introduced with these constraints. %The global constraints that have been used for writing the set partitioning formulation of the problem are the following: \begin{description} \item[CT1] Demand satisfaction: \[ \sum_{j \in J} p_{jts} + \sum_{i \in I} p_{its} = DEM_{t,s} \quad \forall s,t \] \end{description} % as they will be put into the master problem only at the end of the column generation process. %\subsection{Local constraints holding for each plant} %The local constraints that will belong to the pricing subproblem associated with each {\it nuclear} plant in the column generation scheme are the following: \begin{description} \item[CT2] For non-nuclear plants, we only have minimum and maximum production constraints: \[ \underline{P}_{jts} \leq p_{jts} \leq \overline{P}_{jts} \quad \forall t,s,j \] \noindent where $\underline{P}_{jts}$ and $\overline{P}_{jts}$ are the minimum and maximum production for non-nuclear plants $j$ on time step $t$ for scenario $s$. \end{description} \begin{description} \item[CT3] Offline power: if the nuclear power plant is offline the production must be 0 \[ t \in ea_{ik} \Rightarrow p_{its} = 0 \quad \forall t,s,j \] \end{description} \begin{description} \item[CT4] Minimum power: if the nuclear power plant is online the production must be greater or equal than 0 \[ p_{its} \geq 0 \quad \forall i,t,s \] \end{description} \begin{description} \item[CT5] Maximum power before activation of imposition of power profile constraint \[ \text{if} \quad (t \in ec_{ik}) \text{ and } x_{its}\geq BO_{ik} \Rightarrow p_{its} \leq \overline{P}_{it} \quad \forall i,t,s,k \] \noindent where $BO_{ik}$ is a bound on stock of fuel that activates the imposition of power profile during cycle $k$, and $\overline{P}_{it}$ is the maximum power of plant $i$ during time step $t$ \end{description} \begin{description} \item[CT6] Maximum power after activation of imposition of power level constraint : \begin{eqnarray*} \left.\begin{array}[c]{llr} \text{ if } & t\in ec_{ik} \text{ and } x_{its}0$ is a very big parameter so that no bound holds if $Y_{ik}=0$, similarly to \cite{roadef10Danois}. The computation of the final level of fuel at the end of the time horizon is made in constraint (\ref{stock_final}). In these constraints, only one term in the right-hand-side sum will be non-zero by (\ref{eq:y_kk}) and will correspond to the final stock.\\ Constraints (\ref{x_ikts}) define the level of stock during a production time step as the stock level after reload minus the total production of the campaign up to this step. The modulation constraints [CT12] are expressed in (\ref{modulation}). These constraints use variables $z_{ikts}$ that are forced to be 1 as soon as $x_{ikts}\geq BO_{ik}$ by constraints (\ref{z_{ikts}}).\\ Constraints (\ref{BOprofil}) [CT6] impose the piecewise linear decreasing production profile as soon as $PB_{ik}(x_{ikts})\underline{P}_{it}D^t \leq x_{ikts} < BO_{ik}$, which implies $\xi_{ikts}=1$ via constraints (\ref{xi_ikts_1}-\ref{xi_ikts_2}), or no production otherwise. %The definition of $\bar{x}^d_{iks}$ takes into account the fact that when a power plant goes for reload the remaining fuel is lost. Constraints (\ref{eq:limit_amax}), (\ref{eq:limit_smax}) and (\ref{eq:final_stock}) bound the quantity of fuel that can be reloaded in a cycle. %Constraints (\ref{eq:less_stock}) represent the phenomenon that the rest of fuel of the previous cycle $k-1$ reduces before the refueling of cycle $k$. Constraints (\ref{eq:BO_level_uno}) and (\ref{eq:BO_level_due}) state that when the level of fuel is below the value $\mathtt{BO}$ the variable $z_{itks}$ allow the production to follow the piecewise-linear profile. %Constraints (\ref{eq:activer_xi_uno}) and (\ref{eq:activer_xi_due}) impose that the production must follow the piecewise-linear profile if there are enough quantity of stock $x_{itks}$. Finally constraints (\ref{eq:prodImposed_BO}) impose that the production profile must be equal to the imposed $PB_{ik}(x_{itks})\overline{P}_{it}$ value. Constraints (\ref{eq:relationxt_x}) describe the relation between the two values of stock (i.e., $x_{itks}$ and $x^d_{iks}$). The complexity of the model is mainly due to the combinatorial choice for binary variables $y_{i\week k}$, $z_{itks}$ and $\xi_{itks}$, and the nonlinearity of constraints (\ref{x1x2}), (\ref{stock_final}), (\ref{eq:activer_xi_uno}), (\ref{BOprofil}) and (\ref{xi_ikts_2}). %As we said before the first one states that once the fuel stock level in a given nuclear power plant falls below a certain threshold production must follow a piecewise linear decreasing profile, while the second tries to ensure a high utilization of the nuclear power plants by stipulating that the average deviation of the production cannot be more than a certain tolerance from the maximum possible production level (however only prior to the aforementioned thresholds). \section{Decomposition and solution method} \label{sec:decomposition} \subsection{Decomposition} The set of constraints can be decomposed in two sets. The {\it global} constraints are demand constraints (\ref{eq:equil_dem}) and incompatibility constraints (\ref{eq:cliques_CR}) which link all power plants, and constraints (\ref{eq:limit_prod_j}) that bound the production in thermal plants. The {\it local} constraints (\ref{eq:limit_prod_i})-(\ref{xi_ikts_2})) can be decomposed by each nuclear power plant $i$. Assume that we are able to find for each nuclear power plant $i$ the set of all possible production plans (i.e., feasible values of $p_{ikts}, \forall t \in \timeline$, outages and reload quantities). Define $\lc{P}$ as the index set of all feasible production plans for the set of nuclear power plants, and let $\lc{P}_i$ be the subset of production plans associated to nuclear power plant $i$. Therefore, we can write a set partitioning-like reformulation (i.e., a set partitioning model with some more constraints) of the original problem where binary variables $y_\ell$ indicates whether production plan $\ell \in \lc{P}_i$ is selected for plant $i$, and variables $p_{jts}$ represent the thermal energy used for each scenario $s \in S$ and each period of time $t \in T$. \begin{align} (F) \quad z(F) = \min &\sum\limits_{i \in I} \sum\limits_{\ell \in \lc{P}_i} c_{\ell} y_{\ell} + \frac{1}{|S|}\sum_{s \in S}\sum_{t \in \timeline}\sum_{j \in J} C_{jts} p_{jts} & \label{eq:F_objF}\\ \text{s.t.}\quad &\sum\limits_{i \in I} \sum\limits_{\ell \in \lc{P}_i} p_{\ell ts} y_{\ell} + \sum_{j \in J} p_{jts} = DEM_{ts}, & \forall t \in \timeline, \forall s \in S \label{eq:F_dem}\\ &\sum\limits_{\ell \in \lc{P}_i} y_{\ell} = 1 &\forall i \in I &\label{eq:F_setpar}\\ & \underline{P}_{jts} \leq p_{jts} \leq \overline{P}_{jts}, & \forall j \in J, \forall t \in \timeline, \forall s \in S & \label{eq:F_limit_prod_j} \\ &\sum\limits_{\ell \in \lc{P}} a_{\ell G}y_\ell \leq \alpha_G, &\forall G \in \mathcal{G} & & \label{eq:F_cliques}\\ &y_{\ell} \in \{0,1\} & \forall \ell \in \lc{P}, &\label{eq:F_binF}\\ &p_{jts} \geq 0 & \forall j \in J, \forall t \in \timeline, \forall s \in S & \label{eq:F_positF} \end{align} \noindent where $p_{\ell ts}$ is the quantity produced in plan $\ell$ in time step $t$ over scenario $s$, and $a_{\ell G}$ are 0-1 coefficients associated with plan $\ell$ in constraint set $G$. Equation (\ref{eq:F_objF}) is the resulting objective function to be minimized. Constraints (\ref{eq:F_dem}) specify that the energy demand of each period $t \in T$ and for each scenario $s \in S$ must be satisfied. Constraints (\ref{eq:F_setpar}) impose that for each nuclear power plant one and only one column (refuel and production plan) must be chosen. Constraints (\ref{eq:F_limit_prod_j}) bound the quantity of energy that a thermal power plant can produce each time period. Constraints (\ref{eq:F_cliques}) limit the number of nuclear power plants that can stop together. Constraints (\ref{eq:F_binF}) and (\ref{eq:F_positF}) define the range of the variables. %As we said above the problem of generating the production planning (i.e., the set of coefficients associated to the variable $y_\ell$) consists in resolving the following problem: one for each nuclear power plant $i$. %Then, suppose that for each power plant $i$ we are able to solve the subproblem defined by an objective function, described below, and by the constraints (\ref{eq:r_minmax})-(\ref{eq:relationxt_x}) which define the value of the refuel, the production quantities and the outages for refueling all over the time period considered (i.e., a complete production planning). Note that constraints (\ref{eq:F_dem}), (\ref{eq:F_limit_prod_j}) and (\ref{eq:F_cliques}) correspond to constraints (\ref{eq:equil_dem}), (\ref{eq:limit_prod_j}) and (\ref{eq:cliques_CR}). %\begin{align} %(P) \quad z(P) = \sum\limits_{i\in I} \sum_{k \in K} C_{ik} r_{ik} - \sum_{i \in I} C_i^f x_{is}^f \\ %\text{s.t. (\ref{eq:r_minmax})-(\ref{eq:relationxt_x})} %\end{align} %%%%% NB LAURENT : a quoi sert ce probleme P? Coupe la lecture en introduisant des notations inutiles + ambiguite avec pricing %%%%%%%%% %This problem can be decomposed in several subproblems: one for each nuclear production plant $i \in I$. Each subproblem is still very difficult to solve due to the non linear constraints and binary variables. Nevertheless, it contains in some way the lot sizing problem with stock balance equations. We exploit this similarity in the solution approach. \subsection{Solution method} \label{subsec:approach} The efficiency of the approach largely depends on the tractability of problem (F) and the $|I|$ subproblems associated with each nuclear power plant. These problems have very large scale given the number $|S|$ of scenarios and the number $T$ of time steps. To deal with this complexity, we chose to simplify the original set of data by aggregating the $S$ scenarios into a single {\it average} scenario $\bar{s}$, and the time steps $t=1,\ldots,T$ into weeks $\week=1,\ldots,W$. We note for every week $\week \in W$: \begin{eqnarray*} DEM_{\week\bar{s}} &=& \frac{1}{|S|} \sum\limits_{s\in S} \ \sum\limits_{t \in \timesteps{\week}} D^t.\ DEM_{ts}\\ C_{j\week\bar{s}} &=& \frac{1}{|S|} \sum\limits_{s\in S} \frac{\wmax}{\tmax} \sum\limits_{t \in \timesteps{\week}} C_{jts}\\ \underline{E}_{j\week\bar{s}} &=& \frac{1}{|S|} \sum\limits_{s\in S}\ \sum_{t \in \timesteps{\week}} D^t.\ \underline{P}_{jts}\\ \overline{E}_{j\week\bar{s}} &=& \frac{1}{|S|} \sum\limits_{s\in S}\ \sum_{t \in \timesteps{\week}} D^t.\ \overline{P}_{jts}\\ \overline{E}_{i\week}&=& \sum_{t\in \timesteps{\week}} \overline{P}_{it} \ .\ D^t \end{eqnarray*} %%%%%% NB LAURENT -> Daniel : j'ai replace tes notations EMAX etc. en section 9 car on n'utilise pas les notations avant longtemps %%%%%%%%%%% From now on, production quantities stand for each week $w \in W$, and we use notations $e_{\ell w}$ for the production of column $\ell \in \lc{P}_i$ in week $w$, and production variables $e_{jw}$ for thermal power plants. Therefore the problem $(F)$ to solve becomes \begin{align} (\bar{F}) \quad z(\bar{F}) = \min &\sum\limits_{i \in I} \sum\limits_{\ell \in \lc{P}_i} c_{\ell} y_{\ell} + \sum_{\week \in \weeks} C_{j\week\bar{s}} e_{\week{}j} & \label{eq:Fw_objF}\\ \text{s.t.}\quad &\sum\limits_{i\in I} \sum\limits_{\ell \in \lc{P}_i} e_{\week\ell}y_{\ell} + \sum_{j \in J} e_{\week{}j} = DEM_{\week\bar{s}} & \forall \week \in \weeks \label{eq:Fw_dem}\\ &\sum\limits_{\ell \in \lc{P}_i} y_{\ell} = 1 &\forall i\in I &\label{eq:Fw_setpar}\\ & \underline{E}_{j\week\bar{s}} \leq e_{\week{}j} \leq \overline{E}_{j\week\bar{s}} & \forall j \in J, \forall \week \in \weeks & \label{eq:Fw_limit_prod_j} \\ & \sum\limits_{i\in I} \sum\limits_{\ell \in \lc{P}_i} a_{\ell G}y_\ell \leq \alpha_G & \forall G \in \mathcal{G} & \label{eq:Fw_cliques}\\ &y_{\ell} \in \{0,1\} & \forall \ell \in \lc{P} &\label{eq:Fw_binF}\\ &e_{\week{}j} \geq 0 & \forall j \in J, \week \in \weeks & \label{eq:Fw_positF} \end{align} Formulation $\bar{F}$ is impractical to solve even for instances of moderate size, since the number of variables is typically exponential, which justifies a column generation approach. We define a linear relaxation $(L\bar{F})$ of $(F)$: %%% NB Laurent : ce n'est pas LA relaxation lineaire de F on on relache aussi les cliques %%%%% \begin{align} (L\bar{F}) \quad z(L\bar{F}) = \min &\sum\limits_{i \in I} \sum\limits_{\ell \in \lc{P}_i} c_{\ell} y_{\ell} + \sum_{w \in W} C_{jw,\bar{s}} e_{wj} & \label{eq:LF_objF}\\ \text{s.t.}\quad &\sum\limits_{i\in I} \sum\limits_{\ell \in \lc{P}_i} e_{w \ell}y_{\ell} + \sum_{j \in J} e_{wj} = DEM^{w,\bar{s}}, & \forall w \in W \label{eq:LF_dem}\\ &\sum\limits_{\ell \in \lc{P}_i} y_{\ell} = 1 &\forall i\in I &\label{eq:LF_setpar}\\ & \underline{E}_j^{w,\bar{s}} \leq e_{wj} \leq \overline{E}_j^{w,\bar{s}}, & \forall j \in J, w \in W & \label{eq:LF_limit_prod_j} \\ %& \sum\limits_{i=1}^{I} \sum\limits_{\ell \in \lc{P}_i} y_\ell a_{i\week} \leq \alpha_G, & \forall G \in \mathcal{G} & \label{eq:Mh_cliques}\\ & 0 \leq y_{\ell} \leq 1 & \forall \ell \in \lc{P}, &\label{eq:LF_binF}\\ &e_{wj} \geq 0 & \forall j \in J, \forall \week \in \weeks & \label{eq:LF_positF} \end{align} When passing from $(L\bar{F})$ to $(F)$, constraints (\ref{eq:Fw_cliques}) have been removed since their dual variables are hard to compute in the pricing problem within the limited computational time imposed by the competition. For compensation we allow to add to problem $L\bar{F}$ some columns with positive reduced costs as described in section (\ref{sec:diversification}). The solution of problem $L\bar{F}$ provides the dual variables $\mu_w$ and $u_i$ associated with constraints (\ref{eq:Fw_dem}) and (\ref{eq:Fw_setpar}), respectively. The pricing problem exploits them for generating the column(s) with minimum reduced cost: \begin{equation} {\ell}^*= \argmin_{{\ell} \in \lc{P}} \bar{c}_{\ell} \quad \text{ with } \quad c_{\ell}=c_{\ell} - \sum_{w \in W} e_{w \ell} \mu_w - u_i\label{subp} \end{equation} \noindent The pricing subproblem ($P_i$) associated with each nuclear power plant $i \in I$ can be written as follows, given the decomposition previously described and aggregation of time steps into weeks. Indexes of various scenarios have disappeared and the $t$ indexes have been replaced by $w$ indexes in data and variables. We denote by $W'_{ik}$ the set of weeks when production is possible in cycle $k$. \begin{align} %(P) \quad \bar{c}_\ell = \min\limits_{r_{i(K)}e_{(\weeks)\ell}x^f_{i}} & \sum_{k \in K} C_{ik} r_{ik} - C_i^f %x_{i}^f - u_i - \sum\limits_\week \mu_\week e_{\week\ell} %& \label{eq:Pobj_CR}\\ (P_i) \quad \min & \sum_{k \in K_i} C_{ik} r_{ik} - C_i^f x^f_{i} - \sum\limits_{w \in W, k \in K_i(w)} p_{ikw} \mu_{w} - u_i \\ \text{s.t.} \ %%%%%%% local constraints%%%%%%%% & 0 \leq p_{ikw} \leq \overline{P}_{iw} \eta(i,w,k), & k \in K_i, w \in W & \label{SPeq:limit_prod_i} \\ & Y_{ik}= \sum\limits_{w \in W_{ik}} y_{iwk} &\forall k \in K_i & \label{SPeq:Yik}\\ & \underline{R}_{ik} Y_{ik} \leq r_{ik} \leq \overline{R}_{ik} Y_{ik} &\forall k \in K_i &\label{SPeq:r_minmax}\\ & Y_{i,k-1} \geq Y_{i,k} & \forall k \in K_i \setminus \{1\} &\label{SPeq:y_kk}\\ & x^1_{i,k+1} = x^2_{ik} - \sum_{w \in W'_{ik}} p_{ikw} &\forall k \in K_i & \label{SPx1x2} \\ %\mathtt{[CT10,CT9]} & x^2_{ik} = \left[\left( \frac{Q_{i,k} -1}{Q_{i,k}}\right)( x^1_{ik} - \mathtt{BO}_{i,k-1}) + r_{ik} + \mathtt{BO}_{i,k} \right]Y_{ik} + x^1_{ik}(1-Y_{ik}) &\forall k \in K_i & \label{SPx1x2} \\ & x^1_{ik} \leq A_{ik} + (1 - Y_{i,k})(M - A_{ik}) &\forall k \in K_i & \label{SPAmax} \\ & x^2_{ik} \leq S_{ik} + (1 - Y_{i,k})(M - S_{ik}) &\forall k \in K_i & \label{SPSmax} \\ & x^f_{i} = \sum_{k \in K_i} x_{i,k+1,s}(Y_{ik} - Y_{i,k+1}) & & \label{SPstock_final} \\ & x_{ikw} = x^2_{ik}-\sum\limits_{w'\leq w} p_{ikw'} & \forall k \in K_i, w \in W'_{ik}, s \in S & \label{SPx_ikts} \\ %%%%%%% ajout Laurent contraintes Roberto pour [CT12] et [CT6]%%%%% & \sum\limits_{w \in W'_{ik}} (\underline{P}_{iw}-p_{ikw}) z_{ikw} \leq \underline{M}_{ik} & \forall k \in K_i & \label{SPmodulation} \\ & M \ z_{ikws} \geq \frac{x_{ikws}}{BO_{ik}}-1 & \forall i \in I, k \in K_i, t \in T_{ik}, s \in S & \label{SPz_ikts} \\ & (1-\epsilon)\xi_{ikw} \leq \frac{p_{ikw}}{PB_{ik}(x_{ikw})\underline{P}_{iw}} \leq (1+\epsilon) \xi_{ikw} & \forall k \in K_i, w \in W'_{ik} & \label{SPBOprofil} \\ & \xi_{ikw} \geq 1-\frac{x_{ikw}}{BO_{ik}} & \forall k \in K_i, w \in W'_{ik} & \label{SPxi_ikts_1} \\ & \xi_{ikw} \geq 1-\frac{PB_{ik}(x_{ikw})\underline{P}_{iw}}{x_{ikw}} & \forall k \in K_i, w \in W'_{ik} & \label{SPxi_ikts_2} \\ & y_{i\week k}, Y_{ik}, z_{ikw}, \xi_{ikw} \in \{0,1\} & & \label{bin_var} \\ & p_{ikw}, p_{jw}, x^1_{ik}, x^2_{ik}, x^f_{is}, x_{ikw} \geq 0 & & \label{cont_var} \end{align} The Pricing problem (P) consists in optimizing every problem $(P_i)$ for $i \in I$ and output the column with best optimal value (minimum reduced cost). The columns with negative reduced cost output by $(P)$ are added to problem $L\bar{F}$ and the process iterates between $L\bar{F}$ and $P$ until $P$ finds no more negative reduced cost column. As we mentioned in section 1, the proposed column generation method is a heuristic scheme. Indeed, due to the time constraints imposed by the competition, only an approximation of the pricing problem $P$ is solved at each iteration. This holds because the pricing subproblem is a shortest path problem solved on a graph whose paths do not explicitly represent all feasible plans, as explained in section 5. At the end of the column generation, the set of columns forms a relaxed restricted master program which is called the Restricted $L\bar{F}$ ($RL\bar{F}$), and solved to optimality as explained in section \ref{sec:solving_the_master_problem}. The obtained optimal solution is a set of columns and variables $p_{wj}$ satisfying constraints (\ref{eq:LF_dem})-(\ref{eq:LF_limit_prod_j}) at minimum cost. Each column is a production plan containing the dates for the outages and the periods when stock reaches a critical level equal to $BO_{ik}$ (if such periods exist for a given nuclear plant). Therefore the solution obtained at the end of the first phase is used for fixing the decision variables on the outages and the date when the $BO$-threshold is attained. The advantage is to avoid to deal explicitly with non linear constraints, since once one has fixed $BO$, the corresponding decreasing profile is completely defined for the successive periods until the next outage. The only decisions that remain to be made in the second phase are the reload variables $r_{ik}$ and the production variables $p_{ikts}$ , which are easier to settle once the binary outage variables $y_{i\week k}$ have been fixed. For setting the values of $p_{ikts}$, $p_{jts}$ and $r_{ik}$, it is necessary to use the whole model with all time steps $t =1,\ldots,T$ and all scenarios $s \in S$. Since it has no more binary variables, it can be solved efficiently by Linear Programming. %The determination of outage variables $y_{iwt}$ is obtained by solving with column generation (CG) the set partitioning reformulation of the problem. %The reload and production variables are fixed for each plant in a second stage. The two-stage method is stated below. \begin{description} \item[\textbf{Phase 1}] Fixing dates for outages and BO stock-level periods \begin{description} \item[\textit{subphase 1}] Solve the aggregate set partitioning-like formulation $L\bar{F}$ by column generation %\end{description} \begin{itemize} \item Aggregate data, passing from $|\timeline \times S|=(5000 \times 500)$ to $|\weeks \times \{\bar{s}\}|=(260 \times 1)$ \item Create an acyclic graph $G_i$, for each nuclear power plant $i$ \begin{itemize} \item Repeat until a time limit is elapsed or no new column of negative reduced cost exists \begin{itemize} \item Solve a shortest path subproblem in $G_i$, for each nuclear power plant $i$ \item Add the new columns (corresponding to shortest paths) to the $L\bar{F}$ master problem and solve it by the Simplex method \end{itemize} \end{itemize} \end{itemize} %\begin{description} \item[\textit{subphase 2}] Give the final restricted master problem to a MIP solver adding all the clique-like inequalities (\ref{eq:Fw_cliques}), exploiting the fact that only a small number of nuclear power plants can be refueled at the same time. Outputs are the dates for outages and $BO$ stock-level periods. \end{description} %\end{itemize} \item[\textbf{Phase 2}] Solve LP problems for defining the reload level for each date of outages and optimizing the production planning for each power plant and each scenario. \end{description} Phase 2 is described in Section \ref{sec:ReloadProduction}. Sections \ref{sec:pricing} and \ref{sec:solving_the_master_problem} detail the first phase of the method. In phase 1, for every nuclear plant $i \in I$, two decisions are fixed: dates for outages respecting the scheduling constraints on outages (i.e., \ref{eq:Fw_cliques}) and periods when stock reaches the $BO_{ik}$ level, with imposed non linear decreasing power profile on the following periods. This step is the most time-consuming step of the method. %It is divided in two stages. In the first stage, the initial problem is transformed into an {\it aggregate} covering master problem where time steps $t$ are aggregated into weeks $h$ and a single {\it average} scenario called $\bar{s}$ is considered, so that the demand to cover is the average demand over all $S$ scenarios. \section{Solving the pricing problem} \label{sec:pricing} The pricing problem ($P_i$) consists in finding a reload and production plan for each nuclear power plant $i$. The problem is close to a lot sizing problem and can be solved by a shortest path algorithm in a graph. The main point is how to build the graph. %(see for example \cite{IntProg}, \cite{DynLotSize} and \cite{EcoLotSize}). The difficulty of this step is that the decision variables are a combination of continuous and binary variables. The binary variables define the possible outage dates while the continuous variables define the quantity of fuel that must be reloaded and the production for each time step. The further difficulty is due to the $BO$ piecewise non linear function. Last but not least comes from the limited computational time imposed by the competition. Therefore, in the proposed approach, we chose to design a graph in which a production plan would be made of a relatively small number of arcs, each production campaign being represented by at most two arcs, and in such a way that cost computations of each arc could be managed efficiently. For efficiency reasons, we made the simplification of discretizing the possible fuel stock levels at the beginning of outages, so that the stock level could be fixed once for all at both ends of each arc of the graph. In the construction of the graph constraints (\ref{stock_level}) (\ref{eq:activer_xi_due}) are taken into account either when a node is defined or either when an arc is defined. The pricing problem has been implemented as a three-step process: \begin{algorithm}[H] \caption{The pricing algorithm} \label{algo:pricing_problem} \begin{algorithmic} % \STATE set all $e(h, \ell)$ to $EMAX_h$ % \STATE order the set of weeks by decreasing values of % $\left(C_{i,k_n}-\mu_{\week{}_m}\right)$~: $h_1, h_2 \dots $ \FORALL{nuclear power plants $i$} \STATE \textit{\textbf{1)}} create an acyclic static \textit{stock graph} $G_i=(N_i,A_i)_i$; \WHILE{there is a column to add and the time is not elapsed} \STATE add the column to the master and calculate the new reduced cost \STATE \textit{\textbf{2)}} for each arc $a \in A_i$ compute weight $w_a$ (optimal reduced cost on the arc) ; \STATE \textit{\textbf{3)}} generate a new column by finding in $G_i$ a path of minimum weight; \ENDWHILE \ENDFOR % \WHILE{$\left(C_{i,k_n}-\mu_{\week{}_m}\right)>0$ \AND constraints \ref{eq:wadue},\ref{eq:watre},\ref{eq:waquattro},\ref{eq:wacinque} can be satisfied} % \STATE nullify $e_{\week{}_m \ell}$ % \STATE increment $m$ % \ENDWHILE \end{algorithmic} \end{algorithm} Step 1 is quite time consuming %(5 to 30 minutes for all power plants), j'ai enlevé car on est dans la methode generique, pas experiments but occurs only once at the beginning of the whole process. Steps 2 and 3 are light processes %(1 to 10 seconds) idem and can be serialized with the computation of new dual variables $\mu_\week$. %\begin{figure}[H] % \centering % \scalebox{1}{ % \input{column_process}} % \caption{Column generation process} % \label{fig:column_process} %\end{figure} %\bigskip %REMARQUE LAURENT SUR DESSIN: j'ajouterais dans le dessin sur la patate en haut Ă  droite, $\bar{c}_{{\ell}^*}<0$ sur la flèche qui part vers le bas, et ajouter une nouvelle %flèche de cette mĂŞme patate vers une patate "STOP" avec la condition $\bar{c}_{{\ell}^*} \geq 0$ car pas de critère d'arrĂŞt dans le schĂ©ma \subsection{Creation of a static graph}\label{ssec:graph} \paragraph{Description of the graph.} Each node $n \in N_i$ in graph $G_i$ represents a possible state of power plant $i$ and is characterized by a 4-tuple~$(\mathit{type}_n, x_n, \week_n, k_n)$, where: \begin{itemize} \item $\week_n$ is the week associated with $n$~; \item $x_n$ is the fuel stock of power plant $i$ \textit{at the beginning} of week $\week_n$~; \item $k_n$ is the id of the campaign in which the power plant will be during week $\week_n$~; \item $\mathit{type}_n$ is the node type. There are two main types of nodes: \begin{itemize} \item an ``\texttt{OUTAGE}'' node $n$ marks the start of a new production campaign and $\week_n$ is the week during which the outage occurs~; \item a ``\texttt{BO}'' node $n$ marks the start of a time period during which the fuel stock will be lower than $BO_{ik_n}$ and the power level will be imposed by constraint \texttt{[CT6]}~; the stock $x_n$ of a \texttt{BO} node is indeed a little higher than $BO_{ik_n}$, but ensures that $BO_{ik_n}$ can be reached \textit{during} week $\week_n$ by imposing the maximal power production of constraint \texttt{[CT5]} at the beginning of the week. In the model, it corresponds to setting the constraints (\ref{eq:PBO_level_uno})-(\ref{eq:Prelation_pitk_ewl}). \end{itemize} Besides, two nodes with special types are added, corresponding to the \texttt{initial} and \texttt{final} states of a power plant. \end{itemize} Each arc $a=(n_1,n_2) \in A_i$ in graph $G_i$ represents a possible direct transition from state $n_1$ to $n_2$. \begin{itemize} \item If $n_1$ is an \texttt{OUTAGE} node, this means that it is possible to reach state $n_2$ by proceeding to an outage on week $\week_{n_1}$ followed by a production campaign, without need for any other outage and without reaching the $BO$ stock level during the time interval $[\week_{n_1},\week_{n_2}]$. Note that $n_2$ can be either a \texttt{BO} node or an \texttt{OUTAGE} node. If it is a \texttt{BO} node, then $k_{n_2}=k_{n_1}$; else $k_{n_2}=k_{n_1}+1$. Arcs leaving from \texttt{OUTAGE} nodes will be referred to as \texttt{OUTAGE} arcs. \item If $n_1$ is a \texttt{BO} node, this means there exists a production plan leading from $n_1$ to $n_2$ that respects the imposed profile of constraint \texttt{[CT6]} with the acceptable $\epsilon$ tolerance. Note that here $n_2$ has to be an \texttt{OUTAGE} node and thus $k_{n_2}=k_{n_1}+1$. Arcs leaving from \texttt{BO} nodes will be referred to as \texttt{BO} arcs. \end{itemize} Furthermore, arcs are constructed so as to guarantee that technical constraints \texttt{[CT3]} to \texttt{[CT13bis]} (eq.~(\ref{eq:Plimit_prod_i}) to (\ref{eq:Pmodulation})) can be respected for any path from the \texttt{initial} to the \texttt{final} node of $G_i$. \paragraph{Building the graph.} We build graph $G_i$ in the following way: \begin{algorithm}[H] \caption{Creation of graph $G_i$} \label{algo:noeuds} \begin{algorithmic}[1] \STATE create the \texttt{initial} node at week $1$ \FOR{week $\week=1$ \TO $W$} \FORALL{previously created nodes $n$ with $\week_n=\week$} \STATE \label{algostep:noeuds} create new nodes $n'$ with $\week_{n'}>\week$ so that there exists a production plan leading the power plant from state $n$ to state $n'$ without violating constraints \texttt{[CT3]} to \texttt{[CT13bis]} \ENDFOR \ENDFOR \end{algorithmic} \end{algorithm} In order to find new nodes (step \ref{algostep:noeuds} of Algorithm \ref{algo:noeuds}), we follow the beam of possible stock interval(s) over time from a node $n$, as stated in Algorithm \ref{algo:beam}: \begin{algorithm}[H] \caption{Creation of new nodes in $G_i$ from $n=(\mathit{type}_n, x_n, \week_n, k_n)$} \label{algo:beam} \begin{algorithmic}[1] \FOR{week $\week=\week_n+1$ \TO $W-1$} \STATE \label{algostep:possibleStock} compute the possible stock interval $[x_{min_{n,\week}},x_{max_{n,\week}}]$ of power plant $i$ at the beginning of week $\week$ \IF{the outage of campaign $k_n+1$ is allowed at week $\week$ according to constraint \texttt{CT13}} \STATE compute the intersection of $[x_{min_{n,\week}},x_{max_{n,\week}}]$ with $[0,A_{ik_n+1}]$ \STATE \label{algostep:nouvelOutage} create new nodes $n'=(\mathtt{OUTAGE}, x_{n'}, \week, k_n+1)$ with $x_{n'}$ in this intersection (constraint \texttt{CT11}) \ENDIF \IF{$\mathtt{BO}_{ik_n} \in [x_{min_{n,\week}},x_{max_{n,\week}}]$} \STATE \label{algostep:nouveauBO} create new nodes $n'=(\mathtt{BO}, x_{n'}, \week-1, k_n)$ with $x_{n'} \in [x_{min_{n,\week-1}},x_{max_{n,\week-1}}]$ computed so that if the power plant follows the imposed profile of \texttt{CT6} during week $\week-1$, its stock will fall under $\mathtt{BO}_{ik_n}$ at the beginning of week $\week$ \ENDIF \ENDFOR \end{algorithmic} \end{algorithm} For steps \ref{algostep:nouvelOutage} and \ref{algostep:nouveauBO} of algorithm \ref{algo:beam}, we discretize the interval $[x_{\min_{n,\week}}, x_{\max_{n,\week}}]$ of possible stocks of fuel at week $\week$ by a fixed amount of fuel $\Delta_i$. New nodes are only created when needed, so as to ensure that there will be a node in every non-empty intervals $[x_{\min_{n,\week}}, x_{\max_{n,\week}}] \cap [x,x+\Delta_i]$, for all $x \in \mathbb{R}^+$. We compute $\Delta_i$ relatively to the mean capacity of powerplant $i$ after choosing \textit{a priori} a global discretisation parameter $\rho$~: \[ \Delta_i = \rho \ .\ \frac{\sum_{\week \in \weeks} \overline{E}_{i\week}}{\wmax} \] Step \ref{algostep:possibleStock} of algorithm \ref{algo:beam} depends on the node type. \begin{itemize} \item If $type_n=\mathtt{BO}$, we follow the imposed profile of constraint \texttt{[CT6}, i.e. (\ref{eq:PprodImposed_BO})\texttt{]} with the authorized $\epsilon$-deviation, counted positively for $x_{min_{n,\week}}$ and negatively for $x_{max_{n,\week}}$. \item If $type_n=\mathtt{OUTAGE}$, we use \texttt{[CT10}, i.e. (\ref{eq:Pless_stock})\texttt{]} to compute stock variation during the outage, \texttt{[CT7},(\ref{eq:Pr_minmax})\texttt{]} (bounds on refueling) and \texttt{[CT11},(\ref{eq:Plimit_amax})-(\ref{eq:Plimit_smax})\texttt{]} (maximum stock after refueling) to compute the interval after refueling~; when the production campaign starts, $x_{max_{n,\week}}$ is computed using the maximum allowed modulation of \texttt{[CT12},(\ref{eq:Pmodulation})\texttt{]} -- and $x_{min_{n,\week}}$ using no modulation at all. \end{itemize} The node creation process is illustrated by figure \ref{fig:beam} for the second case (starting from an \texttt{OUTAGE} node). \begin{figure}[h] \centering \includegraphics[width=\linewidth]{graph.eps} \caption{Creating new nodes from an \texttt{OUTAGE} node $n$ of campaign $k$ in graph $G_i$ inside the beam of possible stock intervals. Stars inside the beam correspond to the creation of new \texttt{OUTAGE} nodes associated with campaign $k+1$, whereas circles are new \texttt{BO} nodes associated with campaign $k$. Note that no \texttt{BO} node is created at week $T^a_{ik+1}$, because the paths that would go through such a node would violate constraint \texttt{CT13bis}.} \label{fig:beam} \end{figure} When all nodes have been created by the process of Algorithm \ref{algo:noeuds}, a second pass takes place in order to create the arcs of the graph between all nodes that can be connected. \subsection{Weighting the arcs of a graph} The goal of this graph is to compute a production plan $\ell \in \lc{P}_i$ for power plant $i$ with minimal reduced cost $\bar{c}_{\ell}$, which amounts to find according to equation (\ref{subp})~: % Nevertheless, in order to point out that they are variables in this pricing subproblem, $e_{\week{}\ell}$ becomes $e(h,l)$ and $\bar{c}_{\ell}$ becomes $\bar{c}_{\ell} (r,e,x)$ : \[ \min\limits_{r,e,x_i^f} \sum_{k \in K} C_{ik} r_{ik} - C^f_i x^f_i - \sum_{\week \in W} \mu_\week e_{\week} - u_i \] \noindent where \begin{itemize} \item $x^f_i$ is the quantity of remaining fuel at the end of the timeline in the aggregated scenario~; \item $e_{\week}$ is the energy produced at week $\week$ in production plan $\ell$~; \item $\mu_\week$ is the dual variable associated to constraint (\ref{eq:Fw_dem})~; \item $u_i$ is the dual variable associated to constraint (\ref{eq:Fw_setpar}) \end{itemize} Let us recall that a production plan is supported by some unique path ${\ell} \subset A_i$ in graph $G_i$, and that each campaign $k$ in a production plan is carried either by a single \texttt{OUTAGE} arc, or by an \texttt{OUTAGE} arc followed by a \texttt{BO} arc. It comes that the previous expression can be rewritten as: \[ \min \bar{c}_{\ell} = \sum_{a=(n,n') \in {\ell}} \left( C_{ik_n} r_a - C^f_i x^f_a - \sum_{\week \in \weeksof{a}} \mu_\week e_{\week} \right) - u_i \] \noindent where: \begin{itemize} \item $r_a$ is the refuel quantity at the beginning of arc $a$, i.e.~: \[ r_a = \left\lbrace \begin{array}[c]{ll} r_{ik_n} & \text{ if $a$ is an \texttt{OUTAGE} arc}\\ 0 & \text{ if $a$ is a \texttt{BO} arc}\\ \end{array} \right. \] \item $x^f_a$ is the remaining fuel at the end of the timeline if $a$ is the final arc in path $\ell$, i.e.~: \[ x^f_a=\left\lbrace \begin{array}[c]{ll} x^f_i & \text{ if $type_{n'}=$ \texttt{FINAL}}\\ 0 & \text{ otherwise} \end{array} \right. \] \item $\weeksof{a} = [\week_{n},\week_{n'}[$ is the set of weeks spanned by arc $a$. \end{itemize} With this rewriting, we see that each variable that plays a part in the reduced cost of path $\ell$ is \textit{uniquely linked} with a variable associated with some arc of ${\ell}$. Let us call the \textit{reduced cost} of arc $a\in A_i$ the value $c_a$ defined by: \begin{equation} \bar{c}_a=\min\limits_{r,e,x_a^f} C_{ik_n} r_a - C^f_i x^f_a - \sum_{\week \in \weeksof{a}} \mu_\week e_{\week} \end{equation} \noindent While ensuring that constraints of ($P$) hold true for every path in $\ell$, we can easily find the production plan with minimum reduced cost by solving a shortest path in the graph. The reduced cost of this production plan is obtained by substracting $u_i$ from the minimum weight path value. \paragraph{Determining the weight of a \texttt{BO} arc} is straightforward because there is no refuel at the beginning of such an arc ($r_a=0$) and the productions $e_{(\weeksof{a})}$ are fixed for each timestep by constraint \texttt{[CT6}, i.e. (\ref{eq:PprodImposed_BO})\texttt{]}, if we disregard the allowed deviation $\epsilon$. If the arc is \texttt{FINAL}, we can compute $x^f_a$ by substracting the productions from its initial fuel stock. The weight of a \texttt{BO} arc $a$ is given by the following expressions, where everything is known: \begin{equation} \weight_a = \left\lbrace \begin{array}[c]{ll} - \ C^f_i\ .\ \left(x_n - \sum_{\week \in \weeksof{a}} e_{\week{}} \right) - \sum_{\week \in \weeksof{a}} \mu_\week\ .\ e_{\week{}} & \text{ if $a$ is \texttt{FINAL} }\\ - \sum_{\week \in \weeksof{a}} \mu_\week\ .\ e_{\week{}} & \text{ otherwise} \end{array} \right. \label{eq:weightBOarc} \end{equation} \paragraph{Determining the weight of a non final \texttt{OUTAGE} arc} is more intricate because the allowed modulation gives some latitude to productions $e_{(\weeksof{a})}$. However, the campaign $k_n$, the production weeks $\weeksof{a}=[\week_n,\week_{n'}[$ and the levels of fuel $(x_n,x_{n'})$ at the beginning and at the end of the arc are known. Following equation (\ref{eq:Pless_stock}), let $\bar{x}_n=\left( \frac{Q_{ik_n} -1}{Q_{ik_n}} (x_n - \mathtt{BO}_{ik_n}) + \mathtt{BO}_{ik_n-1}\right)$. The problem can then be stated as follows: % the quantity appearing in constraint \texttt{CT10}, in: %\\ %and $\mathit{PMAX}_h = %\sum_{t\in \mathit{timesteps(h)}}\ \mathit{PMAX}_i^t\ .\ D^t$ the %maximal fuel spent during week $\week$ in: \begin{align} %\label{eq:weightOUTAGEarc} (\pla) \quad \weight_a = & \min_{r_a,e_{(\weeksof{a})}} C_{ik_n} r_a - \sum_{\week \in \weeksof{a}} \mu_\week e_{\week{}} &&&\label{eq:waobj}\\ \text{s.t.}\quad & \bar{x}_n + r_a = x_{n'} + \sum_{\week \in \weeksof{a}} e_{\week{}} \quad && \mathtt{[from (\ref{eq:PlotSizing})]} &\label{eq:WAlotSizing}\\ & 0 \leq e_{\week{}} \leq \overline{E}_{i\week} \quad & \forall \week \in \weeksof{a} \quad & \mathtt{[from (\ref{eq:Plimit_prod_i})]} &\label{eq:WAlimit_prod_i}\\ & \underline{R}_{ik_n} \leq r_a \leq \overline{R}_{ik_n} \quad && \mathtt{[from (\ref{eq:Pr_minmax})]} &\label{eq:WAr_minmax}\\ & \bar{x}_n + r_a \leq \overline{S}_{ik_n} \quad && \mathtt{[from (\ref{eq:Plimit_smax})]} &\label{eq:WAlimit_smax}\\ & \sum_{\week \in \weeksof{a}} (\overline{E}_{i\week}-e_{\week{}}) \leq \overline{M}_{ik_n} \quad && \mathtt{[from (\ref{eq:Pmodulation})]} &\label{eq:WAmodulation} \end{align} We can use (\ref{eq:WAlotSizing}) to rewrite $r_a$ and remove it from the problem. It comes that we are looking for the minimum under constraints (\ref{eq:WAlimit_prod_i}) to (\ref{eq:WAmodulation}) of: \[ \weight_a = C_{ik_n}\left(x_{n'}-\bar{x}_n\right) + \sum_{\week \in \weeksof{a}} \left(C_{ik_n}-\mu_\week\right) e_{\week{}} \] As the first part of the sum is \textit{known} for each arc, we only have to find the productions $e_{\week{}}$ that minimize $\sum_{\week \in \weeksof{a}} \left(C_{ik_n}-\mu_\week\right) e_{\week{}}$. We tried two methods to achieve this goal: \paragraph{The exact solution} is obviously given by the Simplex method, since (\pla) is a linear program, or by Algorithm \ref{algo:min_arc_weight_exact} below, where \textit{nullify} has to be understood as ``set as close to zero as possible without violating constraints (\ref{eq:WAlimit_prod_i}) to (\ref{eq:WAmodulation})'': \begin{algorithm}[H] \caption{Exact algorithm for minimizing arc weight} \label{algo:min_arc_weight_exact} \begin{algorithmic} \STATE set all $e_{\week{}}$ to $\overline{E}_{i\week}$ \STATE order the set of weeks by decreasing values of $\left(C_{ik_n}-\mu_{\week_m}\right)$~: $\week_1, \week_2 \dots $ \STATE start with $m=1$ \WHILE{constraints (\ref{eq:WAlimit_prod_i}) to (\ref{eq:WAmodulation}) are not satisfied} \STATE nullify $e_{\week_m}$ \STATE increment $m$ \ENDWHILE \WHILE{$\left(C_{ik_n}-\mu_{\week_m}\right)>0$ \AND constraints (\ref{eq:WAlimit_prod_i}) to (\ref{eq:WAmodulation}) can be satisfied} \STATE nullify $e_{\week_m}$ \STATE increment $m$ \ENDWHILE \end{algorithmic} \end{algorithm} However, we encountered two problems with the exact method: \begin{itemize} \item the overall column generation process was converging too fast, and no new column was generated after about two or three iterations~; this can be overcome however by the \textit{diversification} strategy exposed in section \ref{sec:diversification}~; \item this method was 3 to 10 times slower than the heuristic exposed below (depending on the instances), because it has to walk through the ordered set of weeks for each arc. Thus, on some instances, we were not able to produce enough columns in the imposed time and ended with no valid solution when using this exact algorithm. \end{itemize} \paragraph{We used a heuristic} to obtain a good solution, even if not optimal: we chose to \textit{spread} the modulation on all weeks proportionally to $\overline{E}_{i\week}$, but with a different factor whether the \textit{reduced unit cost} $\bar{c}_{ik\week}=C_{ik}-\mu_\week$ is positive or negative. For that purpose, we split the set of weeks into two sets: $\weeks_k^+ = \{\week\in\weeks \text{ s.t. } C_{ik}-\mu_\week>0\}$ and $\weeks_k^- = \{\week\in\weeks \text{ s.t. } C_{ik}-\mu_\week\leq 0\}$. Then, when looking for the minimal reduced cost of our arc $a$, we consider two variables $\alpha_a$ and $\beta_a$ in the $[0,1]$ interval and we \textit{set} the values of $e_{\week{}}$ from these variables: \begin{eqnarray*} &&\forall \week \in \weeksof{a}\ \cap\ \weeks_{k_n}^+,\ \ e_{\week{}}=\alpha_a\ .\ \overline{E}_{i\week} \\ &&\forall \week \in \weeksof{a}\ \cap\ \weeks_{k_n}^-,\ \ e_{\week{}}=\beta_a\ .\ \overline{E}_{i\week} \end{eqnarray*} The problem becomes: find the factors $\alpha_a$ and $\beta_a$ that minimize the expression: \[ \alpha_a \left(\sum_{\week \in \weeksof{a}\ \cap\ \weeks_{k_n}^+} (C_{ik_n}-\mu_\week)\overline{E}_{i\week} \right) + \beta_a \left(\sum_{\week \in \weeksof{a}\ \cap\ \weeks_{k_n}^-} (C_{ik_n}-\mu_\week)\overline{E}_{i\week} \right) \] The parenthesized expressions do not have to be computed for each arc. Instead, we can compute them \textit{in advance} for each $(k,[\week,\week'[)_{\week<\week'}$ couple. Similarly, the constraints (\ref{eq:WAlimit_prod_i}) to (\ref{eq:WAmodulation}) can be rewritten as linear constraints on $\alpha_a$ and $\beta_a$ and we save a lot of computing overhead by precomputing $\sum_{\week \in [\week,\week'[\ \cap\ \weeks_{k}^-} \overline{E}_{i\week}$ for each couple $(k,[\week,\week'[)_{\week<\week'}$ and their counterparts on $ \weeks_{k}^+$. As it reduces the number of variables to a couple, this heuristic has proven to be more efficient and effective than looking for the exact solution, when placed in the context of this column generation problem~: even if the minimum reduced cost is not reached by this computation, the induced time savings allow the production of a greater set of columns amongst which better combinations can be found. \subsection{Obtaining dual variables $\mu_\week$ and $u_i$}\label{subsec:getDual} At each iteration we solved $L\bar{F}$ to optimality and the dual variables $\mu_\week$ and $u_i$ are respectively associated with constraints (\ref{eq:Fw_dem}) and (\ref{eq:Fw_setpar}). For the first run the variables are arbitrarily chosen as follows: $u_i=0$ for all $i \in I$ and $\mu_\week = \max_{i \in I, k \in [0,K-1] } \{ C_{ik} + 1\}$ so that $\weeks_k^+ = \weeks$ for all power plants. This will create a first pool of production plans where nuclear power plants are used at their maximal capacity. \section{Solving the final master problem}%{Fixing dates for outages and BO stock-level periods solving an aggregate problem by Column Generation} \label{step1} \label{sec:solving_the_master_problem} %\subsection{Computing the lower bound} %The Linear relaxation of the Restricted Master Problem $LRMP$ is to be solved at each iteration of the process. This LP-solving requires to find for every nuclear plant $i$ a production plan verifying local constraints [CT3-CT13] with minimum reduced cost, i.e. to find for every nuclear plant $i=0,\ldots,I-1$: %Moreover, to each production plan $\ell \in \lc{P}_i$ is associated $c_\ell$ which is the total cost of refueling plant $i$ (reduced by the expected value of remaining fuel). Let $\lch{P}$ be a subset of production plans. We call Restricted Master problem ($RM$) associated with $\lch{P}$ the above formulation $F$ restricted to columns of $\lch{P}$ and without constraints (\ref{eq:SP_cliques}). When relaxing the integrality constraints we obtain the linear relaxation ($RM$): %As the original problem can be decomposed by the $I$ nuclear plants we reformulate the problem as a ({\it aggregate}) set partitioning like Master Problem (MP): %\begin{align} %%\label{eqLab:LMRP} %(RM) \quad z(RM) = \min &\sum\limits_{i=0}^{I-1} \sum\limits_{\ell \in \lch{P}_i} c_{\ell} y_{\ell} + \sum\limits_{j=0}^{J-1} \sum\limits_{\week{}=0}^{\wmax-1} C_j^{\week{},\bar{s}} e_{\week{}j}&\label{objLRMP} %\\ %\text{s.t.}\quad &\sum\limits_{i=0}^{I-1} \sum\limits_{\ell \in \lch{P}_i} e_{\week{} \ell} y_\ell + \sum\limits_{j=0}^{J-1} e_{\week{}j} = DEM^{\week{},\bar{s}} \quad &&\forall h=0,\ldots,\wmax-1 \label{covLRMP}\\ %&\sum\limits_{\ell \in \lch{P}_i} y_\ell = 1 \quad &&\forall i=0,\ldots,I-1 \label{1iLRMP}\\ %&EMIN_j^{\week{},\bar{s}} \leq e_{\week{}j} \leq EMAX_j^{\week{},\bar{s}} \quad && \forall j=0,\ldots,J-1, \ h=0,\ldots,\wmax-1 \label{limitLRMP}\\ %%&\text{Constraints [CT14,...,CT21]} && \label{cliquesF}\\ %&y_\ell \geq 0, \quad && \forall \ell \in \lch{P}, \label{binLRMP}\\ %&e_{\week{}j} \geq 0 \quad && \forall j=0,\ldots,J-1, \ h=0,\ldots,\wmax-1 \label{positLRMP} %\end{align} % %The initial subset $\lch{P}$ of production plans of the Restricted Master Problem is empty at the beginning, therefore there are no $y_\ell$ variables, but only $e_{\week{}j}$ variables in the first program. %%The pricing problem can be modeled as a shortest path problem on a specific graph as explained in \cite{EcoLotSize}, \cite{DynLotSize} and \cite{IntProg}. The main difficulty of the method is how to build the graph. % %As we said above the pricing subproblem and the associated graphs for the shortest path computations are described in section \ref{sec:pricing}. Due to the time limit imposed in the competition, an approximated pricing subproblem is solved at each iteration: the columns that are added have to be supported by the specific graph whose construction is explained in subsection \ref{ssec:graph}. Hence $RL\bar{F}$ is not solved to optimality, but the Approximated Linear Relaxed Master Problem ($ARM$) which is the projection of $RM$ on the space where new columns have to be supported by these graphs, is solved to optimality. When no more negative reduced cost columns are found, it means that the $ALRMP$ is solved to optimality. %%At each iteration, columns with negative reduced cost are added to the RMP and the process is naturally stopped when, for all plants $i$, no production plan with negative reduced cost is found, which means that LF has been solved to optimality. % %%\subsection{Fixing the date of outages and BO-periods} %%\label{subsec:MIP} % Once $L\bar{F}$ has been solved to optimality by the Column Generation process, the following last restricted $\bar{F}$ is solved to optimality \begin{align} (R\bar{F}) \quad z(\bar{F}) = \min &\sum\limits_{i \in I} \sum\limits_{\ell \in \lch{P}_i} c_{\ell} y_{\ell} + \sum_{\week \in W} C_{j\week\bar{s}} e_{\week{}j} & \label{eq:RFw_objF}\\ \text{s.t.}\quad &\sum\limits_{i \in I} \sum\limits_{\ell \in \lch{P}_i} y_{\ell} e_{\week\ell} + \sum_{j \in J} e_{\week{}j} = DEM_{\week\bar{s}} & \forall \week \in W \label{eq:RFw_dem}\\ &\sum\limits_{\ell \in \lch{P}_i} y_{\ell} = 1 &\forall i=1,\ldots,I &\label{eq:RFw_setpar}\\ & \underline{E}_{j\week\bar{s}} \leq e_{\week{}j} \leq \overline{E}_{j\week\bar{s}} & \forall j \in J, \forall \week \in \weeks & \label{eq:RFw_limit_prod_j} \\ & \sum\limits_{i=1}^{I} \sum\limits_{\ell \in \lch{P}_i} a_{i\week} y_\ell \leq \alpha_G & \forall G \in \mathcal{G} & \label{eq:RFw_cliques}\\ &y_{\ell} \in \{0,1\} & \forall \ell \in \lch{P} &\label{eq:RFw_binF}\\ &e_{\week{}j} \geq 0 & \forall j \in J, \forall \week \in \weeks & \label{eq:RFw_positF} \end{align} \noindent where $\lch{P}$ contains the subset of production plans generated from the beginning. This $MIP$ includes constraints (\ref{eq:RFw_cliques}) which are the incompatibility constraints for some subsets $\lch{P}$ of production plans, corresponding to constraints [CT14-CT21] on the scheduling of outages. These constraints are generated dynamically during the column generation process. All of them are based on the same principle: given a set of columns, it is possible to build an incompatibility graph and we look for the maximal clique on this graph. The graph has a node for each column and an edge between any pair of incompatible columns (i.e. columns that cannot be at the same time in the optimal solution due to constraints [CT14-CT21]). If $\ell \in \lc{P}_i$ denotes the column selected for each nuclear plant $i$ in the optimal $MIP$ solution, then the variables $y_{i\week k}$, $z_{itks}$ and $\xi_{itks}$, representing the dates of outages, the dates when the stock reaches the $BO$-level and the corresponding imposed decreasing profile for the following weeks, are fixed for the rest of the whole method. \section{Columns diversification}\label{sec:diversification} The main difficulty encountered with the method described so far comes when solving the $MIP$ for column selection since constraints \texttt{[CT14-21]} are not part of the column generation process. We are not able to find a proper subset of columns that satisfies these constraints, most of the time. To overcome this problem, we used an alternate scheme for weighting the columns: instead of using weights $\weight_a$ given by equations \eqref{eq:weightBOarc} and (\ref{eq:waobj}), we compute and use weights $\weight'_a$ in the following way: \begin{equation} \label{eq:diversificationWeights} \weight'_a=\left\{ \begin{array}[c]{ll} 0 & \text{ for \texttt{BO} arcs}\\ {nbConflicts(a,col)} & \text{ for \texttt{OUTAGE} arcs}\\ \end{array}\right. \end{equation} \noindent where $nbConflicts(a,cols)$ heuristically figures out the current number of (column, constraint) couples that would conflict with a new column going through arc $a$. More precisely, we compute $\weight'_a$ with the following algorithm: \begin{algorithm}[H] \label{algo:diversification} \caption{Alternate weight computation during a \emph{diversification} pass} \begin{algorithmic} \STATE $\weight'_a \leftarrow 0$ \FORALL{existing column $\spi \in \union_i\lc{P}_i$} \FORALL{constraint $m$ of type \texttt{[CT14]} to \texttt{[CT18]}} \IF{column $\spi$ conflicts with an outage of power plant $i_a$ at week $\week_a$ by constraint $m$} \STATE $\weight'_a \leftarrow \weight'_a + 1$ \ENDIF \ENDFOR \FORALL{constraint $m$ of type \texttt{[CT19]}} \IF{for some outage $k$ of $\spi$, $(i_a,k_a)$ and $(i_\spi,k)$ are part of set $A_m$ \AND $\week_a$ is interval $[\week_k+L_{i_\spi k m}-L_{i_a k_a m},\week_k+L_{i_\spi k m}-L_{i_a k_a m}+TU_{i_\spi k m}]$} \STATE $\weight'_a \leftarrow \weight'_a + \frac{1}{Q_m}$ \ENDIF \ENDFOR \FORALL{constraint $m$ of type \texttt{[CT20]} concerning week $\week_m$} \IF{$\week_m \in [\week_a,\week_a+DA_{i_a k_a}] \cap [\week_k,\week_k+DA_{i_\spi k}]$ for some outage $k$ of column $\spi$} \STATE $\weight'_a \leftarrow \weight'_a + \frac{1}{N_m(\week_m)}$ \ENDIF \ENDFOR \ENDFOR \FORALL{constraint $m$ of type \texttt{[CT21]} concerning interval $IT_m$ such that $i_a \in C_m$} \STATE $\weight'_a \leftarrow \weight'_a + \max_{\week\in[\week_a,\week_a+DA_{i_a k_a}]\cap IT_m} \left[ \sum_{(\spi,k) \in \union_i\left(\lc{P}_i\times K_{i}\right) \ /\ i_a\neq i_\spi\ \land\ i_\spi\in C_m\ \land\ \week \in [\week_k,\week_k+DA_{i_\spi k}]\cap IT_m} \ \frac{\overline{P}_{i_\spi \week}}{IMAX_m} \right]$ \ENDFOR \end{algorithmic} \end{algorithm} Note that $\weight'_a$ only depends on the power plant $i_a$, the week of outage $\week_a$ and the campaign $k_a$~; furthermore, when new columns are generated, they only modify weights $\weight'_a$ by an additive value, independently of other columns. These facts allow major factorizations in the algorithm. This method is called each time negative reduced cost columns are found through finding the shortest path on the former graph. Then, when no more negative reduced cost columns are found with the shortest path algorithm, only the diversification method is processed and the minimum path algorithm is not ran unless zero reduced cost columns are added by diversification. Indeed, their addition could modify the value of the dual variables though not decreasing the value of the dual problem. Note that this diversification method creates positive reduced cost columns that are added if not already present. \section{Calculate the reload and the production quantities} \label{sec:ReloadProduction} Once the dates for outages and BO-profiles have been fixed as described in section \ref{sec:solving_the_master_problem}, it is still possible to improve the value of the optimal solution by looking for the optimal value of the remaining variables which are continuous. In other words, we have to determine the optimal reload quantities $r_{ik}$ during each outage $k$ for every nuclear plant $i$, and the optimal production values $p_{jts}$ and $p_{its}$ for non-nuclear and nuclear plants, and finally the optimal stock values $p_{its}$ for time steps $t$ and scenario $s$ such that the stock level is over the BO threshold. \subsection{LP-solving for fixing reload quantities} The reload quantities $r_{ik}$ are determined by solving a partially-disaggregated LP where every scenario $s$ is considered instead of the average scenario $\bar{s}$, but time steps $t$ are still aggregated into weeks $\week$ so as to accelerate solving. The optimal reload quantities $r_{ik}$ are found by solving a linear program ($\plu$) with CPLEX, where non-negative decision variables are reload variables $r_{ik}$, production variables $e_{j\week s}$ and $e_{i\week s}$, and stock variables $x_{i\weeks s}$. The constraints in this LP are all constraints [CT1-CT12], but for weeks $\week$ instead of time steps $t$. Let us define a binary value $O_{iw}$ for power plant $i \in I$ over a week $\week \in \weeks$ which is equal to $0$ if week $w$ is an outage period for power plant $i$, $1$ otherwise, and $Y^{\#}_{ik}$ is the value of $Y_{ik}$ fixed at the end of Phase 1. The two former values are known now outage periods have been fixed. Let us also define the maximum production for power plant in $i \in I$ over week $\week \in \weeks$ as \[ EMAX_{\week i}= \sum_{t\in \mathit{\timesteps{\week}}} \mathit{PMAX}_i^t \ .\ D^t \] \noindent and the maximum and minimum productions for power plant in $j \in J$ and scenario $s \in S$ over the week $w \in \weeks$ as \[ EMAX_{\week j s}= \sum_{t\in \mathit{\timesteps{\week}}} \overline{P}_{jts} \ .\ D^t \] \noindent and \[ EMIN_{\week j s}= \sum_{t\in \mathit{\timesteps{\week}}} \underline{P}_{jts} \ .\ D^t \] \noindent and \[ C_{jws} = \frac{W}{T} \sum\limits_{t \in w} C_{jts}\\ \] the mean cost of thermal energy production over a week $w \in W$ for $j \in J$ and $s \in S$. The objective function is thus the following: % \begin{align} % &&z(\plu) = & \min \sum_{i \in I} \sum_{k \in K} C_{ik} r(i,k) + \frac{1}{S} \sum_s ( \sum_h \sum_j C_{j}^{\week{}s} e(j,h,s) -\sum_i C_{i,\wmax+1} x(i,\wmax{},s)) \nonumber % \end{align} \begin{align} z(\plu) = & \min \sum_{i \in I} \sum_{k \in K} C_{ik} r_{ik} + \frac{1}{S} \sum_{s \in S} ( \sum_{w \in W} \sum_{j \in J} C_{j\week s} e_{j\week s} -\sum_{i \in I} C_i^f x_{is}^f ) \nonumber\\ %\quad z = \min &\sum\limits_{i\in I} \sum_{k \in K} C_{ik} r_{ik} + \frac{1}{|S|}\sum\limits_{s \in S}\sum_{t \in T}\sum_{j \in J} C_{jts} p_{jts} - \sum_{i \in I} C_i^f x_{is}^f & \label{eq:obj_CR}\\ \text{s.t.} \quad & \sum\limits_{i \in I} \sum\limits_{k \in K_i(\weekof{t})} p_{ikts} + \sum\limits_{j \in J} p_{jts} = DEM_{ts}, &\forall t \in T, \forall s \in S & \label{eqRel:equil_dem} \\ & 0 \leq e_{i\week s} \leq EMAX_{wi}O_{iw}, &\forall i \in I, \week \in \weeks, \forall s \in S & \label{eqRel:limit_prod_i} \\ & EMIN_{wjs} \leq e_{j\week s} \leq EMAX_{wjs}, &\forall j \in J, t \in T, s \in S & \label{eqRel:limit_prod_j} \\ %& \sum\limits_{(i\week k) \in G} y_{i\week k} \leq \alpha_G, &\forall G \in \mathcal{G} & \label{eq:cliques_CR}\\ %& Y_{ik}= \sum\limits_{w \in W_{ik}} y_{iwk} &\forall i \in I, k \in K_i & \label{eq:Yik}\\ & \underline{R}_{ik} Y^_{ik} \leq r_{ik} \leq \overline{R}_{ik} Y^{\#}_{ik} &\forall i \in I, k \in K_i &\label{eqRel:r_minmax}\\ %& Y^\#_{i,k-1} \geq Y_{i,k} &\forall i \in I, \forall k \in K_i \setminus \{1\} &\label{eq:y_kk}\\ & x^1_{i,k+1,s} = x^2_{iks} - \sum_{w \in W'_{ik}} e_{i\week s} &\forall i \in I, k \in K_i s \in S & \label{eqRel:xx'} \\ %\mathtt{[CT10,CT9]} & x^2_{iks} = \left[\left( \frac{Q_{i,k} -1}{Q_{i,k}}\right)( x^1_{iks} - \mathtt{BO}_{i,k-1}) + r_{ik} + \mathtt{BO}_{i,k} \right]Y^{\#}_{ik} + x^1_{iks}(1-Y^{\#}_{ik}) &\forall i \in I, k \in K_i, s \in S & \label{eqRel:x'x} \\ & x^1_{iks} \leq A_{ik} + (1 - Y^{\#}_{i,k})(M - A_{ik}), &\forall i \in I, k \in K_i, s \in S & \label{eqRel:Amax} \\ & x^2_{iks} \leq S_{ik} + (1 - Y^{\#}_{i,k})(M - S_{ik}), &\forall i \in I, k \in K_i, s \in S & & \label{eqRel:Smax} \\ & x^f_{is} = \sum_{k \in K_i} x_{i,k+1,s}(Y^{\#}_{ik} - Y^{\#}_{i,k+1}), &\forall i \in I, s \in S & \label{eqRel:stock_final} \\ %&(BO_{ik} - x_{itks})z_{ikts} \geq 0, & \forall i \in I, \forall k \in K_i,\forall t \in \timeline, \forall s \in S & \label{eq:BO_level_uno} \\ %&(x_{itks} - BO_{ik})(1-z_{ikts}) \geq 0, & \forall i \in I, \forall k \in K_i,\forall t \in \timeline, \forall s \in S & \label{eq:BO_level_due} \\ %&(PB_{ik}(x_{itks})\overline{P}_{it}D_t - x_{itks})\xi_{itks} \leq 0 & \forall i \in I, \forall k \in K_i, \forall t \in \timeline, \forall s \in S & \label{eq:activer_xi_uno} \\ %&(x_{itks} - PB_{ik}(x_{itks})\overline{P}_{it}D_t)(1-\xi_{itks}) \leq 0 & \forall i \in I, \forall k \in K_i, \forall t \in \timeline, \forall s \in S & \label{eq:activer_xi_due} \\ %&(1 - y_{i\week k})(p_{itks}-PB_{ik}(x_{itks})\overline{P}_{it}\xi_{itks})z_{itks} = 0 & \forall i \in I, \forall k \in K_i,\forall t \in \timeline, \forall s \in S & \label{eq:prodImposed_BO} \\ %& x_{itks} = x^d_{iks} + \sum\limits_{\tau=1}^{t} p_{i\tau ks} &\forall i \in I, \forall k \in K_i, \forall t \in \timeline, \forall s \in S & \label{eq:relationxt_x} \\ %& y_{i\week k}, Y_{ik} \in \{0,1\} & & \label{eq:binary_variable} \\ & \sum\limits_{\week\in W'_{ik}} ( \overline{E}_{i\week}-e_{i\week s} ) \leq \overline{M}_{ik} & \forall i \in I, k \in K_i, s \in S & \label{eqRel:Pmodulation} \\ & e_{j\week s},e_{i\week s},r_{ik}, x^1_{iks}, x^2_{iks}, x^f_{is} \geq 0 & & \label{eqRel:continuous_var} \end{align} This linear problem $\plu$ is solved for the different combinations of columns that produce the optimal value of the Master Problem $F$ and the combination that gives the lowest value is kept to go to the next step where the original problem will be solved with outage periods and reload quantities fixed. \subsection{Disaggregate LP-solving for every scenario and initial time steps} At this final step, for every scenario $s \in S$, we go back to time steps $t$. Therefore, $S$ Linear Programs ($\pld$) are solved where decision variables are: \begin{itemize} \item[(i)] the production variables $p_{jts}$ and $p_{its}$ for non-nuclear and nuclear plants, respectively, \item[(ii)] the stock variables $x_{its}$ for time steps $t$ such that the stock level is over the $BO$ threshold (for the other time steps, the stock level is already fixed as imposed by the decreasing profile) \end{itemize} Constraints are the same as in the previous step, i.e. constraints [CT1-CT12], but for time steps $t$ now, i.e. original constraints [CT1-CT12]. The objective functions are the total production cost for each scenario $s$, reduced by the objective value of remaining fuel: \begin{align} && z(\pld) = \min & \sum_t \sum_j C_{jts} p_{jts} D^t -\sum_i C^f_{i} x^f_{is} \nonumber \end{align} The final objective value $z$ is the refueling cost and the average total production cost reduced by the objective value of remaining fuel over all $S$ scenarios: \begin{equation} \label{eq:z} z = \min \sum_i \sum_k C_{ik} r_{ik} + \frac{1}{S} \sum_s \left( \sum_t \sum_j C_{jts} p_{jts} D^t -\sum_i C^f_{i} x^f_{is}\right) \end{equation} %\begin{align} %&& z = \min & \sum_i \sum_k C_{ik} r(i,k) + \frac{1}{S} \sum_s ( \sum_t \sum_j C_{j}^{ts} p(j,t,s) D^t -\sum_i C_{i,T+1} x(i,T,s)) \nonumber %\end{align} \section{Computational Results} \label{sec:ComputationalResults} All algorithms described in this paper were coded in C++ and compiled with the gcc (Debian 4.4.4-8) 4.4.5 compiler. CPLEX 12.1 was used as the LP solver and as the integer programming solver. The experiments were performed on a Squeeze-x64-64 (image based on Debian version sid for AMD64/EM64T) Bi AMD Opteron Dual core 3,2 Ghz with 16 GB of physical memory and 16 GB of swap memory. In our computational experiments, we considered 16 instances given by EDF during the Challenge (see \cite{sujetOnline}). Almost all of them ask for a planification of $6$ campaigns for a time period of $5$ years. The instances have been partitioned into the following three classes. \textit{Class 1.} A set of $6$ instances (A0-A5) used for the qualification phase. In these instances the number of different scenarios is at most $30$ for about $20$ nuclear plants and $30$ other thermic power plants. The time step is the day and there are less than $27$ constraints linking the nuclear power plants. \textit{Class 2.} A set of 5 instances (B6-B10) used for the second phase. In these instances the number of different scenarios considered can rise up to $121$. There are about $50$ nuclear plants and $20$ other thermic power plants. The time step is $8$ hours and there are about $140$ constraints linking the nuclear power plants. There are between $114$ and $235$ constraints [CT13] on the dates of the outages. \textit{Class 3.} A set of 5 instances (X11-X15) used for the final ranking with the $5$ former instances. They were not released until the end of the challenge. In these instances the number of different scenarios are $50$ or $121$. There are about $50$ nuclear plants and $20$ other thermic power plants. The time step is $8$ hours and there are about $140$ constraints linking the nuclear power plants. There are between $207$ and $260$ constraints [CT13] on the dates of the outages, which is the only noticeable difference with the \textit{class 2} instances. \begin{table}[tbph] \begin{center} \input{ROADEF_results.tex} \end{center} \caption{Results of the ROADEF challenge} \label{tableROADEF} \end{table} In table \ref{tableROADEF}, the best results that were obtained and published on the website of the Challenge are given in the first column whereas the second column displays the ones obtained by the current method. All teams were ranked by the Challenge organizers and with the former results our method ranked $4^{th}$ although we were penalized by the fact that the method was not able to give answers for three of the instances. It is interesting to point out that among the first five methods proposed, our method is the only one, as far as we know, having some mathematical programming foundations. The third column gives the best results we found since then, after the few bugs were corrected. The last column gives the parameter $\rho$ that was chosen to obtain the former result, i.e., the parameter used to build the graphs to solve the $L\bar{F}$. % Indeed the pricing subproblem is heuristically solved in fact, so the found lower bound that appears in the following table is not the real lower bound on the problem (\label{MP}). Several computational tests have been done with different value for this parameters and sometime one instance cannot be solved at the end of the method: it is the case in table \ref{table05} where data8 is not solved by the method. %Recall here that the problem solved at each iteration - if necessary - is the linear relaxation of the Master Problem (\ref{MP}) %where the constraints (\ref{binCMP}) are replaced by the following ones: %$$ x_r \geq 0 \mbox{ } \forall r \in \cup_{i \in I} R_i \label{relCMP}$$ Tables \ref{tableResults} and \ref{table05} show several results about the method: \begin{itemize} \item The first column specifies the considered instance~; \item The second column shows the optimal result of the $LR\bar{F}$~; \item The third column gives the iteration at which the former value is reached: the following iterations do not decrease the former value~; \item The fourth column shows the total number of iterations of the algorithm~; \item The fifth column gives the final number of plans that have been computed at the end~; \item The sixth column gives the final number of constraints [CT14,...,CT21] that have been added to the \mip~; \item The seventh column gives the value of MIP when it is called at the end of the column generation process~; \item The eigth column gives the value of $\plu$~; \item The ninth column gives the result of $z$ defined in equation (\ref{eq:z})~; \item The last column gives the elapsed time for the whole algorithm \end{itemize} \begin{table}[tbph] \begin{center} \input{tableauFinal.tex} \end{center} \caption{Performance of the overall method with $\rho = 0.9$ } \label{tableResults} \end{table} \begin{table}[tbph] \begin{center} \input{tableauFinal30.tex} \end{center} \caption{Performance of the overall method with $\rho = 3$ } \label{table05} \end{table} The value of the relaxation of the Master Problem decreases as negative reduced cost columns are added. As the lowest value is reached in a few dozens of iterations and the columns added afterwards do not change this value, the optimal value of the Approximated Relaxed Pricing Subproblem is reached then. Note that the number of lines (Cs Nb) and variables (Pl Nb) has the same order of magnitude. It means that both the rows and columns generation processes are necessary to find quickly the integer solution of the $ALRMP$. Table \ref{tableTime} gives CPU time in seconds of the method: \begin{itemize} \item The first column specifies the considered instance~; \item The second column shows the time to read the input and to write the output~; \item The third column gives the time to build the graph for the pricing problem~; \item The fourth column shows the total time of column generation~; \item The fifth column gives the time to solve the integer master problem obtained at the end of the column generation process~; \item The sixth column gives the time of the reload phase~; \item The seventh column gives the time of the production phase~; \item The last column gives the elapsed time for the whole algorithm \end{itemize} \begin{table}[tbph] \begin{center} \input{tableauTemps.tex} \end{center} \caption{CPU Time (s) } \label{tableTime} \end{table} Actually the whole method never reaches the set time limit ($3600$ seconds). It can be used, for example, as a very good initial solution which could be a starting point for a local search procedure. The most time-consuming phase is the MIP phase at the end of the column generation process. The reload phase is more time consuming than the production phase due to the multiple resolutions for the different combinations of columns that produce the \mip. We can notice that the time taken by the input and output is not negligible due to the large size of the files to read and write. \section{Conclusion} \label{sec:conclusion} We presented a mathematical formulation for the EDF (\'Electricit\'e De France) electricity production planning problem submitted as subject of the ROADEF/EURO 2010 Challenge, and a heuristic method based on column generation. The method proceeds in two stages~: first we fix dates for outages for each nuclear plant on aggregated data, and in a second step we determine the reload and production quantities for each plant and each scenario. Computational results show the effectiveness of the proposed approach, given that the method was far from using all the time allowed in the Challenge. Therefore further work could consist in generating more sets of schedules for outages using local search on the set of outage dates found in the first phase of the method. \section*{Acknowledgment} The authors of this paper would like to thank Aristide Mingozzi and Guillaume Turri for their useful help and/or comments they gave at the beginning of this project. \bibliographystyle{plain} \bibliography{Paintshop} \end{document} %By this way, it is expected that this heuristic fixing of dates for outages will provide production and refueling campaigns whose total cost is relatively close to the optimal cost. %\begin{comment} %We first give in subsection 1.1 the context and objective of the problem provided by EDF as it is presented by the company in the official document \cite{EDF}, %then we review the research literature about electricity production planning. %\subsection{State of the art of electricity production planning} %\end{comment} % Optimization in energy is a large subject. Several different problems arise in different topics. With the increasing price of raw materials and the growing awareness of environment concerns energy consumption must be thought in a less greedy way in order to save resources. For years several researchers have been trying to work on that topic as gathered in the book \cite{NRJbook}. Subjects can be as various as wind energy optimization (\cite{wind}), smart grid design or production planning (\cite{multistageOpt}). \begin{comment} \subsection{Computing new dual variables $\mu_\week$} At each iteration when we solve the restricted MP, the Simplex method provides, at optimality, dual variables $\mu_\week$ associated with each week $\week$, and dual variables $u_i$ associated with every nuclear plant $i$. These dual variables are used to update the reduced costs on the arcs $a_i \in A_i$ of every graph $G_i$. $-u_i$ is added once for all at the sink node of graph $G_i$. For an arc $a_i = (\week{}xk,\week{}'x'k') \in A_i$ linking a state node $\week{}xk$ to another state node $\week'x'k'$, we have \begin{itemize} \item[(i)] If $(\week{}xk,\week'x'k')$ is an 'outage' arc then the reduced cost of arc $(\week{}x,\week'x'k')$ is $c(\week{}xk,\week'x'k') = -\sum_{l=\week+1}^{\week{}'} P^*_{l} \mu_l + C_{ik}(x_{\week{}+1}-C10(x))$ where: \begin{itemize} \item values $P^*_{l}$ represent the production on weeks $l=\week,\ldots,\week'-1$ \item $x_{\week{}+1}$ represents the stock after the reload \item $k$ is the cycle number of the outage \item and $C_{ik}$ is the proportional refueling cost. \end{itemize} $P^*_{l}$ and $x_{\week{}+1}$ are computed so as to minimize $c(\week{}xk,\week'x'k')$ over all possible production values $P_l$, $l=\week,\ldots,\week'-1$ and reload quantities $x_{\week{}+1}-C10(x)$ that achieve an overall stock variation of $x-x'$. This is fastly done by sorting the $\mu_l$ dual variables in decreasing order and setting the production values to the maximum allowed production in that order. \item[(ii)] If $(\week{}xk,\week'x'k')$ is a 'BO' arc then the reduced cost of arc $(\week{}xk,\week'x'k')$ is computed as $-\sum\limits_{l=\week}^{\week{}'-1} P_{l} \mu_l$, taking the imposed production profile as $P_{l}$ without optimizing with the allowed $\epsilon$ deviation. \end{itemize} By construction, the sum of the costs on all the arcs of a path is the current reduced cost of the path as computed in (\ref{subp}) (the expected value of remaining fuel is taken into account at the end of the path). \end{comment} % is very intricate because of the numerous possible stock levels. % We use an approximation then by discretizing the stock levels and % solve the pricing subproblem for a plant $i$ in a graph % $G_i=(N_i,A_i)$ where {\it state} nodes in $N_i$ represent pairs $hx$ % meaning that the stock is at level $x$ at the start of week $\week$ % (i.e. before any refill or production). The construction of graphs $G_i$ satisfies the following properties: % \begin{itemize} % \item[(i)] There are two categories of arcs $(hx,h'x')$: 'BO' arcs % modeling a sequence of production timesteps for which the remaining % fuel is under the 'BO' level, and 'outage' arcs modeling an % outage followed by a sequence of production timesteps for which the % remaining fuel is above 'BO'. % \item[(ii)] An 'outage' arc between two state nodes $hx$ and $h'x'$ necessarily satisfies the maximum modulation constraint [CT12]. % \item[(iii)] 'BO' arcs follow the imposed profile with the $\epsilon$ tolerance. % \item[(iv)] The set of discretized levels of stock $x$ is generated % during the process, so that unreachable levels are not generated for % a specific week, in a way that is described in what follows. % \end{itemize} % Let us describe more precisely item (iv). Let $hx$ be a state node. Then:\\ % - if 'outage' arcs can leave from $hx$ (depending on $x$ and $\week$) we generate a % lower bound $\underline{x}= CT10(x)+RMIN$ and an upper bound % $\bar{x}=\min(SMAX, CT10(x)+RMAX)$ for the start of the production % period~; then, we generate a lower bound % $\underline{x'} = \max(\underline{x}-\sum_{t}PMAX(t)*dt,BO)$ and an % upper bound $\bar{x'}=\max(\bar{x}-(\sum_{t}PMAX(t)*dt) + MMAX, BO)$ % where $t$ runs over the production period % and we discretize interval $[\underline{x'},\bar{x'}]$ by steps of an % order of magnitude $PMAX(t)*dt$ for generating new nodes $h'x'$.\\ % - for 'BO' arcs, we follow the required production profile with the % $\epsilon$ tolerance from $\week$ to $h'$, which can lead to two possible intervals % $[\underline{x_1'},\bar{x_1'}]$ and $[0,\bar{x_2'}]$. These two % intervals then get discretized. % {\bf Antoine : revoir cette partie + ins\'erer graphe partiel + faisceau}\\ %The program %authorizes to add positive reduced cost columns if not already present in the pool of columns. The process of adding columns %stops if it is unable to add new columns that are not present in the pool or if a time limit is reached. %Lastly the combination that has the least cost is selected with the outages and reloads fixed. The optimal production is then %computed on each scenario and the average value obtained is the one of the sixth column. %\bigskip %Dataset B %\begin{tabular}[c]{||r|c|l||c|l||} % \hline % & \multicolumn{2}{|c||}{1 thread} & \multicolumn{2}{|c||}{4 threads} \\\hline % data & Objective value & Running time & Objective value & Running time \\\hline %data6 & $8.47103*10^{10}$ & 00:28:50 & $8.470769*10^{10}$ & 00:18:00 \\ %data7 & $8.18527*10^{10}$ & 00:55:20 & $8.184168*10^{10}$ & 00:29:55 \\ %data8 & $8.30129*10^{10}$ & 00:51:11 & $8.299788*10^{10}$ & 00:34:09 \\ %data9 & $8.36119*10^{10}$ & 00:47:45 & $8.367746*10^{10}$ & 00:43:25 \\ %data10 & $7.87675*10^{10}$ & 00:44:27 & $7.873125*10^{10}$ & 00:46:52 \\ %\hline %\end{tabular} % \begin{tabular}[c]{||r|c|l|c||} % \hline % & \multicolumn{2}{|c||}{4 threads} & \\\hline % data & Objective value & Running time & Best\\\hline % % data0 & $8.73626*10^{10}$ & ?? & $8.7309853670*10^{10}$ \\ % data1 & $1.70030*10^{9}$ & ?? & $1.6953838690*10^{9}$ \\ % data2 & $1.46409*10^{9}$ & ?? & $1.4604843397*10^{9}$ \\ % data3 & $1.54710*10^{9}$ & ?? & $1.5442988894*10^{9}$ \\ % data4 & $1.12808*10^{9}$ & ?? & $1.1159135051*10^{9}$ \\ % data5 & $1.28366*10^{9}$ & ?? & $1.2582223638*10^{9}$ \\ % data6 & $8.4709341172*10^{10}$ & ?? & $8.3424716217*10^{10}$ \\ % data7 & $8.1839081672*10^{10}$ & ?? & $8.1174243138*10^{10}$ \\ % data8 & $8.2977018880*10^{10}$ & ?? & $8.1926206073*10^{10}$ \\ % data9 & $8.3339264015*10^{10}$ & ?? & $8.1750858197*10^{10}$ \\ % data10 & $7.8730244463*10^{10}$ & ?? & $7.7767024999*10^{10}$ \\ % data11 & $8.0100972234*10^{10}$ & ?? & $7.9116772289*10^{10}$ \\ % data12 & $7.8193054033*10^{10}$ & ?? & $7.7589910940*10^{10}$ \\ % data13 & $7.6855673485*10^{10}$ & ?? & $7.6449207715*10^{10}$ \\ % data14 & $7.7084931871*10^{10}$ & ?? & $7.6172998633*10^{10}$ \\ % data15 & $7.5037618871*10^{10}$ & ?? & $7.5037618871*10^{10}$ \\ % % % \hline % % \end{tabular} % & \textit{s.t. } & [CT1-CT13bis] & for $i, j, s, h$ \nonumber\\ % $\sum_j p(j,h,s) + \sum_i p(i,h,s) = DEM^{\week{},s}$ & $\text{ for } s,h$ \nonumber\\%\label{dem}\\ % % & & $e'_{j,o} + \sum_{\week{} \in o} p^s_{j,h} = e_{j,o} $ & $\forall o \in O$ \nonumber\\ % & & $e_{j,o} + { F_{j,o+1}} = e'_{j,o+1}$ & $\forall j \in J, \forall o \in O$ \nonumber\\ %\label{mf:P.CapC}\\ % % & & $\sum_{\ell \in \lc{R}_c} x_\ell \leq b$ & $\forall c \in C$, \label{clique}\\ % & & $max^s_{i,h} \geq y^s_{i,h} \geq min^s_{i,h}$ & $\forall i \in \tpp, \forall h \in H, \forall s \in \sce$ \nonumber\\ % & & $max^s_{j,h} \geq p^s_{j,h} \geq min^s_{j,h}$ & $\forall j \in \npp, \forall h \in H, \forall s \in \sce$ \nonumber % %\label{mf:P.IntX}\\ %\sum_{j \in \npp} C_{j} \sum_{\week{} \in H} x^s_{j,h} + \sum_{i \in \tpp} \sum_{\week{} \in H} c^s_{i,t} y^s_{i,t} \nonumber\\ % \begin{tabular}[l]{llll} % z(\pld)= & $\min$ & $\sum_{s \in \sce} \sum_{i \in \tpp} \sum_{t \in \timeline} c^s_{i,t} { y^s_{i,t}}$ & \nonumber\\ % & \textit{s.t. } & $\sum_{j \in \npp} { p^s_{j,t}} + \sum_{i \in \tpp} { y^s_{i,t}} = d^s_t$ & $\forall t \in \timeline, \forall % s \in \sce$ \nonumber\\ %\label{dem}\\ % & & $e'_{j,o} + \sum_{t \in o} { p^s_{j,t}} = e_{j,o}$ & $\forall o \in O$ \nonumber\\ % & & $e_{j,o} + F_{j,o+1} = e'_{j,o+1}$ & $\forall j \in J, \forall o \in O$ \nonumber\\ %\label{mf:P.CapC}\\ % % & & $\sum_{\ell \in \lc{R}_c} x_\ell \leq b$ & $\forall c \in C$, \label{clique}\\ % & & $max^s_{i,t} \geq { y^s_{i,t}} \geq min^s_{i,t}$ & $\forall i \in \tpp, \forall t \in \timeline, \forall s \in \sce$ \nonumber\\ % & & $max^s_{j,t} \geq { p^s_{j,t}} \geq min^s_{j,t}$ & $\forall j \in \npp, \forall t \in \timeline, \forall s \in \sce$ % \nonumber%\label{mf:P.IntX}\\ % \end{tabular} %, which is defined as follows %\[ %e_{\week{}\ell}=\frac{1}{S} \sum\limits_{s=0}^{S-1} \sum_{t \in week \ h} D^t.\ p_{\ell}(i,t,s)\\ %\] %REMARQUE LAURENT : enlevĂ© car notation non cohĂ©rente: p(j,t,s) = variables, de plus e_{\week{}l} n'a pas de lien avec les scenarios s et n'est pas calcule en fonction des p(i,t,s) %dans le pricing %Finally we can define also %\[ %e_{\week{}j}=\frac{1}{S} \sum\limits_{s=0}^{S-1} \sum_{t \in week \ h} D^t.\ p(j,t,s) %\] %IDEM LAURENT : enlevĂ© : e_{\week{}j} est une variable du maitre agrĂ©gĂ© sur 1 scĂ©nario ne dĂ©pendant pas des variables p(j,t,s) par scenarios.