%\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} %%%%% 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}{LP_R} %%%%% Problem name \newcommand{\pld}{LP_s} % ----------------------------------------------------------------------------- \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} } \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 of the ROADEF/EURO 2010 Challenge. This is to our knowledge the first-ranked method among column generation methods proposed in the Challenge, and was ranked $4^{th}$ 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 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 by column generation, producing outage dates. 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 scenarii. Finally, the production quantities between two outages are optimized for each plant and each scenario by solving independent Linear Programs. \end{abstract} {\bf Keywords}: Energy, electricity production, outage, column generation, shortest path, linear programming. \section{Introduction} Optimization in the energy sector is a large subject. With the raise of raw materials prices and the growing awareness of environment concerns, energy consumption must be considered in a less greedy way in order to save resources. For years several researchers have worked on that topic as gathered in book \cite{NRJbook}. The increasing ratio of renewable energies in the global energy production enhances this trend. We report here a sample of works within that huge area of research. Bruglieri and Liberti in \cite{bioM} propose a model to improve the efficiency of biomass energy production plants. Mustakerov and Borissova in \cite{eol} give tools to determine the optimal deployment of a grid of wind turbines on a park knowing the wind direction. Finally several works have also been done on optimally planning energy production when several means can be used to produce 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 studied topics. Electricity companies often seek for performant optimization tools because, since total costs are very large, even small percentage gains through optimization often translate into large savings. Indeed in most countries energy can be produced through different ways: nuclear and other thermal energies provide 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. This subject has several difficulties, one of which is that the exact amount of the future demand and the exact price of raw materials are unknown. Moreover, there are many constraints on the energy production plants: for example, the level of water in dams has to stay above a certain threshold or maintenance operations have to be executed periodically on nuclear plants, periods during which they are not producing. 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 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. 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. The research project consists of 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. In medium-term electricity management, 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), which leads to consider multiple scenarios and ask 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 scnearios. In this paper a column generation method (see for example in \cite{CGbook} for a survey on this approach) is proposed to solve it. 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}) (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 (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 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. The method stops when no more negative reduced cost column is output from the pricing subproblem, the linear relaxation of the original master problem is solved then 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 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 approach developed in this paper is a column generation method, where only a subset of the whole set of possible columns are implicitely considered due to the time limit imposed by the competition. Several approximations have been done in the Pricing Subproblem and therefore the optimality of the final solution is not guaranteed. The main contributions of this paper are the following: original set partitioning-like reformulation of the problem, decomposition of the whole problem in two subproblems (master and slave), pricing subproblem solving based on shortest path computations in a particular state graph, since the slave problem is close to the Uncapacitated Lot Sizing problem, diversification method obtained by changing the objective function and adding different columns which may have a positive reduced cost. Moreover, the method is not limited to the idea of using columns generation but also uses rows generation. 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. This paper is organized as follow. The electricity production planning problem is presented in Section \ref{sec:problemstatement}. Section \ref{sec:approach} describes the overall method used to solve it. Section \ref{sec:solving_the_master_problem} gives the model of the master problem that is solved for fixing dates of outages and critical stock-level periods as first stage. Section \ref{sec:pricing} describes the pricing subproblem and its resolution method. Section \ref{sec:diversification} shows to improve the pricing subproblems in order to generate a greater number of columns assuring diversity. Section \ref{sec:ReloadProduction} exposes the second stage of the method for fixing reload and production quantities for each plant. Section \ref{sec:ComputationalResults} gathers the results obtained on the test instances. \section{Problem statement} \label{sec:problemstatement} The model proposed by the EDF company is the following. For the sake of consistency, the notations used in the ROADEF/EURO 2010 Challenge document \cite{sujetOnline} 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 introduced in the Challenge document are: \begin{itemize} \item $ha(i,k)$ : week of cycle $k$ during which nuclear plant $i$ goes offline (i.e. week of decoupling) \item $p(j,t,s)$ : production of non-nuclear plant $j$ during time step $t$ over scenario $s$ \item $p(i,t,s)$ : production of nuclear plant $i$ during 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} Variables $ha$ are integer, all other variables but $ec$ and $ea$ are continuous variables.\\ 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} \noindent 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$, $C_{i,T+1}$ is the selling price (or the gain) for the remaining fuel 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 have been used for writing the set partitioning formulation of the problem are the following: \begin{itemize} \item[(i)] Demand satisfaction: \[ \text{CT1} \quad \sum_j p(j,t,s) + \sum_i p(i,t,s) = DEM^{t,s} \quad \forall 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 during 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. Nevertheless, they have been used for improving the column generation process and they have been used for the row 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. % 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{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 \quad \forall s,t,i,k \] \noindent 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*} \forall s,t,i,k, &\text{ if }& t\in ec(i,k) \text{ and } x(i,t,s)h$ 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, 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 power plant $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 power plant 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 $[x_{\min_{n,h}}, x_{\max_{n,h}}]$ of possible stocks of fuel at week $h$ 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,h}}, x_{\max_{n,h}}] \cap [x,x+\Delta_i]$, for all $x \in \mathbb{R}^+$. Let us define the maximum production for power plant $i \in I$ over the week $h \in H$ as \[ EMAX_{h,i}= \sum_{t\in \mathit{timesteps(h)}} \mathit{PMAX}_i^t \ .\ D^t \] \noindent and \[ \overline{weekProd_i} = \frac{\sum_{c \in H} EMAX_{h,i}}{H} \] \noindent the mean maximum production per week for power plant $i \in I$, $\Delta_i$ is then computed as follows: \[ \Delta_i = \rho \ .\ \overline{weekProd_i} \] \noindent where $\rho$ is a value called the discretisation parameter that is chosen. 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} 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_{n,h}}$ is computed using the maximum allowed modulation of \texttt{CT12} and $x_{min_{n,h}}$ 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$ in graph $G_i$ inside the beam of possible stock intervals. Stars correspond to \texttt{OUTAGE} nodes, whereas 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 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}$ as explained by equation (\ref{subp}). Nevertheless, in order to point out that they are variables in this pricing subproblem, $e_{h\ell}$ becomes $e(h,l)$ and $\bar{c}_{\ell}$ becomes $\bar{c}_{\ell} (r,e,x)$ : \[ \bar{c}_{\ell}(r,e,x) = \sum_k C_{i,k} r(i,k) - C_{i,T+1} x(i,T) - \sum_h \mu_h e(h, \ell) - u_i \] \noindent where: \begin{itemize} \item $C_{i,k}$ is the cost of a unit of fuel for campaign $k$ and power plant $i$ and $C_{i,T+1}$ is the price for power plant $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{covLRMP}~; \item $x(i,T)$ is the quantity of remaining fuel at the end of the timeline~; \item $r(i,k)$ is the refuel quantity of campaign $k$~; \item $e(h, \ell)$ is the energy produced at week $h$~; \item $u_i$ is the dual variable associated to constraint \ref{1iLRMP} \end{itemize} As a production plan is represented by some path $\ell \subset A_i$ in graph $G_i$, we can rewrite this expression as: \[ \bar{c}_{\ell}(r,e,x) = \sum_{a=(n,n') \in {\ell}} \left( C_{i,k_n} r(a) - \delta_{type_{n'}=\mathtt{FINAL}} C_{i,T+1} x(i,T) - \sum_{h \in [h_{n},h_{n'}[} \mu_h e(h, \ell) \right) - u_i \] \noindent where: \begin{itemize} \item $r(a)$ is the refuel quantity at the beginning of arc $a$~; \item $\delta_{type_{n'}=\mathtt{FINAL}}$ is $1$ if $n'$ is the final node of $G_i$, $0$ otherwise. \end{itemize} Let us call the expression in parenthesis the \textit{reduced cost} $\bar{c}_a(r,e,x)$ of arc $a$. This decomposition on the arcs of the graph means that if we weigh each arc $a$ with its minimal reduced cost \[ w_a = \min_{r,e,x}(\bar{c}_a(r,e,x)) \] \noindent we can easily find the production plan with minimal reduced cost by solving a shortest path in the graph. The reduced cost of this production plan is obtained by substracting $u_i$ to 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 and the productions $e(h, \ell)$ 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'}[} e(h, \ell) \right) - \sum_{h \in [h_{n},h_{n'}[} \mu_h\ .\ e(h, \ell) \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(h, \ell)$. Only the weeks of production $(h_n,h_{n'})$, and the levels of fuel $(x_n,x_{n'})$ at the beginning and at the end of the arc, are known. 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 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 $h$ in: \begin{align} %\label{eq:weightOUTAGEarc} w_a = & \min_{r,e,x} C_{i,k_n} r(a) - \sum_{h \in [h_{n},h_{n'}[} \mu_h e(h, \ell) &&\label{eq:waobj}\\ \text{s.t.}\quad & \bar{x}_n + r(a) = x_{n'} + \sum_{h \in [h_{n},h_{n'}[} e(h, \ell) \quad & \mathtt{[CT10,CT9]} &\label{eq:wauno}\\ & 0 \leq e(h, \ell) \leq EMAX_h \quad \forall h \in [h_{n},h_{n'}[ \quad & \mathtt{[CT4,CT5]} &\label{eq:wadue}\\ & RMIN_{k_n} \leq r(a) \leq RMAX_{k_n} \quad & \mathtt{[CT7]} &\label{eq:watre}\\ & \bar{x}_n + r(a) \leq SMAX_{i,k_n} \quad & \mathtt{[CT11]} &\label{eq:waquattro}\\ & \sum_{h \in [h_{n},h_{n'}[} e(h, \ell) \geq -MMAX_{k_n} + \sum_{h \in [h_{n},h_{n'}[} EMAX_h \quad & \mathtt{[CT12]} &\label{eq:wacinque} \end{align} We can use (\ref{eq:wauno}) to rewrite $r(a)$ and remove it from the problem. It comes that we are looking for the minimum under constraints (\ref{eq:wadue}) to (\ref{eq:wacinque}) 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) e(h, \ell) \] As the first part of the sum is \textit{known} for each arc, we only have to find the productions $e(h, \ell)$ that minimize $\sum_{h \in [h_{n},h_{n'}[} \left(C_{i,k_n}-\mu_h\right) e(h, \ell)$. We tried two methods to achieve this goal: \paragraph{The exact solution} is obviously given by the simplex method, since (\ref{eq:wadue}) to (\ref{eq:wacinque}) is a linear model 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:wadue}) to (\ref{eq:wacinque})": \begin{algorithm}[H] \caption{Exact algorithm for minimizing arc weight} \label{algo:min_arc_weight_exact} \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_{h_m}\right)$~: $h_1, h_2 \dots $ \STATE start with $m=1$ \WHILE{constraints \ref{eq:wadue},\ref{eq:watre},\ref{eq:waquattro},\ref{eq:wacinque} are not satisfied} \STATE nullify $e_{h_m \ell}$ \STATE increment $m$ \ENDWHILE \WHILE{$\left(C_{i,k_n}-\mu_{h_m}\right)>0$ \AND constraints \ref{eq:wadue},\ref{eq:watre},\ref{eq:waquattro},\ref{eq:wacinque} can be satisfied} \STATE nullify $e_{h_m \ell}$ \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 can be overcome however 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 through the ordered set of weeks. \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 $EMAX_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 $e(h, \ell)$ from these variables: \begin{eqnarray*} &&\forall h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^+,\ \ e(h, \ell)=\alpha_a\ .\ EMAX_h \\ &&\forall h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^-,\ \ e(h, \ell)=\beta_a\ .\ EMAX_h \end{eqnarray*} 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)EMAX_h \right) + \beta_a \left(\sum_{h \in [h_n,h_{n'}[\ \cap\ \mathcal{H}_{k_n}^-} (C_{i,k_n}-\mu_h)EMAX_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