%\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{subeqnarray} \usepackage{graphics} \usepackage{tikz} \usepackage{algorithmic,algorithm} \usepackage{pstricks} \usepackage{a4wide} \usepackage{verbatim} \usepackage{graphicx} \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} %%%%% Time step \newcommand{\ts}{T} %%%%% Fuel cost \newcommand{\fc}{r} %%%%% Thermal variable cost \newcommand{\vc}{c} %%%%% Set-partitioning production planning index \newcommand{\spi}{\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}{LP1} %%%%% Problem name \newcommand{\pld}{LP2} % ----------------------------------------------------------------------------- \def \union{\mathop{\mathrm{\cup}}} \def \equi{\Longleftrightarrow} % Insert the name of "your journal" with \journalname{journal of scheduling} % \begin{document} \title{Solving the 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 Letocart } %\authorrunning{Short form of author list} % if too long for running head \institute{F. Author \at first address \\ Tel.: +123-45-678910\\ Fax: +123-45-678910\\ \email{fauthor@example.com} % \\ % \emph{Present address:} of F. Author % if needed \and S. Author \at second address } \date{Received: date / Accepted: date} % The correct dates will be entered by the editor \maketitle %\item[${}^1$] {\it ESSEC, Av. Bernard Hirsch, BP105 95021, Cergy Pontoise, France} \\ % Email: {\it alfandari@essec.fr} %\item[${}^2$] {\it LAMSADE, Universit\'e Paris-Dauphine, 75116 Paris, France} \\ % Email: {\it jerome.monnot@lamsade.dauphine.fr} %\end{itemize} \begin{abstract} This paper presents the column generation method used by the authors to solve the EDF (Electricit\'e De France) electricity production planning problem of the ROADEF/EURO 2010 Challenge, which is to our knowledge the first-ranked method among column generation methods proposed in the Challenge, and was ranked \#4 overall. The problem consists in determining, for the set of all thermal power plants of the French electricity company, a schedule of plant outages, which are necessary for refueling and maintenance operations, and for each plant, a production plan over the whole time horizon. The average cost of the global outage and production planning, computed over a set of demand scenarios, is to be minimized. The method proceeds in three main steps. In the first step which is the most time-consuming, dates for outages are fixed for each nuclear plant by transforming the initial problem into an aggregate covering master problem with a single average scenario and reduced time steps, and solving that reformulated problem by column generation, where the pricing problem associated with each nuclear plant is a shortest path problem in an appropriate graph. In a second step, the reload level is determined at each date of an outage, considering all scenarios in the optimization at that stage. Finally in step three, production quantities between two outages are optimized for each plant and each scenario by solving independent Linear Programs. The numerical results reported are the ones presented in the Challenge. \end{abstract} {\bf Keywords}: Energy, electricity production, outage, column generation, shortest path, linear programming. \section{Introduction} %\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}). Optimization in energy is a large subject. Several problems arise in different topics. With the raise in the prices 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}. The increasing part of the energy production that is to be done through renewable energy enhances this trend. As the area of research is huge, different works were done. Bruglieri and Liberti in \cite{bioM} propose a modelisation to improve the efficiency of biomass energy production plants. Mustakerov and Borissova in \cite{eol} give tools to help to optimally deploy a grid of wind turbines on a park knowing the wind direction. Multiple works have also been done on optimally planning the energy production when several means can be used to produced it. This has been done with both a heuristic method in \cite{heuristicHybride} by Gupta {\textit{et al.}} and an exact method by Morais {\textit{et al.}} in \cite{hybrideProd}. Energy production planning is one of the most common studied topics because a lot of improvements can be achieved. Indeed in most countries energy can be produced through different ways: nuclear and other thermal energies produce the major part of the production but hydraulic energy or other renewable energies have a significant role in several countries. Producers could play with these different production means in order to satisfy the energy demand though decreasing their costs or their impact on the environment. This subject has several difficulties one of which is that the exact amount of the future demand or the exact price of raw materials is not known. Moreover there are a lot of constraints on the energy production plants: the level of water in dams has to stay above a certain threshold or manutention operations have to be executed periodically on nuclear plants, periods during which they are not producing \ldots. This problem that has been studied for years now (\cite{EDF1} for instance) was the starting idea of the EURO/ROADEF 2010 challenge. EDF power generation facilities in France stand for a total of 98.8 GW of installed capacity. The varied range of EDF facilities mixes all forms of energy: thermal (nuclear, coal, fuel oil and gas), hydraulic and other renewable energies. Most of the electricity EDF generates in France is produced by thermal power plants: 90 \% in 2008, among which 86 \% by nuclear power plants. This subject is focused on thermal power plants and especially on those that have to be repeatedly shut down for refueling and maintenance, e.g. nuclear ones. The scheduling of these outages has to comply with various constraints, regarding safety, maintenance, logistics and plant operation while it must lead to production programs with minimum costs. The proposed subject consists of modeling the production assets and finding an optimal outage schedule and includes of two dependent sub problems: determining a schedule of plant outages. The schedule must satisfy constraints in order to comply with limitations on resources which are necessary to perform refueling and maintenance operations and determining an optimal production plan to satisfy demand, i.e. a quantity of energy to produce by each plant at each time step for each scenario. Production must satisfy some technical constraints. In medium-term electricity management numerous uncertainties have to be taken into account (demand, generation units availability, spot market prices, quantities that can be bought or sold,...), which leads to the need to consider multiple scenarios. The objective is to minimize the expected cost of production. Here a column generation method (\cite{CGbook}) is proposed to solve it. %The problem proposed in the challenge is to plan the periods where nuclear plants will be stopped for several years in order to respect the demand and the constraints on nuclear plants but also on the other thermic power plants. % On energy production planning the subject has several difficulties such that the fact the exact amount of the demand is not known. Another one is inherent to the use of nuclear power plants. Indeed there are producing energy for a very law cost but there are several constraints to be respected while using them. Column generation is often used to solve large-scale optimization problems when the problem can be formulated as an integer linear program with a huge number of variables (or \textit{columns}), generally called \textit{Master Problem} (see \cite{CGsurvey} for a set of problems where column generation was applied). The column generation resolution method iteratively creates 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 is considered. If the decision variables are integer, this constraint is relaxed to have continuous decision variables. Each restricted master problem is solved to optimality and iteratively enriched with new columns output from the resolution of a \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 short or polynomial time. %This method enables when trying to solve a large-scale problem to only generate columns that could be part of the optimal one. Indeed only negative reduced cost column are added to the master problem. The method stops when no more negative reduced cost column is output from the pricing subproblem and so the original master problem is solved to optimality. The electricity production planning problem is presented in Section \ref{}. Section \ref{} describes the overall method used to solve it. Section \ref{} gives the model of the master problem that is solved for fixing dates of outages and BO stock-level periods as first stage. Section \ref{} describes the pricing subproblem and its resolution method. Section \ref{} exposes the second stage of the method for fixing reload and production quantities for each plant. Section \ref{} gathers the results obtained on the test instances. Section \ref{} gives some perspectives about this method. \section{Problem statement} \label{sec:problemstatement} The model proposed by the EDF company is the following. For sake of consistency, the notations used in the ROADEF/EURO 2010 Challenge document \cite{EDF} are almost totally kept in this paper. The various indexes used are the index of scenarios $s=0,\ldots,S-1$, the index of time steps $t=0,\ldots,T-1$, the index of weeks $h=0,\ldots,H-1$, the index of non-nuclear plants $j=0,\ldots,J-1$, the index of nuclear plants $i=0,\ldots,I-1$ and the index $k=1,\ldots,K-1$ of production cycles for nuclear plants, $K-1$ being the maximum number of cycles that can start during the considered period. The decision variables are: \begin{itemize} \item $ha(i,k)$ : week of cycle $k$ during which the nuclear plant $i$ goes offline (i.e. week of decoupling) \item $p(j,t,s)$ : production of non-nuclear plant $j$ during the time step $t$ over scenario $s$ \item $p(i,t,s)$ : production of nuclear plant $i$ during the time step $t$ over scenario $s$ \item $r(i,k)$ : reload performed during the outage of cycle $k$ of nuclear plant $i$ \item $x(i,t,s)$ : stock of fuel of nuclear power plant $i$ at the beginning of time step $t$ over scenario $s$ \item $ec(i,k)$ : set of time steps composing the production campaign of cycle $k$ of nuclear plant $i$ \item $ea(i,k)$ : set of weeks composing the outage of cycle $k$ of plant $i$ \end{itemize} The objective function is: \begin{equation} \min \sum_i \sum_k C_{ik} r(i,k) + {1 \over 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)) \end{equation} where $C_{ik}$ is the proportional cost of fuel during cycle $k$ for nuclear plant $i$, $C_j^{ts}$ is the proportional production cost for time $t$ and scenario $s$ for non-nuclear plant $j$, and $D^t$ is the length of time step $t$. As the framework of this paper is column generation, we decompose the set of constraints into the {\it global} constraints that link different plants $i$ and $i'$, $i \not= i'$, and the {\it local} constraints, that are specific to a single plant $i$. \subsection{Global constraints linking different plants} The global constraints that will be belong to the master problem in the column generation scheme are the following: \begin{itemize} \item[(i)] Demand satisfaction: $$\text[CT1] \sum_j p(j,t,s) + \sum_i p(i,t,s) = DEM^{t,s} \text{ for } s,t $$ \item[(ii)] Constraints on outages of different plants $i$ and $i'$, concerning variables $ha(i,k)$: \begin{itemize} \item[[CT14]] Minimum spacing and maximum overlapping between outages of $i$ and $i'$ \item[[CT15]] Minimum spacing and maximum overlapping between outages between a specific period \item[[CT16]] Minimum spacing between decoupling dates \item[[CT17]] Minimum spacing between coupling dates \item[[CT18]] Minimum spacing between coupling and decoupling dates \item[[CT19]] Resource constraints on outages \item[[CT20]] Maximum number of overlapping outages during a week \item[[CT21]] Maximum offline power capacity of a set of power plants during a time period \end{itemize} \end{itemize} All constraints [CT14-CT21] involving variables of dates of outages $ha(i,k)$ are not detailed in this paper, as they will be put into the master problem only at the end of the column generation process. We rather report the local constraints on production and stock as we shall refer to them in the definition of the pricing subproblem and its associated graph. \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{itemize} \item[[CT3]] Offline power \item[[CT4]] Minimum power \item[[CT5]] Maximum power before activation of imposition of power profile constraint $(t \in ec(i,k)) \text{ and } x(i,t,s)\geq BO_{ik} \Rightarrow p(i,t,s) \leq PMAX_i^t \text{ for } s,t,i,k $ where $BO_{ik}$ is a bound on stock of fuel that activates the imposition of power profile during cycle $k$, and $PMAX_i^t$ is the maximum power of plant $i$ during time step $t$ \item[[CT6]] Maximum power after activation of imposition of power level constraint : \begin{eqnarray*} \text{ for } s,t,i,k, &\text{ if }& t\in ec(i,k) \ \and x(i,t,s)h$ so that there exists a production plan leading the powerplant from state $n$ to state $n'$ without breaking 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, h_n, k_n)$} \label{algo:beam} \begin{algorithmic}[1] \FOR{week $h=h_n+1$ \TO $H-1$} \STATE \label{algostep:possibleStock} compute the possible stock interval $[x_{min_{n,h}},x_{max_{n,h}}]$ of the powerplant $i$ at the beginning of week $h$ \IF{the outage of campaign $k_n+1$ is allowed at week $h$ according to constraint \texttt{CT13}} \STATE compute the intersection of $[x_{min_{n,h}},x_{max_{n,h}}]$ with $[0,\mathtt{AMAX}_{i,k_n+1}]$ \STATE \label{algostep:nouvelOutage} create new nodes $n'=(\mathtt{OUTAGE}, x_{n'}, h, k_n+1)$ with $x_{n'}$ in this intersection (constraint \texttt{CT11}) \ENDIF \IF{$\mathtt{BO}_{i,k_n} \in [x_{min_{n,h}},x_{max_{n,h}}]$} \STATE \label{algostep:nouveauBO} create new nodes $n'=(\mathtt{BO}, x_{n'}, h-1, k_n)$ with $x_{n'} \in [x_{min_{n,h-1}},x_{max_{n,h-1}}]$ computed so that if the powerplant follows the imposed profile of \texttt{CT6} during week $h-1$, its stock will fall under $\mathtt{BO}_{i,k_n}$ at the beginning of week $h$ \ENDIF \ENDFOR \end{algorithmic} \end{algorithm} For steps \ref{algostep:nouvelOutage} and \ref{algostep:nouveauBO} of algorithm \ref{algo:beam}, we discretize the interval of possible stocks of fuel at week $h$ $[x_{\min_{n,h}}, x_{\max_{n,h}}]$ by a fixed amount of fuel $dx_{i,h}$~; new nodes are only created when needed, so as to ensure that there will be a node in every non-void intervals $[x_{\min_{n,h}}, x_{\max_{n,h}}] \cap [x,x+dx_{i,h}]$, for all $x \in \mathbb{R}^+$ Step \ref{algostep:possibleStock} of algorithm \ref{algo:beam} depends on the type of the node. \begin{itemize} \item If $type_n=\mathtt{BO}$, we follow the imposed profile of constraint \texttt{CT6} with the authorized $\epsilon$-deviation, counted positively for $x_{min_{n,h}}$ and negatively for $x_{max_{n,h}}$. \item If $type_n=\mathtt{OUTAGE}$, we use \texttt{CT10} to compute stock variation during the outage, \texttt{CT7} (bounds on refueling) and \texttt{CT11} (maximum stock after refueling) to compute the interval after refueling~; when the production campaign starts, $x_{max}$ is computed with using the maximum allowed modulation of \texttt{CT12} and $x_{min}$ with using no modulation at all. \end{itemize} The node creation process is illustrated by figure \ref{fig:beam} in 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$ in graph $G_i$ inside the beam of possible stock intervals. Stars figure out \texttt{OUTAGE} nodes, circles are \texttt{BO} nodes. Note that no \texttt{BO} node is created at week $\mathtt{TA}_{k+1}$, because the paths that would go through such a node would break 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} Let us recall the goal of this graph~: we want to produce a production plan for powerplant $i$ with minimal reduced cost~: $\bar{c}(r,p,x) = \sum_k C_{i,k}\ .\ r(k) - C_{i,T+1}\ .\ x(T) - \sum_h \mu_h\ .\ p(h) - u_i$ where~: \begin{itemize} \item $C_{i,k}$ is the cost of a unit of fuel for campaign $k$ and powerplant $i$ and $C_{i,T+1}$ is the price for powerplant $i$ of a unit of fuel that can be sold at the end of the timeline~; \item $\mu_h$ is the dual variable associated to constraint \ref{covCMP}~; \item $x(T)$ is the quantity of remaining fuel at the end of the timeline~; \item $r(k)$ is the refuel quantity of campaign $k$~; \item $p(h)$ is the fuel spent at week $h$. \item $u_i$ is the dual variable associated to constraint \ref{1iCMP} \end{itemize} As we only consider production plans that can be \textit{supported} by some path $S$ in graph $G_i$, we can rewrite this expression in the following way~: $\bar{c}(r,p,x,S) = \sum_{a=(n,n') \in S} \left( C_{i,k_n}\ .\ r(a) - \delta_{type_{n'}=\mathtt{FINAL}}\ . \ C_{i,T+1}\ .\ x(T) - \sum_{h \in [h_{n},h_{n'}[} \mu_h\ .\ p(h) \right) - u_i$ where~: \begin{itemize} \item $S$ is the support of the production plan, \textit{i.e.} a set of arcs $a=(n,n')$ in graph $G_i$, with $n$ and $n'$ being the outgoing and incoming nodes of $a$; \item $r(a)$ is the refuel quantity at the beginning of $a$~; \item $\delta_{type_{n'}=\mathtt{FINAL}}$ is $1$ if $n'$ is the final node of $G_i$, $0$ else. \end{itemize} Let us call the expression in parenthesis the \textit{reduced cost $\bar{c}_a(r,p,x)$ of arc $a$}. This decomposition on the arcs of the graph means that if we weight each arc $a$ with its minimal reduced cost $\min_{r,p,x}(\bar{c}_a(r,p,x))$, we can easily find the production plan with minimal reduced cost by a standard minimal weight path search in the graph. The reduced cost of this production plan is obtained by substracting $u_i$ to the minimum weight path value. \paragraph{Weighting a \texttt{BO} arc} is straight forward because there is no refuel at the beginning of such an arc and the productions $p(h)$ are fixed for each timestep by constraint \texttt{CT6}, if we neglect the allowed deviation $\epsilon$. The weight of a \texttt{BO} arc $a$ is given by the following expression, where everything is known~: \begin{equation} w_a = \left( - \delta_{type_{n'}=\mathtt{FINAL}}\ . \ C_{i,T+1}\ .\ \left(x_n - \sum_{h \in [h_{n},h_{n'}[}p(h) \right) - \sum_{h \in [h_{n},h_{n'}[} \mu_h\ .\ p(h) \right)\label{eq:weightBOarc} \end{equation} \paragraph{Weighting a non final \texttt{OUTAGE} arc} is more intricate because the allowed modulation gives some latitude to the productions $p(h)$. Only are known the weeks of productions $(h_n,h_{n'})$ and the levels of fuel $(x_n,x_{n'})$ at the beginning and at the end of the arc. The problem can be stated as follows~:\\ Let $\bar{x}_n=\left( \frac{Q_{i,k_n} -1}{Q_{i,k_n}}\ .\ x_n - \mathtt{BO}_{i ,k_n-1} + \mathtt{BO}_{i ,k_n}\right)$ the quantity appearing in constraint \texttt{CT10},\\ and $\mathit{PMAX}_h = \sum_{t\in \mathit{timesteps(h)}}\ \mathit{PMAX}_i^t\ .\ D^t$ the maximal fuel spent during week $h$ in~: \begin{subeqnarray} \label{eq:weightOUTAGEarc} \emph{Find } & w_a = \min_{r,p,x} C_{k_n}\ .\ r(a) - \sum_{h \in [h_{n},h_{n'}[} \mu_h\ .\ p(h) \slabel{eq:waobj}\\ \emph{under constraints } & \left\{ \begin{array}[c]{lr} \mathbf{(1)}\ \bar{x}_n + r(a) = x_{n'} + \sum_{h \in [h_{n},h_{n'}[} p(h)\ \ & \mathtt{[CT10,CT9]}\\ \mathbf{(2)}\ \forall h \in [h_{n},h_{n'}[,\ \ 0 \leq p(h) \leq PMAX_h\ \ & \mathtt{[CT4,CT5]}\\ \mathbf{(3)}\ RMIN_{k_n} \leq r(a) \leq RMAX_{k_n} & \mathtt{[CT7]}\\ \mathbf{(4)}\ \bar{x}_n + r(a) \leq SMAX_{i,k_n} & \mathtt{[CT11]}\\ \mathbf{(5)}\ \sum_{h \in [h_{n},h_{n'}[} p(h) \geq -MMAX_{k_n} + \sum_{h \in [h_{n},h_{n'}[} PMAX_h & \mathtt{[CT12]}\\ \end{array}\right. \slabel{eq:wact} \end{subeqnarray} We can use \eqref{eq:wact}.(1) to rewrite $r(a)$ and remove it from the problem. It comes that we are looking for the minimum under constraints \eqref{eq:wact}.(2) to \eqref{eq:wact}.(5) of~: $w_a = C_{i,k_n}\left(x_{n'}-\bar{x}_n\right) + \sum_{h \in [h_{n},h_{n'}[} \left(C_{i,k_n}-\mu_h\right)\ .\ p(h)$ As the first part of the sum is \textit{known} for each arc, we only have to find the productions $p(h)$ that minimize $\sum_{h \in [h_{n},h_{n'}[} \left(C_{i,k_n}-\mu_h\right)\ .\ p(h)$. We tried two methods to achieve this goal~: \paragraph{The exact solution} is obviously given by the simplex method 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 breaking constraints (2) to (5) of \eqref{eq:wact}''~: \begin{algorithm}[H] \caption{Exact algorithm for minimizing arc weight} \label{algo:min_arc_weight_exact} \begin{algorithmic} \STATE set all $p(h)$ to $PMAX_h$ \STATE order the set of weeks by decreasing values of $\left(C_{i,k_n}-\mu_{h_m}\right)$~: $h_1, h_2 \dots $ \STATE start with $m=1$ \WHILE{constraints (2-5) of \eqref{eq:wact} are not satisfied} \STATE nullify $p(h_m)$ \STATE increment $m$ \ENDWHILE \WHILE{$\left(C_{i,k_n}-\mu_{h_m}\right)>0$ \AND constraints (2-5) of \eqref{eq:wact} can be satisfied} \STATE nullify $p(h_m)$ \STATE increment $m$ \ENDWHILE \end{algorithmic} \end{algorithm} However, we encountered two problems with the exact solution: \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 is not redhibitory though, because this can be countered by the \textit{diversification} strategy exposed in section \ref{sec:diversification}~; \item considering there are hundreds of thousands of arcs in the graph, the method is quite time consuming because it has to walk along the ordered set of weeks. \end{itemize} \paragraph{We used a heuristic} to give a good solution, even if not optimal: we chose to \textit{spread} the modulation on all weeks proportionally to $PMAX_h$, but with a different factor whether the \textit{reduced unit cost} $\bar{c}_{i,k,h}=C_{i,k}-\mu_h$ is positive or negative. For that purpose, we split the set of weeks into two sets: $\mathcal{H}_k^+ = \{h\in[0,H] \text{ s.t. } C_{i,k}-\mu_h>0\}$ and $\mathcal{H}_k^- = \{h\in[0,H] \text{ s.t. } C_{i,k}-\mu_h\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 $p(h)$ from these variables: $\begin{array}[c]{l} \forall h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^+,\ \ p(h)=\alpha_a\ .\ PMAX_h \\ \forall h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^-,\ \ p(h)=\beta_a\ .\ PMAX_h \\ \end{array}$ The problem becomes: find the factors $\alpha_a$ and $\beta_a$ that minimize the expression: $$\alpha_a\ .\ \left(\sum_{h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^+} (C_{i,k_n}-\mu_h)PMAX_h \right) + \beta_a\ .\ \left(\sum_{h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^-} (C_{i,k_n}-\mu_h)PMAX_h \right)$$ The parenthesized expressions do not have to be computed for each arc. Instead, we can compute them \textit{in advance} for each $(k,h,h')_{h