Crucial to freight carriers is the tactical planning of the service network. The aim is to obtain a cyclic plan over a given tactical planning horizon that satisfies predicted demand at a minimum cost. A central input to the planning process is the periodic demand, that is, the demand expected to repeat in every period in the planning horizon. We focus on large-scale tactical planning problems that require deterministic models for computational tractability. The problem of estimating periodic demand in this setting broadly present in practice has hitherto been overlooked in the literature. We address this gap by formally introducing the periodic demand estimation problem and propose a two-step methodology: Based on time series forecasts obtained in the first step, we propose, in the second step, to solve a multilevel mathematical programming formulation whose solution is a periodic demand estimate that minimizes fixed costs, and variable costs incurred by adapting the tactical plan at an operational level.

We report results in an extensive empirical study of a real large-scale application from the Canadian National Railway Company. We compare our periodic demand estimates to the approach commonly used in practice which simply consists in using the mean of the time series forecasts. The results clearly show the importance of the periodic demand estimation problem. Indeed, the planning costs exhibit an important variation over different periodic demand estimates and using an estimate different from the mean forecast can lead to substantial cost reductions. For example, the costs associated with the period demand estimates based on forecasts were comparable to, or even better than those obtained using the mean of actual demand.

Keywords  Freight transportation, tactical planning, large-scale, periodic demand, forecasting  demand.

§ INTRODUCTION

Freight transportation is essential to society and its economic development. In order to satisfy demand in a cost effective way, freight carriers are faced with a multitude of planning problems. In this context, Service Network Design (SND) is an important class of problems. Consider, for example, the Multicommodity Capacitated Fixed-charge Network Design (MCND) problem [Magnanti and Wong, 1984]. The objective is to design a capacitated network – a tactical plan – that allows to transport demand for a set of commodities between different origin-destination pairs at a minimum cost. The latter is given by the sum of fixed and variable costs. The tactical plan is defined over a given period (e.g., a week) and is repeated over a planning horizon (e.g., a few months). Given this cyclic nature of the tactical plan, it relies on an accurate representation of periodic demand.

In any realistic setting, demand for commodities is subject to uncertainty. This has naturally led to stochastic SND formulations <cit.>. As even deterministic SND problems are NP-hard, stochastic formulations are limited to fairly small size problems and cannot yet be used in most real large-scale applications. Hence, a wealth of practical applications rely on deterministic formulations and point estimates of periodic demand. In turn, periodic demand has an important impact on the resulting tactical plan and the associated costs. Despite its importance, there is no study in the literature focused on the periodic demand estimation problem linking time series forecasts to the tactical planning problem of interest. Our work addresses this gap and we use a MCND formulation for illustration purposes.

The impact of demand forecast errors on revenue has been studied in the context of airline revenue management. Through simulation analysis in a simplified setting, [Weatherford and Belobaba, 2002] show that reducing demand forecast errors by 25% increase revenue by a minimum of 1-2% which is a significant number. [Fiig et al., 2019] confirm those findings in a more complex airline revenue management setting and show that reduced forecast errors lead to increased revenue. As opposed to passenger transportation, freight carriers typically have flexibility regarding the routing of demand as long as it respects certain constraints, such as delivery time. Moreover, the freight demand origin-destination matrices are often unbalanced, meaning that there can be excess supply in certain directions. It implies that the cost associated with demand forecast errors can vary over commodities. This further motivates the importance of linking the periodic demand estimation and the corresponding SND problem.

Our proposed methodology proceeds in two steps: First we forecast demand for a given set of commodities for each period of the planning horizon. This corresponds to a multivariate multistep time series forecasting problem. Then, we define mappings from the time series forecasts to periodic demand. The different mappings lead to different periodic demand estimates. We solve a mathematical program explicitly linking these mappings and the MCND formulation. It selects the periodic demand that minimizes the fixed and variable costs over the planning horizon.

Brief Background on Time Series Forecasting. The periodic demand estimates are based on time series forecasts. There is an extensive related body of literature in statistics and machine learning. Our problem is particularly challenging for a number of reasons. First, the historical data on which the forecasting models rely, are results of operational decisions that are constrained by available capacity. Second, there are a large number of commodities and relatively long forecasting horizon which can lead to spatiotemporal correlations. Third, the demand varies over time and long-term dependencies and seasonality can be specific to each commodity. This highlights the potential need for modeling both commodity specific behavior and correlation between commodities while classic time series models and exponential smoothing methods assume independence across time series.

Freight demand forecasting works mainly focus on small networks of port terminals [Milenković et al., 2019] with either statistical models [Schulze and Prinz, 2009] or neural networks [Tsai and Huang, 2017]. With the recent development in intelligent transport systems and availability of large sources of data, forecasting methods have been shifting from model-based statistical models to data-driven machine learning approaches and more specifically deep learning models [Karlaftis and Vlahogianni, 2011]. Neural networks challenge the statistical models such as AR processes and Holt Winters method [Holt, 2004, Winters, 1960] with their augmented capacity to model non-linearities [Hornik et al., 1989]. The capacity of neural networks to model complex data to forecast traffic flows is increasingly exploited [Nguyen et al., 2018]. The Long Short-Term Memory (LSTM) recurrent neural network [Hochreiter and Schmidhuber, 1997, Sutskever et al., 2014] is a successful architecture to model both short-term and long-term dependencies [Längkvist et al., 2014].
Nevertheless, empirical evidence shows that it is still hard to achieve a level of accuracy comparable to that of classic time series models <cit.>.

Contributions. The paper offers both methodological and empirical contributions. First, we formally introduce the periodic demand estimation problem and propose a two-step methodology. Based on time series forecasts obtained in the first step, we propose a multilevel mathematical programming formulation whose solution is a periodic demand estimate that minimizes fixed and variable costs. The formulation hence explicitly links the periodic demand estimates to the tactical planning problem of interest. It is computationally tractable as it can be solved sequentially to optimality. Second, we describe a real large-scale application at the Canadian National Railway Company (CN). We present an extensive empirical study that clearly shows the importance of the periodic demand estimation problem. In this context, we compare different forecasting models from the statistics and deep learning literature. In turn we analyze the impact of the definition of periodic demand distinguishing between time series forecast errors and the errors introduced by different periodic demand estimates. Moreover, we benchmark against an approach used in practice that consists in averaging time series forecasts.

Paper Organization. The remainder paper is structured as follows. Next we formally introduce the periodic demand estimation problem. In Section <ref> we describe the proposed two-step methodology. We then focus on empirical results, first introducing our application in Section <ref>, followed by the results in Section <ref>. Finally, Section <ref> concludes and outlines some directions for future research.

§ PROBLEM DESCRIPTION

We start by briefly summarizing the planning process we consider: We take the point of view of a freight carrier that wishes to define a tactical plan. First, the carrier estimates the periodic demand for each commodity over the tactical planning horizon. Second, the periodic demand estimates are used as an input to solve the tactical planning problem of interest. The latter involves design decisions that are fixed over the tactical planning horizon, and flow decisions.
At the operational level, the tactical plan is adjusted – i.e., the flow decisions – according to actual demand realizations. In addition to these adjustments, other decisions could be taken to cope with demand fluctuations, such as outsourcing. There are hence two sources of costs to consider in the tactical planning process: the fixed cost of the tactical plan (design decisions) and the variable cost (flow and outsourcing decisions) resulting from the adjustments in each period.

We attend to large-scale problems that require a deterministic formulation to be tractable. Therefore, at the tactical planning level, the demand is treated as fixed and known while, in reality, it varies in each period. In this section we describe in detail the problem of estimating the periodic demand so as to minimize fixed and variable costs. We first introduce notation related to tactical and operational planning time horizons. We then define the various concepts of demand we encounter, followed by a description of an MCND formulation.
Finally, we describe the link between observed demand, periodic demand and the MCND formulation, which formally introduces our problem.

For the demand forecasting problem, it is important to distinguish the tactical and operational planning horizons. We therefore introduce two different notations related to time. First, the tactical planning horizon $\mathcal{T}$ can be divided into periods of equal length $t = 1, \ldots, T$. Second, each period $t$ can be further divided into $D$ time periods, $d = 1, \ldots,D$ that we here refer to as the operational horizon. Figure <ref> provides an illustration where the tactical horizon is composed of $T = 4$ weeks, and a week $t$ is composed of $D = 7$ days.

[thick, ->] (0.2,0) – (12,0) node [below] ;

in 1,...,4
(3*-2, 0.4) – node[pos=0.5] (point)  (3*-2, -0.4);

in 1,...,6
(3+3**0.14-2, 0.2) – node[pos=0.5] (pointD)  (3+3**0.14-2, -0.2);

in 1,...,6
(6+3**0.14-2, 0.2) – node[pos=0.5] (pointE)  (6+3**0.14-2, -0.2);

in 1,...,6
(9+3**0.14-2, 0.2) – node[pos=0.5] (pointF)  (9+3**0.14-2, -0.2);

(point1) node [above = 0.45cm of point1] $t=1$;
(point2) node [above = 0.45cm of point2] $t=2$;
(point3) node [above = 0.45cm of point3] $t=3$;
(point4) node [above = 0.45cm of point4] $t=T=4$;

(pointE1) node [ above = 0.10cm of pointE1] $d=1$;
(pointE6) node [ above = 0.10cm of pointE6] $d=6$;

(pointF2) node [black, above = 0.10cm of pointF2] ;
(pointF3) node [black, above = 0.10cm of pointF3] ;

(point2) edge[bend right = 35, <->, line width=0.25mm] node [below = 0.15cm of point1] Period (point3) ;

(point1) edge[bend right = 40, <->, line width=0.25mm] node [anchor=north] Tactical horizon $\mathcal{T}$ (point4) ;

(pointF2) edge[bend right = 80, <->, line width=0.10mm] node [anchor=north] Operational (pointF3) ;

Time scales for planning of a freight carrier

Let $\mathbf{y}_t$ be the demand vector of period $t$, $\mathbf{y}_t = ( y_{t1}, \dots, y_{tK})^{\top}$ where $y_{tk}$ is the quantity of commodity $k$ to be transported during period $t$. In this context, a commodity $k$ is characterized by its origin $o_k$, destination $d_k$ and type $\gamma_k$. We denote the set of commodities $\mathcal{K}$ and its cardinality $K$.
Let $\mathbf{y}_d^t$ be the demand vector for each operational time $d=1,\ldots,D$ within period $t$, $\mathbf{y}_d^t = ( y_{d1}^t, \dots, y^t_{dK})^{\top}$ where $y_{dk}^t$ is the demand for commodity $k$ to be carried at time $d$ in period $t$. The demand for a period $t$ is hence

\begin{equation}
y_{tk} = \sum_{d =1}^D y_{dk}^t, \quad k \in \mathcal{K}.
\label{eq:demand_relation}
\end{equation}

We introduce the demand matrix for horizon $\mathcal{T}$, $\mathbf{Y}^{\mathcal{T}} \in \mathbb{R}_+^{T \times K}$ with $[\mathbf{Y}^{\mathcal{T}}]_{tk} = y_{tk}$.
For a given tactical planning horizon $\mathcal{T}$, the plan is repeated at each $t = 1, \ldots, T$.

Let $\mathbf{y}^{\text{p}\mathcal{T}}$ be the periodic demand vector for tactical horizon $\mathcal{T}$,  $\mathbf{y}^{\text{p}\mathcal{T}} = ( y_1^{\text{p}\mathcal{T}}, \dots, y_K^{\text{p}\mathcal{T}})^{\top}$ where $y_k^{\text{p} \mathcal{T}}$ is the periodic demand for commodity $k$. To simplify the notation, we henceforth remove the superscript $\mathcal{T}$ but recall that the periodic demand and the demand matrix always refer to a given horizon.

We focus on estimating the periodic demand $\mathbf{y}^{\text{p}}$.
In this context it is important to note that time series forecasting models produce demand forecasts for each commodity in each period $t$. That is, at period $t_0$, the forecasting models output an estimate of $\mathbf{Y}$ denoted $\hat{\mathbf{Y}}$ which consists in $T$ point estimates $\hat{\mathbf{y}}_{t_0+1}, \dots, \hat{\mathbf{y}}_{t_0+T}$.
The periodic demands $y_k^{\text{p}}$ are then estimated from these forecasts $\hat{\mathbf{Y}}$, or, for validation or analysis, from $\mathbf{Y}$. Let $h$ denote the mapping of $\mathbf{Y}$ to a periodic demand vector $\mathbf{y}^{\text{p}}$:

\begin{equation}
\begin{split}
h \colon \quad \mathbb{R}_+^{T \times K} & \to \mathbb{R}_+^{K}  \\
\mathbf{Y} &\mapsto \mathbf{y}^{\text{p}} = h(\mathbf{Y}).
\end{split}
\label{eq:h_mapping}
\end{equation}

When the periodic demand is a mapping of the forecasts, we use the notation $\hat{\mathbf{y}}^{\text{p}} = h(\hat{\mathbf{Y}})$. Our objective is to define the mapping $h$ minimizing fixed and variable planning costs.

We provide an illustrative example in Figure <ref>. The two graphs show two demand distributions (shown with one dot per period) with identical mean (depicted with a solid line). The mean represents one particular mapping $\mathbf{y}^{\text{p}} = h(\mathbf{Y})$. In the example, we illustrate three other possible mappings: the maximum, median and third quartile. As opposed to the mean, these other mappings do not result in the same periodic demand estimates in the right and left-hand graphs.

Illustration of a periodic demand from point estimates for $T$ periods

We now introduce a path-based MCND formulation [Crainic, 2000] that we use for illustrating our methodology. An arc-based formulation can be found, e.g., in [Chouman et al., 2017]. Let $\mathcal{G} = (\mathcal{N}, \mathcal{A})$ denote a space-time graph where $\mathcal{N}$ is the set of nodes and $\mathcal{A}$ is the set of arcs.
Commodity $k$ uses a path $p$, i.e., a sequence of arcs in $\mathcal{G}$. The source node of the first arc is $o_k$, and the sink node of the last arc is $d_k$. Let $\mathcal{P}$ denote the set of paths. In the case of insufficient capacity, demand is outsourced and we denote $\mathcal{P}^{\text{out}}$ the paths corresponding to outsourcing options, such that $\mathcal{P}^{\text{out}} \subset \mathcal{P}$. Furthermore, let $\mathcal{P}_k$ denote the set of paths for commodity $k$, $\mathcal{P}_k^{\text{out}}$ the outsourcing paths for commodity $k$ such that $\mathcal{P}_k^{\text{out}} \subset \mathcal{P}_k$ and $\mathcal{K}_p$ the set of commodities that can use $p$. Note that here we refer to outsourcing in a broad sense. It could mean outsourcing to a third party, or making use of additional capacity from the same carrier that was not originally part of the plan. For example, in our intermodal rail transporation application (Section <ref>), outsourcing means using capacity from non-intermodal trains.

The MCND problem consists in satisfying demand at minimum cost. It has two categories of decision variables: Binary design variables $z_p, ~\forall p\in \mathcal{P},$ equal to one if path $p$ is used and zero otherwise, and flow variables $x_{pk} \geq 0, ~\forall k\in \mathcal{K}, p\in \mathcal{P}_k$. Depending on the type of freight (bulk versus containers, for instance), $x_{pk}$ is either continuous or integer. The path-based mixed integer linear programming formulation is:

MCND      min_z, x  ∑_p ∈𝒫 C_p^design z_p   +  ∑_k ∈𝒦 ∑_p ∈𝒫_k ∖𝒫_k^out C_p^flow x_pk +  ∑_k ∈𝒦 ∑_p ∈𝒫_k^out  C_p^out x_pk

s.t.  ∑_p ∈𝒫_k x_pk = y_k^p,  k ∈𝒦,

∑_k ∈𝒦_p x_pk ≤u_p z_p,  p ∈𝒫,

x_pk ≥0 ,  k ∈𝒦, p ∈𝒫_k,

z_p ∈{0,1},   p ∈𝒫.

The objective function (<ref>) includes a fixed design cost
$C_p^{\text{design}}\geq0$ for the paths built to transport demand. The second cost is the variable flow cost $C_p^{\text{flow}}\geq0$ which accounts for satisfied demand and the third cost $C_p^{\text{out}}\geq0$ is the flow cost of outsourced demand.
Constraints (<ref>) ensure that the periodic demand is satisfied for each commodity. Constraints (<ref>) enforce flows on selected paths only, and that the flow does not exceed the path capacity, $u_p$.

MCND is solved to obtain a tactical plan based on a given periodic demand. However, in practice, demand varies from one period to another. The tactical plan is therefore adjusted at the operational planning level. That is, in each period the commodity flows can be adjusted to satisfy the actual demand value of this period also taking into account other uncertain aspects, such as schedule delays. The observed data $\mathbf{y}_d^t$, typically used for training forecasting models, result from this operational planning process. Consequently, data at this level of detail can be constrained by the available services and
the observed demand may therefore not correspond to the true demand. This is known as censored data in the literature <cit.>.

In summary, we focus on estimating periodic demand $\mathbf{y}^{\text{p}}$ at a given time $t_0$ for a tactical planning horizon $\mathcal{T}$. The demand forecasts $\hat{\mathbf{y}}_{t_0+1},\ldots,\hat{\mathbf{y}}_{t_0+T}$ are obtained using historical data of demand $\{\mathbf{y}_d^s, s=t_0-1,t_0-2,\ldots,t_0-H, d=1,\ldots,D_s\}$ where $H$ is the number of periods in the historical data and $D_s$ is the number of operational time intervals in each period $s$.
The periodic demand should be defined such that it minimizes fixed costs, as well as variable costs associated with adapting the plan over the tactical planning horizon. It is hence necessary to link the mapping $h$ to the tactical planning problem of interest. In the following section, we propose a formulation for this purpose using MCND as an example of tactical planning problem formulation.

§ PERIODIC DEMAND ESTIMATION

Each time a tactical plan is to be defined, our approach proceeds in two steps. First, we use a time series forecasting model to predict demand for each period
in the planning horizon. Second, we solve a multilevel formulation for the joint periodic demand estimation and tactical planning problem. In the following subsection we describe assumptions and their implications on the time series forecasting problem. In Section <ref>, we delineate the mathematical programming formulation.

§.§ Time Series Forecasting

We use time series forecasting methods to predict, at period $t$, the $T$ point estimates $\hat{\mathbf{y}}_{t+1}, \dots, \hat{\mathbf{y}}_{t+T}$. As we highlight in the previous section, historical data captures operational flows which can be constrained by the supply and, therefore, may not correspond to actual demand. Time series forecasting with censored data is challenging and difficult to validate. In this section we introduce two weak assumptions that allow us to work with historical uncensored data with which we use time series forecasting methods.

Recall from the problem description that demand which cannot be satisfied by the planned capacity is outsourced. This is typically the case for carriers
as unsatisfied demand would otherwise accumulate over time periods. Taking into account the outsourcing, demand is hence assumed satisfied in each time period $t$. This leads us to the following assumption on the historical data.

Historical data is uncensored when aggregated over time periods. That is, $\mathbf{y}_{s}=\sum_{d=1}^{D_s} \mathbf{y}_d^s, s=t-1,t-2,\ldots,t-H$ are uncensored.

While this assumption simplifies the forecasting problem we note that the aggregation results in fewer data points to learn from.

The tactical planning problem formulation is based on a space-time graph. The departure and arrival times of a commodity $k$ are hence implicitly given by $o_k$ and $d_k$. We assume that the arrival and departure times are endogenous decisions, stated in other words in the following assumption.

Predicted demand per time period for each commodity, $\hat{\mathbf{y}}_t,~t=1,\ldots,T$, are sufficiently precise for tactical planning.

This is a weak assumption considering that, if exogenous predictions of $\hat{\mathbf{y}}_d^t,~d=1,\ldots,D_t$ are required for each $t=1,\ldots,T$, it is possible to define a model (different from the time series one) that projects $\hat{\mathbf{y}}_t$ down to that level.

§.§ A Multilevel Formulation

We define the feasible set of periodic demand vectors

\begin{equation}
\mathcal{Y} = \{ \hat{\mathbf{y}}^{\text{p}} = h_i(\hat{\mathbf{y}}_1, \ldots \hat{\mathbf{y}}_T), i=1,\ldots,I \}
\label{eq:periodic_set}
\end{equation}

by a finite set of mappings $h_i$, $~i=1,\ldots,I$  (<ref>). We propose the following multilevel formulation PDE for the periodic demand estimation problem.

PDE   min_𝐲̂^p       C^PDE = ∑_t=1^T [ ∑_p ∈𝒫 C_p^design z_p  +  ∑_k ∈𝒦 ∑_p ∈𝒫_k ∖𝒫_k^out C_p^flow x_tpk +  ∑_k ∈𝒦 ∑_p ∈𝒫_k^out  C_p^out x_tpk ]

s.t.      𝐲̂^p ∈𝒴

MCND      min_z, x ∑_p ∈𝒫 C_p^design z_p   +  ∑_k ∈𝒦 ∑_p ∈𝒫_k ∖𝒫_k^out C_p^flow x_pk +  ∑_k ∈𝒦 ∑_p ∈𝒫_k^out  C_p^out x_pk

s.t.  ∑_p ∈𝒫_k x_pk = ŷ_k^p,  k ∈𝒦,

∑_k ∈𝒦_p x_pk ≤u_p z_p,  p ∈𝒫,

x_pk ≥0 ,  k ∈𝒦, p ∈𝒫_k,

z_p ∈{0,1},   p ∈𝒫,

wMCND      min_x_1,…,x_T ∑_t=1^T [ ∑_k ∈𝒦 ∑_p ∈𝒫_k ∖𝒫_k^out C_p^flow x_tpk +  ∑_k ∈𝒦 ∑_p ∈𝒫_k^out  C_p^out x_tpk ]

s.t.  ∑_p ∈𝒫_k x_tpk = ŷ_tk,  t=1,…,T, k ∈𝒦,

∑_k ∈𝒦_p x_tpk ≤u_p z_p,  t=1,…,T, p ∈𝒫,

x_tpk ≥0 ,  t=1,…,T, k ∈𝒦, p ∈𝒫_k.

The upper level selects $\hat{\mathbf{y}}^{\text{p}}$ that minimizes the total fixed and variable costs over the whole tactical planning horizon. The objective function (<ref>) hence depends on the design and flow variables from the lower levels MCND and wMCND, respectively.

In wMCND we introduce the flow variables $x_{tpk}$ for commodity $k$ on path $p$ in period $t$ and determine flows for each period minimizing variable cost (<ref>) for a fixed design solution $z$ given by MCND. Constraints (<ref>) ensure that the demand is satisfied for each commodity in each period. The set of paths $\mathcal{P}_k$ includes outsourcing paths $\mathcal{P}_k^{\text{out}}$ for a commodity $k$, so constraints (<ref>) and (<ref>) can always be satisfied. Constraints (<ref>) enforce flows in each period to be only on selected paths and smaller than the capacity of the path. We draw the attention to the time series forecasts that occur in wMCND while the periodic demand estimates occur in MCND.

The decision variables of wMCND do not occur in the objective function (<ref>) of MCND. Thus, for $(z^*, x^*)$ an optimal solution of MCND, if $z^*$ is feasible for wMCND, then $z^*$ is an optimal solution to MCND-wMCND. We can therefore make the following claim.

If $z^*$ is feasible for wMCND, then MCND-wMCND can be solved sequentially to optimality for a fixed $\hat{\mathbf{y}}^{\text{p}}$.

In this work, we consider the four following mappings:

\begin{equation}
\mathbf{y}_{\text{max}}^{\text{p}} = h_1(\mathbf{y}_1, \ldots, \mathbf{y}_T)= \max_{t=1,\ldots,T} \mathbf{y}_t,
\label{eq:max_periodic}
\end{equation}

\begin{equation}
\mathbf{y}_{\text{mean}}^{\text{p}} =  h_2(\mathbf{y}_1, \ldots, \mathbf{y}_T) =  \frac{1}{T} \sum_{t=1}^T \mathbf{y}_t,
\label{eq:mean_periodic}
\end{equation}

\begin{equation}
\mathbf{y}_{\text{q2}}^{\text{p}} =  h_3(\mathbf{y}_1, \ldots, \mathbf{y}_T) =  Q_2(\mathbf{y}_{t}, t = 1 \ldots, T),
\label{eq:q2_periodic}
\end{equation}

\begin{equation}
\mathbf{y}_{\text{q3}}^{\text{p}} =  h_4(\mathbf{y}_1, \ldots, \mathbf{y}_T) =  Q_3(\mathbf{y}_{t}, t = 1 \ldots, T)
\label{eq:q3_periodic}
\end{equation}

which represent the maximum, mean, second quartile $Q_2$ and third quartile $Q_3$, respectively. The corresponding estimates from the forecasts are denoted
$\hat{\mathbf{y}}_{\text{max}}^{\text{p}}$, $\hat{\mathbf{y}}_{\text{mean}}^{\text{p}}$, $\hat{\mathbf{y}}_{\text{q2}}^{\text{p}}$ and $\hat{\mathbf{y}}_{\text{q3}}^{\text{p}}$. Given that we consider a discrete $\mathcal{Y}$ of small cardinality, we find the solution to PDE by solving MCND-wMCND for each $\mathbf{y}^{\text{p}}\in\mathcal{Y}$. In the following section, we describe a specific instance of (<ref>)-(<ref>) in the context of tactical planning for intermodal rail transportation.

§ APPLICATION

We illustrate our approach on the intermodal network of CN, composed of 24 main intermodal terminals and 133 origin-destination (OD) pairs. Figure <ref> depicts a map of the network. The railtracks extend from East to West of Canada and from Canada to South of the United States and gather 25 intermodal terminals. The railroad carries a variety of container types (20, 40, 45, 48 and 53-feet long), yet for tactical planning purposes they can be aggregated into either 40-feet or 53-feet containers. Indeed, 20-feet containers can be considered as half-40-feet containers. Other sizes (45, 48 and 53-feet) occupy the same locations on railcars <cit.> as 53-feet containers. A commodity is defined by an origin, a destination and a type of container, and we consider a total of $K = 170$ commodities. Two commodities can hence have the same OD pair but differ by the type of container. The tactical period is a week and a tactical horizon lasts $T = 10$ weeks. The train schedule is repeated each week and CN operates so that demand over a week is satisfied.  Assumption <ref> therefore holds. More precisely, in case of insufficient capacity on intermodal trains, they use general cargo trains or additional ad-hoc intermodal trains to satisfy demand.

Intermodal Network of the Canadian National Railway Company. Source: www.cn.ca

The specific MCND problem to the intermodal network of CN is the tactical block planning problem [Morganti et al., 2020]. A block refers to a consolidation of railcars. In this context, a set of railcars flowing as a single unit between a given OD pair and where containers loaded on the railcars have the same OD.
[Morganti et al., 2020] introduce a path-based Block Planning formulation (BP). It is defined using a space-time graph generated based on a schedule of intermodal trains. The graph contains 28,854 arcs and 15,269 nodes. A block is a path in this graph and the set is denoted $\mathcal{B}$ with $|\mathcal{B}|=2,208$. We keep this notation to be consistent with [Morganti et al., 2020] but note that the set $\mathcal{B}$ corresponds to the set of paths $\mathcal{P}$ in MCND. The set $\mathcal{B}$ contains a subset of artificial blocks $\mathcal{B}^{\text{artif}}$ whose role is to transport demand exceeding capacity. They hence correspond to the outsourcing paths. They are built without design cost, i.e., $C_b^{\text{design}} = 0, \forall b \in \mathcal{B}^{\text{artif}}$. Similarly to MCND, $\mathcal{B}_k$ and $\mathcal{B}_k^{\text{artif}}$ denote respectively the set of blocks and the set of artificial blocks for commodity $k$,  and $\mathcal{K}_b$ the set of commodities that can use $b$.

Below we briefly describe BP and refer to [Morganti et al., 2020] for more details. There are three categories of decision variables. First, the design variables $z_b, b \in \mathcal{B}$ where $z_b$ equals one if block $b$ is built, and zero otherwise. Second, integer flow variables $x_{bk},k\in\mathcal{K},b\in\mathcal{B}_k$ that equal the number of containers for commodity $k$ transported on block $b$. Third, auxiliary variables for the number of 40-feet $v_b^{40}$ and 53-feet $v_b^{53}$ double-stack platforms to carry the containers assigned to block $b\in\mathcal{B}$.

BP      min_x, z ∑_b ∈ℬ ∖ℬ^artif  C_b^design z_b +  ∑_k ∈𝒦 ∑_b ∈ℬ_k ∖ℬ_k^artif  C_bk^flow x_bk + ∑_k ∈𝒦 ∑_b ∈ℬ_k^artif  C_bk^out x_bk

s.t.  ∑_b ∈ℬ_k x_bk = y_k^p,    k ∈𝒦,

x_bk ≤y_k^p z_b,   k ∈𝒦, b ∈ℬ_k,

v_b^53 =  max[ 0, ⌈1/2 ( ∑_k ∈𝒦_b, τ_k = 53 x_bk - ∑_k ∈𝒦_b, τ_k=40  x_bk ) ⌉],

b ∈ℬ,

v_b^40 = ⌈1/2 ( ∑_k ∈𝒦_b x_bk  ) ⌉- v_b^53,   b ∈ℬ,

∑_b ∈ℬ_a ( L^40 v_b^40 + L^53 v_b^53 ) ≤u_a,   a ∈𝒜^TM,

z_b ∈{0,1},   b ∈ℬ,

v_b^40, v_b^53 ∈ℕ,   b ∈ℬ,

x_bk ∈ℕ,   k ∈𝒦, b ∈ℬ_k.

The objective function (<ref>) minimizes fixed and variable costs as well as a variable cost associated with outsourced demand (flow on artificial blocks). Constraints (<ref>) ensure that the demand is satisfied by either the network capacity or outsourcing. Constraints (<ref>) enforce flows to be on selected blocks only. Constraints (<ref>) and (<ref>) fix the number of platforms required to transport the demand. These constraints take into account how containers of different sizes can be double stacked. Since 40-feet platforms use less train capacity than 53-feet platforms, they are used whenever there are less 53-feet containers than 40-feet ones (40-feet container stacked in the bottom position and 53-feet container on top), and 53-feet platforms are used otherwise. Constraints (<ref>) ensure that the train capacity, expressed in number of feet, is not exceeded. The platform lengths are denoted $L^{40}$ and $L^{53}$, respectively. The train capacity, $u_a, a\in \mathcal{A}^{TM}$, is defined for the set of arcs in the space-time graph that represent moving trains, $\mathcal{A}^{TM}$. We denote by $\mathcal{B}_a$ the set of blocks that use train moving arc $a \in \mathcal{A}^{TM}$.

Finally, we give an order of magnitude of the size of the formulation. In the case of the instances we solve in this paper, there are over 386,000 variables and some 18,000 constraints.

§.§ Block Generation For Weekly Demand Inputs

While the tactical plan is computed for a weekly schedule, [Morganti et al., 2020] assume that the time at which the demand arrives to the system within the week is given exogenously. We provide an illustrative example in Figure <ref>. Demand enters the network via a node noted DIN which is associated to one admissible train departure node. Containers are either assigned to a block, or wait at the terminal, represented by flow on arcs called ContainersWaiting.

Under Assumption <ref>, we propose a slightly different block generation so that the model optimally distributes the weekly demand over the train departures.
For this purpose, we introduce a new set of nodes $\mathcal{N}^{\text{WIN}}$ such that there is one node WIN $n_{\theta}^{\text{WIN}} \in \mathcal{N}^{\text{WIN}}$ per terminal $\theta \in \Theta$, where $\Theta$ is the set of terminals in the network. We also introduce a new set of arcs
$\mathcal{A}^{\text{WIN}} = \{ (n_{\theta}^{\text{WIN}}, j) \mid \theta \in \Theta, j \in \mathcal{N}^{\text{DIN}}_{\theta}\}$, where $\mathcal{N}^{\text{DIN}}_{\theta}$ is the set of DIN nodes for terminal $\theta$.

We illustrate this change compared to [Morganti et al., 2020] in Figure <ref>. For each terminal, the node WIN receives the weekly demand input and splits the commodity flow to the different DIN nodes on the arcs $\mathcal{A}^{\text{WIN}}$ at no cost. Instead of arriving at different points in time, demand now arrives in one source node and the model selects the optimal distribution over the week.

[shape=circle,draw=black, fill = gray!100, scale = 0.8] (A) at (1,1) DIN;
[shape=circle,draw=black, draw = white, fill = white] (B) at (1,2.3) ;
[shape=circle,draw=black, draw = white, fill = white] (C) at (5, 2.3) ;
[shape=circle, draw=black, fill = gray!100, scale = 0.8] (D) at (5, 1) DIN';
[shape=circle,draw=black, draw = white, fill = white] (E) at (-0.1, -0.5) ;
[shape=circle,draw=black, draw = white, fill = white] (F) at (3.5, -0.5) ;
[densely dashed, ->] (A) edge node[left]  (B);
[->](A) edge node[below] ContainersWaiting (D);
[densely dashed, ->](D) edge node[left]  (C);
[->] (E) edge[ultra thick] node[right] Daily
Demand In
[->] (F) edge[ultra thick] node[right] Daily
Demand In

[Morganti et al., 2020]

[shape=circle,draw=black, fill = gray!100, scale = 0.8] (A) at (1,1) DIN;
[shape=circle,draw=black, draw=black, draw = white] (B) at (1,2.3) ;
[shape=circle,draw=black, draw=black, draw = white] (C) at (5, 2.3) ;
[shape=circle, draw=black, fill = gray!100, scale = 0.8] (D) at (5, 1) DIN';
[shape = circle, draw = black, fill = green!20, scale = 0.8] (G) at (1, -1)$n_{\theta}^{\text{WIN}}$;
[shape=circle,draw=black, draw = white, fill = white] (E) at (-0.1, -2.5) ;

[densely dashed, ->] (A) edge node[left]  (B);
[densely dashed, ->](D) edge node[left]  (C);
[->] (G) edge node[right] $a \in \mathcal{A}^{\text{WIN}}$ (A);
[->] (G) edge node[right] $a' \in \mathcal{A}^{\text{WIN}}$ (D);
[->] (E) edge[ultra thick] node[right] Weekly
Demand In

Our model

Illustration of the difference between [Morganti et al., 2020] and our model

§.§ Periodic Demand Estimation Problem

We present below the formulation PDE specific to our application. The MCND formulation is replaced by BP and we introduce a weekly BP formulation, wBP. For the latter, we introduce flow variables and auxiliary platform variables for each week $t$, $x_{tbk}, v_{tb}^{40}, v_{tb}^{53}, t \in \mathcal{T},  k \in \mathcal{K},  b \in \mathcal{B}$.

PDE BP min_𝐲̂^p      C^PDE = ∑_t=1^T [ ∑_b ∈ℬ ∖ℬ^artif  C_b^design z_b + ∑_k ∈𝒦 ∑_b ∈ℬ_k ∖ℬ_k^artif  C_tbk^flow x_tbk + ∑_k ∈𝒦 ∑_b ∈ℬ_k^artif  C_tbk^out x_tbk ]

s.t.      𝐲̂^p ∈𝒴,

BP      min_x, z ∑_b ∈ℬ ∖ℬ^artif  C_b^design z_b + ∑_k ∈𝒦 ∑_b ∈ℬ_k ∖ℬ_k^artif  C_bk^flow x_bk + ∑_k ∈𝒦 ∑_b ∈ℬ_k^artif  C_bk^out x_bk

s.t.  ∑_b ∈ℬ_k x_bk = ŷ_k^p,    k ∈𝒦,

x_bk ≤ŷ_k^p z_b,   k ∈𝒦,  b ∈ℬ_k,

v_b^53 =  max[ 0, ⌈1/2 ( ∑_k ∈𝒦_b, τ_k = 53 x_bk - ∑_k ∈𝒦_b, τ_k=40  x_bk ) ⌉],

b ∈ℬ,

v_b^40 = ⌈1/2 ( ∑_k ∈𝒦_b x_bk  ) ⌉- v_b^53,   b ∈ℬ,

∑_b ∈ℬ_a ( L^40 v_b^40 + L^53 v_b^53 ) ≤u_a,   a ∈𝒜^TM,

z_b ∈{0,1},   b ∈ℬ,

v_b^40, v_b^53 ∈ℕ,   b ∈ℬ,

x_bk  ∈ℕ,   k ∈𝒦, b ∈ℬ_k,

wBP      min_x_1, …, x_T ∑_t = 1^T ∑_k ∈𝒦 [ ∑_b ∈ℬ_k ∖ℬ_k^artif  C_tbk^flow x_tbk + ∑_b ∈ℬ_k^artif  C_tbk^out x_tbk ]

s.t.  ∑_b ∈ℬ_k x_tbk = ŷ_tk,    t ∈𝒯,  k ∈𝒦,

x_tbk ≤ŷ_tk z_b,   t ∈𝒯,  k ∈𝒦,  b ∈ℬ_k,

v_tb^53 =  max[ 0, ⌈1/2 ( ∑_k ∈𝒦_b, τ_k = 53 x_tbk - ∑_k ∈𝒦_b, τ_k=40  x_tbk ) ⌉],

t ∈𝒯,  b ∈ℬ,

v_tb^40 = ⌈1/2 ( ∑_k ∈𝒦_b x_tbk  ) ⌉- v_tb^53,   t ∈𝒯,   b ∈ℬ,

∑_b ∈ℬ_a ( L^40 v_tb^40 + L^53 v_tb^53 ) ≤u_a,   t ∈𝒯,  a ∈𝒜^TM,

v_tb^40, v_tb^53 ∈ℕ,   t ∈𝒯, b ∈ℬ ,

x_tbk ∈ℕ,   t ∈𝒯,  k ∈𝒦,  b ∈ℬ_k.

The objective function (<ref>) has the same structure as (<ref>), with design costs, flow costs and outsourcing costs.
We note that we define BP (<ref>)-(<ref>) using periodic demand $\hat{\mathbf{y}}^{\text{p}}$, while we define wBP (<ref>)-(<ref>) for fixed design variables and weekly demand $\hat{\mathbf{y}}_t$.
Furthermore, there is no fixed cost associated with $b \in \mathcal{B}^{\text{artif}}$. Therefore, a solution for BP is always feasible for wBP. Using Claim <ref>, we can solve BP-wBP sequentially.

§ COMPUTATIONAL RESULTS

Our dataset contains the observed daily container shipments of all types of containers on each origin-destination pair of the rail network, collected over 6 years, from December 2013 to November 2019. The tactical plan is weekly. To ensure the observed demand is not constrained by the supply, we do a weekly aggregation of 2,226 observations of daily demand for each commodity. This results in 318 observations per commodity. Hence, $y_{tk}$ refers to the number of containers of type $\tau_k$ to be carried from origin $o_k$ to destination $d_k$ during week $t$.

In this section, we first provide a descriptive analysis of the data. We report the results of the time series forecasting models in Section <ref>. Finally, in Section <ref> we report results for the periodic demand estimation problem. Note that, for confidentiality reasons, we only report relative numbers in all results.

§.§ Descriptive Analysis

The large size of the transportation network and the large number of commodities suggest that they can be spatiotemporal correlated. We provide here an analysis of the different types of correlations we identified: between commodities and between commodities and weather.

§.§.§ Correlation Between Commodities

We start by analyzing interweek correlation by computing, for each commodity, estimates of the Pearson correlation coefficient between weekly shipments over successive weeks. Figure <ref> presents the distribution of the coefficient over the commodities for lags from 1 to 10 weeks. It shows that there are substantial positive correlations between weeks for all commodities. Correlations are strong between successive weeks and are slightly weaker for longer lags.

We now turn our attention to intercommodity correlations.
For each pair of commodities, we compute estimates of the Pearson correlation coefficient between weekly shipments over successive weeks. We present in Figure <ref> the distribution of those correlations for lags from 1 to 10 weeks. There are 170 commodities, hence 28,730 intercommodity correlation coefficients to examine at each lag. Correlations vary across pairs of commodities. Some are large (positive or negative) for short and longer lags, while 50% of the pairs have a weak correlation between -0.2 and 0.2. Outliers at each time lag represent 9% of the significant coefficients.

[t]@l@ Distribution of interweek correlations
over commodities.

[t]@l@ Distribution of intercommodity correlations
over commodities.

Distribution of the Pearson correlation coefficient for the two different types of correlations. We consider only the coefficients which met the 95% confidence level threshold and we indicate their proportion in parenthesis.

In summary, the data shows evidence of two types of correlation: strong positive interweek correlation between the weekly shipments over successive weeks, and strong intercommodity correlation between the weekly shipments of different commodities. Interweek correlation highlights the potential need for an autoregressive model while intercommodity correlation highlights a potential need for a model able to learn various dependence structures from the data.

§.§.§ Correlation Between Demand and Weather

Weather can be an important aspect for rail freight transportation, especially in North America where railways extend over the subcontinent which is subject to major weather disruptions such as snowstorms. To assess the importance of weather on shipments, we use meteorological data and estimate the Pearson correlation coefficient between weekly shipments and several weather indicators. The data comes from National Oceanic and Atmospheric Administration (NOAA, https://www.ncdc.noaa.gov/) for terminals in the United States and from Statistics Canada (https://www.statcan.gc.ca) for terminals in Canada. More precisely, we use the average daily temperature and total daily snowfall in centimeters for the main 17
terminals for the complete time range covered by the data. We compute the weekly temperature as the average temperature over the week and the accumulated snow (cm) as the sum of the daily values.
For each terminal, at each week $t$, we sum the total departing and arriving demand over all commodities.

Figure <ref> shows the distribution over the terminals of the Pearson correlation coefficient between accumulated snow over week $t$ and departing and arriving demand at week $t+\text{lag}$. We note a negative correlation for both arriving and departing demand for successive weeks. It is weaker for longer lags.

[t]@l@ Departing Demand

[t]@l@ Arriving Demand

Distribution of the Pearson correlation coefficient between accumulated snow over week $t$ and demand (arriving and departing) at week $t+\text{lag}$ over all the terminals. We consider only the coefficients which met the 95% confidence level threshold and we indicate their proportion in parenthesis.

We present in Figure <ref> the distribution over the terminals of the Pearson correlation coefficient between the average temperature and the demand arriving or departing over terminals. It shows an average positive correlation between demand and temperature which is weaker for longer lags. Some terminals have, however, negative correlations. This can be explained by a seasonality effect. On the one hand, summer and spring are busier periods for most ODs and temperature are higher than the rest of the year. On the other hand, for import terminals, fall and the Chinese New Year are busier periods when temperatures are lower.

[t]@l@ Departing Demand

[t]@l@ Arriving Demand

Distribution of the Pearson correlation coefficient between average temperature at week $t$ and demand (arriving and departing) at week $t+\text{lag}$ over all the terminals. We consider only the coefficients which met the 95% confidence level threshold and we indicate their proportion in parenthesis.

Correlation between demand and weather highlights the potential need for a forecasting model which can include weather features. In the following section we present forecasting results for different forecasting models, ranging from very simple ones based on strong independence assumptions to neural networks that relax some of those assumptions.

§.§ Time Series Forecasting Results

We divide the dataset into a training, validation and test sets. Time series forecasting models require the last seen observed data before predicting the demand to come. Thus, we use the first 5 years of data for training (December 2013 - December 2018), the next 4 months of data for validation (January 2019 - April 2019) and the last 7 months (May 2019 - November 2019) for testing. At each week $t_0$ in the dataset, we forecast demand for all commodities for $t = t_0+1, \dots ,t_0+T$. We consider $T=10$ weeks. To simplify the notation, we assume $t_0 = 0$ and refer to the estimates by $t = 1, \dots, T$.

We compare four types of models detailed below: a simple model that uses last observed values as prediction (CONSTANT) as well as autoregressive (AR), feedforward neural network (FFNN) and recurrent neural network (RNN) models.

CONSTANT This is the simplest possible model. It is based on the assumptions that commodities $k \in \mathcal{K}$ are independent and that the demand from observed week $t_0$ is the forecasts for the next $T$ weeks:

\begin{equation}
\label{eq:naive}
\hat{\mathbf{y}}_1 = \dots = \hat{\mathbf{y}}_T = \mathbf{y}_0.
\end{equation}

AR For each commodity, we fit an autoregressive $AR(p)$ process on the training data. This implies that the commodities are treated as independent. We use the estimated coefficients $\hat{\bm{\phi}} = (\hat{\phi}_1, \dots, \hat{\phi}_p)^T$ to compute the forecasts on the test set. The multistep forecasts for $t > 1$ are obtained using

\begin{equation}
\hat{y}_{tk} = \hat{\phi}_{1k} \hat{y}_{t-1, k} + \dots + \hat{\phi}_{tk} y_{0k} + \dots + \hat{\phi}_{pk} y_{t-p, k}.
\label{eq:multistep_ar}
\end{equation}

FFNN and RNN We leverage the capacity of these models to forecast demand for multiple commodities simultaneously. As opposed to CONSTANT and AR, they relax the independence assumption on the commodities. We build one main neural network architecture described in Figure <ref>. It is composed of 2 inputs layers, a stack of hidden dense or recurrent layers and one output layer. We consider several variants of this architecture. When the layers in the grey square are dense (or feed-forward), it forms the feed-forward architecture (FFNN). When the layers are LSTM, it forms the recurrent architecture (RNN). The dimension of the output layer is equal to the number of commodities that we predict simultaneously.

node distance = 0mm and 5mm,
block/.style = draw, rounded corners, minimum height=5em, text width=5em, align=center,
block2/.style = draw, rounded corners, minimum height=5em, text width=4em, align=center,
halfblock/.style = block, minimum height=2em,
halfblock2/.style = block, minimum height=2em, fill=black!10,
line/.style = draw, -Latex

(n1) [block]  HIDDEN LAYERS;
(n2) [block2, right=of n1] OUTPUT
LAYER;  library

(n1a) [halfblock, below left=of n1.north west]  LINEAR LAYER;
(n1b) [halfblock2, above left=of n1.south west]  ;

(n1a1) [halfblock, left=of n1a]  INPUT LAYER 1;
(n1b1) [halfblock, left=of n1b]  INPUT LAYER 2;

[line]    (n1)  edge (n2)
(n1a) edge (n1.west |- n1a)
(n1b)  –  (n1.west |- n1b)
(n1a1) edge (n1a)
(n1b1) edge (n1b);

Neural network architecture

At each period $t_0$, the output layer computes the forecast for the next time step $t_0+1$. Input Layer 1 is dedicated to external data such as weather features to model the correlation between demand and weather. Input Layer 2 is dedicated to the autoregressive modeling. It takes as input the demand of the previous weeks of all commodities, either observed or forecasted. When we generate the multistep demand forecasts, for Input Layer 2, we use $\hat{\mathbf{y}}_{t_0+1}$ to forecast demand at $t_0+2$, $\hat{\mathbf{y}}_{t_0+1}$ and $\hat{\mathbf{y}}_{t_0+2}$ for $t_0+3$, until $t_0+T$ with $\hat{\mathbf{y}}_{t_0+1}, \ldots, \hat{\mathbf{y}}_{t_0+T-1}$.

We evaluate several variants of the architecture to model the spatiotemporal correlations and the correlation with weather. They differ depending on the considered input data: lagged observed/forecasted demand and/or weather features. To train, validate and test the model, we use real observed weather data described in Section <ref>. At prediction time, such information is not available and we would then rely on weather forecasts. These results are hence designed to assess the potential for weather related features in an optimistic setting in regards to their accuracy.

Forecasting 170 commodities simultaneously requires a rich and large dataset. Our dataset is fairly limited as it contains only 318 weekly demand data points for each commodity. To facilitate the forecasting task, we create a partition of $\mathcal{K}$ and train a neural network for each set of the partition. For this purpose, we split the set of commodities into 2 subsets: $\mathcal{K}_{53}$ which contains commodities of 53-feet container type and $\mathcal{K}_{40}$ which contains commodities of 40-feet container ($ |\mathcal{K}_{53}| = 58$ and $|\mathcal{K}_{40}| = 112$).
Table <ref> summarizes the model variants.
Their name include a letter “W” if weather features are used to train the model, and “SPLIT” to indicate a partition of $\mathcal{K}$.

Commodities     Weather features     Autoregressive Features
RNN               $\mathcal{K}$

FFNN               $\mathcal{K}$

RNN-W             $\mathcal{K}$

FFNN-W            $\mathcal{K}$

RNN-W-SPLIT1      $\mathcal{K}_{53}$

RNN-W-SPLIT2      $\mathcal{K}_{40}$

FFNN-W-SPLIT1     $\mathcal{K}_{53}$

FFNN-W-SPLIT2     $\mathcal{K}_{40}$
Features and set of commodities for each variant of the neural network architecture evaluated

The neural networks are trained with the backpropagation algorithm and the stochastic gradient descent using a Mean Squared Error (MSE) loss. For each model, we do a hyperparameter search with the Tree of Parzen Estimators implemented in the python library Hyperopt [Bergstra et al., 2013]. We select the set of hyperparameters which minimizes the MSE on the validation dataset. We report the detailed input features and the chosen set of hyperparameters for each trained model in the Appendix (Table <ref>).

§.§.§ Forecast Accuracy Measures

We measure the accuracy on the test set with two metrics: the Weighted Absolute Percentage Error (WAPE)

\begin{equation}
\text{WAPE}_k = \frac{\sum_{t_0 \in \mathcal{D}_{\text{test}}} \sum_{t = 1}^T |y_{t_0+t,k} - \hat{y}_{t_0+t,k}|}{\sum_{t_0 \in \mathcal{D}_{\text{test}}} \sum_{t = 1}^T y_{t_0+t,k}} \times 100, ~k=1,\ldots,K,
\label{eq:wape}
\end{equation}

where $\mathcal{D}_{\text{test}}$ is the test set of size $N$, and the Root Mean Squared Error (RMSE)

\begin{equation}
\text{RMSE}_k = \sqrt{ \frac{1}{N \times T} \sum_{t_0 \in \mathcal{D}_{\text{test}}} \sum_{t=1}^T \left( y_{t_0+t,k} - \hat{y}_{t_0+t,k} \right)^2 }, ~k=1,\ldots,K.
\label{eq:rmse}
\end{equation}

The WAPE is a weighted version of the Mean Absolute Percentage Error (MAPE) which handles small or zero demand values. Low demands are frequent in our data for some commodities with sparse demand over the year. Hence the importance of having a metric independent to the scale of the time series such as WAPE. The RMSE puts a high weight on large errors, which is also an important metric to consider.

§.§.§ Results

Table <ref> reports the performance metrics averaged over all commodities. We note that the metrics for models based on a partition of the commodities (SPLIT) are averaged over the partitions.
The results show that the AR has the best performance and it is considerably better than the NN models. The CONSTANT baseline has a WAPE close to AR but a considerably worse RMSE. The descriptive statistics in Section <ref> show strong intercommodity correlations. Nevertheless, the AR model, based on the assumption that demands for commodities are independent, performs better than the NN models where this assumption is relaxed. Neural networks have more parameters to fit on the same limited data. While they have the capacity to model non-linear relationships, they also require more data to be trained. We believe that our data containing only 318 observations for each commodity is too limited for training the NN models. Moreover, we note that these findings are consistent with other studies <cit.>. That is, basic time series models outperform deep learning models on difficult time series forecasting tasks, such as this one.

Model     RMSE     WAPE
CONSTANT     86.0     34.7%

AR     78.0     34.0%

FFNN     105.2     38.7%

FFNN-W     84.8     37.1%

FFNN-W-SPLIT     90.4     37.2%

RNN     105.1     37.8%

RNN-W     86.3     37.8%

RNN-W-SPLIT     85.2     38.6%

Performance metrics of the forecasting models averaged over all commodities. The best and second best metric values are highlighted in bold.

The results for the deep learning models confirm that adding weather features and considering all commodities simultaneously help to improve the performance. This is consistent with the descriptive statistics reported in Section <ref>. Previous demands for all commodities constitute relevant information to consider.

Henceforth, we keep two forecasting models when analyzing periodic demand results: the overall best performing model (AR) as well as the best deep learning model (FFNN-W).

§.§ Periodic Demand Estimation

We divide the results related to periodic demand estimation into two parts. The purpose is to disentangle the errors associated with the periodic demand estimation from those associated with the demand forecasts. The two parts are briefly described in the following:

* Analysis 1: the impact of periodic demand estimation. We assume we have no forecast errors, i.e., the carrier knows perfectly the demand to come for the planning horizon. The periodic demand is estimated with the mappings (<ref>)-(<ref>) from historical data (ground truth values) and we compute the tactical costs generated by those periodic demands.

* Analysis 2: the impact of imperfect demand forecasts. We estimate the periodic demands from the forecasts obtained with the AR and FFNN-W models. We compute the associated tactical costs and compare them to Analysis 1.

We consider two demand instances – $I1$ and $I2$ – from the test set of the forecasting models. They are from two distinct periods: $I1$ corresponds to end of spring and beginning of summer, from May 6th to July 14th, 2019 and I2 corresponds to end of summer and beginning of fall, from July 29th to October 6th, 2019. Both instances have $K=170$ commodities and a tactical planning horizon of $T=10$ weeks. Instance $I2$ corresponds to a busier period for the carrier: Figure <ref> shows the difference of the total demand summed over all commodities of $I2$ at each week relative to the total demand of $I1$. We note that the total demand for $I2$ can be up to 20% higher than its $I1$ counterpart.

All the results are generated by fixing
the variable $\hat{\mathbf{y}}^{\text{p}}$ in PDE and sequentially solving BP-wBP. We recall that the periodic demand $\mathbf{y}^{\text{p}}$ is a mapping from the real demand values,  $\hat{\mathbf{y}}^{\text{p}}$ is a mapping from the demand forecasts, $\mathbf{Y}$ is the matrix of real demand values and $\hat{\mathbf{Y}}$ is the matrix of demand forecasts. We denote by  $\textbf{BP}(\mathbf{y}^\text{p})$ and $\textbf{BP}(\hat{\mathbf{y}}^\text{p})$ when BP is solved with $\mathbf{y}^\text{p}$ and $\hat{\mathbf{y}}^\text{p}$, respectively, in constraints (<ref>). Furthermore, we denote by $\textbf{wBP}(\mathbf{Y})$ and $\textbf{wBP}(\hat{\mathbf{Y}})$ when wBP is solved with demand values $\mathbf{Y}$ and $\hat{\mathbf{Y}}$, respectively, in constraints (<ref>).

Percentage difference of the total demand at each week of $I2$ relative to  $I1$

§.§.§ Analysis 1: The Impact of Periodic Demand Estimation

The analysis in this section is based on two sets of results for demand instances $I1$ and $I2$:

* Reference: We solve $\textbf{BP}(\mathbf{y}_t)$ for $t=1,\ldots, T$. In other words, we do not restrict the demand to be periodic and we obtain a tactical cost $C_{\text{ref}}^{\text{PDE}}$ (<ref>).

* Periodic: We solve  $\textbf{BP}(\mathbf{y}_{.}^{\text{p}})-\textbf{wBP}(\mathbf{Y})$ with the four periodic demand vectors $\mathbf{y}^{\text{p}}_{\text{mean}}$, $\mathbf{y}^{\text{p}}_{\text{max}}$, $\mathbf{y}^{\text{p}}_{\text{q2}}$ and $\mathbf{y}^{\text{p}}_{\text{q3}}$. For each solution, we compute the percentage gap to the reference cost $C_{\text{ref}}^{\text{PDE}}$.

Table <ref> reports the results where each number is the percentage gap to the reference cost. The cost $C^{\text{PDE}}$ is a non-trivial trade-off between the three components $C^{\text{design}}$, $C^{\text{flow}}$ and $C^{\text{out}}$, that are linked with demand. Hence we cannot analyze each component independently of the others. The parameter values in the objective function (<ref>) have been chosen with CN to best represent their costs. We note that $C_{\text{ref}}^{\text{PDE}}$ (<ref>) constitutes a lower bound: at each week, the blocks built are specific for the demand to best exploit the network capacity.

Definition $\mathbf{y}^{\text{p}}$
in $\textbf{BP}$      4cPercentage difference of costs

1c     $C^{\text{PDE}}$

$\mathbf{y}^{\text{p}}_{\text{max}}$      77.3%     -2.8%       -3.4%     84.0%

$I1$     $\mathbf{y}^{\text{p}}_{\text{mean}}$     93.3%     -16.9%      -4.4%     101.3%

$\mathbf{y}^{\text{p}}_{\text{q2}}$       93.8%     -16.8%      -3.3%     101.8%

$\mathbf{y}^{\text{p}}_{\text{q3}}$       28.4%     -7.5%       -1.2%     30.8%

$\mathbf{y}^{\text{p}}_{\text{max}}$      73.0%     -14.1%     -5.7%     76.6%

$I2$     $\mathbf{y}^{\text{p}}_{\text{mean}}$     40.6%     -20.0%     -2.8%     42.6%

$\mathbf{y}^{\text{p}}_{\text{q2}}$       41.6%     -23.2%     -2.8%     43.6%

$\mathbf{y}^{\text{p}}_{\text{q3}}$       23.1%     -14.8%     -2.6%     24.3%

Tactical costs relative to the reference. We use historical data for the demand at each week, that is $\textbf{wBP}(\mathbf{Y})$. The minimum value for $C^{\text{PDE}}$ indicates the optimal periodic demand estimate in the subset $\mathcal{Y}$ we consider.

Two important findings emerge. First, the tactical cost has an important variation over the different periodic demand estimates. This underlines the importance of the periodic demand estimation problem. Second, using another estimate than the commonly used mean periodic demand $\mathbf{y}^{\text{p}}_{\text{mean}}$ can lead to an important cost reduction. In our case using the third quartile $\mathbf{y}^{\text{p}}_{\text{q3}}$ reduces the total costs by 33.6% in $I1$, and by 12.4% in $I2$. When a smaller periodic demand estimate is used, $\mathbf{y}^{\text{p}}_{\text{mean}}$ and $\mathbf{y}^{\text{p}}_{\text{q2}}$ for instance, less blocks are built, at the expense of outsourced demand.

We analyze the total periodic demand relative to the total actual demand to explain the results. We report the values in  Table <ref>. The total periodic demand of $\mathbf{y}^\text{p}_{\text{max}}$ is almost three times larger than the one of $\mathbf{y}^{\text{p}}_{\text{q3}}$ in both instances which explains the large gap in their costs. The design made from solving $\textbf{BP}(\textbf{y}_{\text{q3}}^{\text{p}})$ is based on a lower estimation of the periodic demand than $\textbf{BP}(\textbf{y}_{\text{max}}^{\text{p}})$ which generates fewer blocks in the network. The advantage
is that they are better used: the design from solving $\textbf{BP}(\textbf{y}_{\text{max}}^{\text{p}})$ is made of a large number of small blocks which cannot accommodate all demands. Hence the increase in outsourced demand. We highlight that while the increase in outsourced demand seems large, it represents only few percent of the total demand.

$I1$     $I2$
$\sum_{t=1}^T\sum_{k \in \mathcal{K}} y_{tk}$     -     -

$\mathbf{y}^\text{p}_{\text{max}}$     34.8%     31.8%

$\mathbf{y}^{\text{p}}_{\text{mean}}$     0.0%     0.0%

$\mathbf{y}^\text{p}_{\text{q2}}$     -1.0%     -0.5%

$\mathbf{y}^{\text{p}}_{\text{q3}}$     11.8%     11.7%

Total periodic demand. The point of reference is indicated by a dash.

§.§.§ Analysis 2: The Impact of Imperfect Demand Forecasts

In practice, carriers rely on demand forecasts to estimate a periodic demand and build the design for the tactical planning horizon. Then, at each week of the horizon, they adapt operationally the tactical plan based on observed demand. In this section, we follow our methodology and analyze the quality of the solution a posteriori.

We proceed in two steps: First, we estimate the periodic demand by solving $\textbf{BP}(\hat{\mathbf{y}}^{\text{p}})-\textbf{wBP}(\hat{\mathbf{Y}})$ for each mapping $h$ and select the one minimizing $C^{\text{PDE}}$ (<ref>). Second, we assess the tactical cost associated with the periodic demand estimate. For this purpose we solve $\textbf{BP}(\hat{\mathbf{y}}^{\text{p}})-\textbf{wBP}(\mathbf{Y})$. Note that we use real demand values $\mathbf{Y}$ to accurately assess this cost. Hence, it is affected by two combined sources of error: the one of the periodic demand estimation and the forecast error, both discussed separately in previous sections.

Step 1: Periodic demand estimation We report results in Table <ref> and indicate by a dash the point of reference ($\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$). We note that the value of $\hat{\mathbf{Y}}$ depends on both the forecasting model and the instance. In other words, we can only compare results from the same forecasting model and the same instance. The results show that, for both forecasting models and both instances, the tactical costs are minimized with $\hat{\mathbf{y}}_{\text{max}}^{\text{p}}$. In Analysis 1 with historical data, tactical costs were minimized with ${\mathbf{y}}_{\text{q3}}^{\text{p}}$. This is because forecasting models smooth demand and struggle to forecast accurately the peaks. With historical data, we have access to the maximum of demand which can be an outlier, thus expensive.
Following the methodology, we choose the periodic demand $\hat{\mathbf{y}}_{\text{max}}^{\text{p}}$ for Step 2 for both forecasting models.

2cForecasting Model

Demand $\mathbf{y}^{\text{p}}$
in BP           1c|AR     FFNN-W

$\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$         1c|-38.4%     -32.2%

$I1$     $\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$          1c|-      -

$\hat{\mathbf{y}}^{\text{p}}_{\text{q2}}$          1c|17.5%     28.0%

$\hat{\mathbf{y}}^{\text{p}}_{\text{q3}}$          1c|-30.8%      -18.8%

$\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$          1c|-7.4%     -28.4%

$I2$     $\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$           1c|-     -

$\hat{\mathbf{y}}^{\text{p}}_{\text{q2}}$           1c|12.7%     29.2%

$\hat{\mathbf{y}}^{\text{p}}_{\text{q3}}$           1c|8.8%     -11.7%

Percentage difference of tactical cost $C^{\text{PDE}}$ resulting from solving $\textbf{BP}(\hat{\mathbf{y}}^{\text{p}})-\textbf{wBP}(\hat{\mathbf{Y}})$ with forecasts of demand from two models: AR and FFNN-W. For each model and each instance, we compare the value with the one from $\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$.

Step 2: Assessment of the tactical costs Table <ref> reports the relative costs resulting from solving $\textbf{BP}(\hat{\mathbf{y}}^{\text{p}})-\textbf{wBP}(\mathbf{Y})$. We use the same reference as in Analysis 1, that is, the lower bound on $C^{\text{PDE}}$. In addition to the best periodic demand identified by our methodology, we report the result for $\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$. Consistent with the findings in Analysis 1, we note that using the latter leads to a large increase in costs. Despite the relatively large forecast errors reported in Table <ref>, the results show that estimating periodic demand using our methodology can even reduce the costs compared to $\mathbf{y}^{\text{p}}_{\text{mean}}$ computed on perfect information (historical data).

Demand $\mathbf{y}^{\text{p}}$
in $\textbf{BP}$      4cPercentage difference of costs

1c     $C^{\text{PDE}}$

Historical data        $\mathbf{y}^{\text{p}}_{\text{mean}}$      93.3%      -16.9%     -4.4%     101.3%

Historical data     	$\mathbf{y}^{\text{p}}_{\text{q3}}$	      28.4%      -7.5 %     -1.2%     30.8%

$I1$     AR                   	$\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$       55.9%      -17.8%     -1.9%     60.6%

AR                   	$\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$      125.1%     -14.1%     -3.9%     135.8%

FFNN-W               	$\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$       48.2%      -12.3%     -2.5%     52.4%

FFNN-W               	$\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$      117.4%     -16.2%     -6.0%     127.5%

Historical data     	$\mathbf{y}^{\text{p}}_{\text{mean}}$	         40.6%      -20.0%     -2.8 %     42.6%

Historical data     	$\mathbf{y}^{\text{p}}_{\text{q3}}$	             23.1%      -14.8%     -2.6 %     24.3%

$I2$     AR                   	$\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$       50.1%      -23.6%     -10.0%     52.8%

AR                   	$\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$      137.6%     -31.8%     -12.1%     144.4%

FFNN-W               	$\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$       80.0%      -25.3%     -7.1 %     83.9%

FFNN-W               	$\hat{\mathbf{y}}^{\text{p}}_{\text{mean}}$      164.4%     -31.1%     -17.8%     172.6%

Percentage difference of tactical costs resulting from solving $\textbf{BP}(\hat{\mathbf{y}}^{\text{p}})-\textbf{wBP}(\mathbf{Y})$. The point of reference is the reference used in Analysis 1 ($\textbf{BP}(\mathbf{y}_t)$ for $t=1,\ldots, T$).

We report the total demand values in Table <ref>. For demand instance $I1$, both periodic demand $\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$ overestimate the reference demand. The design built is capable of handling more demand, which results in less outsourced demand. However, they underestimate the periodic demand $\mathbf{y}^{\text{p}}_{\text{q3}}$, which lead to a less important decrease in outsourced demand. In instance $I2$, the periodic $\hat{\mathbf{y}}^{\text{p,AR}}_{\text{max}}$ overestimates the total demand yet it leads to an increase in outsourced demand. This is because $\hat{\mathbf{y}}^{\text{p,AR}}_{\text{max}}$ either overestimates demand for commodities that already lack capacity with $\mathbf{y}^{\text{p}}_{\text{mean}}$, or it underestimates demand for some  commodities for which blocks are then not built in the design and are consequently outsourced at each week.

$I1$     $I2$
Historical data     $\mathbf{y}^{\text{p}}_{\text{mean}}$     -     -

Historical data     $\mathbf{y}^{\text{p}}_{\text{q3}}$     11.9%     11.7%

AR     $\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$     10.2%     0.7%

FFNN-W     $\hat{\mathbf{y}}^{\text{p}}_{\text{max}}$     10.0%     -1.4%

Total periodic demand summed over commodities. The point of reference is indicated by a dash.

§ CONCLUSION AND FUTURE RESEARCH

Tactical planning is essential to freight carriers as it allows to, e.g., design the service network to meet expected demand while minimizing cost. In this work we focused on large-scale tactical planning that is restricted to deterministic models for the sake of computational tractability. Even though estimates of periodic demand is a central input to such models, the associated estimation problem has not been studied in the literature. In this paper we addressed this gap: We formally introduced the periodic demand estimation problem and we proposed a methodology that proceeds in two steps. The first step consists in using a time series forecasting model to predict demand for each period in the tactical planning horizon. The second step defines periodic demand as a solution to a multilevel mathematical program that explicitly connects the estimation problem to the tactical planning problem of interest. This allows to estimate periodic demand such that the costs are minimized. Since the origin-destination demand matrices typically are unbalanced, this can be of importance as the cost of forecast errors is not evenly distributed across commodities.

We reported results for a real large-scale application at the Canadian National Railway Company. The results clearly showed the importance of the periodic demand estimation problem when compared to the approach commonly used in practice. The latter consists in averaging the time series forecasts over the tactical planning horizon. Compared to this practice, the results showed that using another estimate can lead to a substantial reduction in cost.
As expected, the results also showed that the time series forecasting problem is difficult and the forecast errors hence are relatively large.
Nevertheless, the periodic demand estimates that resulted from the proposed methodology still led to costs that were comparable, or even better, than those obtained by using the average demand baseline computed on perfect information (i.e., no forecast error). Moreover, the costs were substantially reduced compared to averaging the forecasts.

In terms of exposition, we chose to limit the methodology to the MCND formulation. However, the methodology applies to other cyclic network design formulations. Similarly, adaptation of the tactical plan in each period (wMCND and wBP formulations) can also be represented differently from this paper. The methodology hinges on the separation between the design variables that are fixed for all periods in the tactical planning horizon while the flow decisions are not. The adaptation of the flow decisions serves as a proxy for operational costs.

Given that we introduced a new problem in this paper, it opens up a number of directions for future research. First, improving the time series forecasts (step one in the methodology).
Second, extending the feasible set of periodic demand values to more general mappings and devise an effective solution approach for this case. The work reported in this paper constituted a first step in addressing the periodic demand estimation problem that hitherto has been overlooked in the literature. We showed that adequately addressing it can lead to important cost reductions.

§ ACKNOWLEDGMENTS

This research was funded by the Canadian National Railway Company (CN) Chair in Optimization of Railway Operations at Université de Montréal and a Collaborative Research and Development Grant from the Natural Sciences and Engineering Research Council of Canada (CRD-477938-14). Moreover, we gratefully acknowledge the close collaboration with personnel from different divisions of CN.

urlstyle[Bergstra et al., 2013]
Bergstra, J., Yamins, D., and Cox, D. D.
Making a science of model search: Hyperparameter optimization in
hundreds of dimensions for vision architectures.
ICML'13: Proceedings of the 30th International Conference on
International Conference on Machine Learning, 28:0 115–123, 2013.

[Chouman et al., 2017]
Chouman, M., Crainic, T. G., and Gendron, B.
Commodity representations and cut-set-based inequalities for
multicommodity capacitated fixed-charge network design.
Transportation Science, 510 (2):0 650–667,

[Crainic, 2000]
Crainic, T. G.
Service network design in freight transportation.
European Journal of Operational Research, 1220
(2):0 272–288, 2000.

[Crainic et al., 2020]
Crainic, T. G., Hewitt, M., Maggioni, F., and Rei, W.
Partial benders decomposition: General methodology and application to
stochastic network design.
Transportation Science, 2020.

[Fiig et al., 2019]
Fiig, T., Weatherford, L. R., and Wittman, M. D.
Can demand forecast accuracy be linked to airline revenue?
Journal of Revenue and Pricing Management, 180
(4):0 291–305, 2019.

[Hochreiter and Schmidhuber, 1997]
Hochreiter, S. and Schmidhuber, J.
Long short-term memory.
Neural computation, 90 (8):0 1735–1780, 1997.

[Holt, 2004]
Holt, C. C.
Forecasting seasonals and trends by exponentially weighted moving
International Journal of Forecasting, 200
(1):0 5–10, 2004.

[Hornik et al., 1989]
Hornik, K., Stinchcombe, M., and White, H.
Multilayer feedforward networks are universal approximators.
Neural networks, 20 (5):0 359–366, 1989.

[Karlaftis and Vlahogianni, 2011]
Karlaftis, M. G. and Vlahogianni, E. I.
Statistical methods versus neural networks in transportation
research: Differences, similarities and some insights.
Transportation Research Part C: Emerging Technologies,
190 (3):0 387–399, 2011.

[Längkvist et al., 2014]
Längkvist, M., Karlsson, L., and Loutfi, A.
A review of unsupervised feature learning and deep learning for
time-series modeling.
Pattern Recognition Letters, 42:0 11–24, 2014.

[Magnanti and Wong, 1984]
Magnanti, T. L. and Wong, R. T.
Network design and transportation planning: Models and algorithms.
Transportation Science, 180 (1):0 1–55, 1984.

[Makridakis et al., 2018]
Makridakis, S., Spiliotis, E., and Assimakopoulos, V.
Statistical and machine learning forecasting methods: Concerns and
ways forward.
PlOS ONE, 130 (3):0 1–26, 2018.

[Mantovani et al., 2018]
Mantovani, S., Morganti, G., Umang, N., Crainic, T. G., Frejinger, E., and
Larsen, E.
The load planning problem for double-stack intermodal trains.
European Journal of Operational Research, 2670
(1):0 107 – 119, 2018.

[Milenković et al., 2019]
Milenković, M., Milosavljevic, N., Bojović, N., and Val, S.
Container flow forecasting through neural networks based on
Operational Research,  1–33, 2019.

[Morganti et al., 2020]
Morganti, G., Crainic, T. G., Frejinger, E., and Ricciardi, N.
Block planning for intermodal rail: Methodology and case study.
Transportation Research Procedia, 47:0 19 – 26, 2020.

[Nguyen et al., 2018]
Nguyen, H., Kieu, L.-M., Wen, T., and Cai, C.
Deep learning methods in transportation domain: A review.
IET Intelligent Transport Systems, 120 (9):0
998–1004, 2018.

[Park et al., 2007]
Park, J. W., Genton, M. G., and Ghosh, S. K.
Censored time series analysis with autoregressive moving average
Canadian Journal of Statistics, 350 (1):0
151–168, 2007.

[Schulze and Prinz, 2009]
Schulze, P. M. and Prinz, A.
Forecasting container transshipment in germany.
Applied Economics, 410 (22):0 2809–2815,

[Sutskever et al., 2014]
Sutskever, I., Vinyals, O., and Le, Q. V.
Sequence to sequence learning with neural networks.
In Advances in Neural Information Processing Systems,  27,
3104–3112, 2014.

[Tsai and Huang, 2017]
Tsai, F.-M. and Huang, L. J.
Using artificial neural networks to predict container flows between
the major ports of Asia.
International Journal of Production Research, 550
(17):0 5001–5010, 2017.

[Weatherford and Belobaba, 2002]
Weatherford, L. and Belobaba, P.
Revenue impacts of fare input and demand forecast accuracy in airline
yield management.
Journal of the Operational Research Society, 530
(8):0 811–821, 2002.

[Winters, 1960]
Winters, P. R.
Forecasting sales by exponentially weighted moving averages.
Management Science, 60 (3):0 324–342, 1960.

§ APPENDIX

§.§ Input features of neural networks

Input Layer 1 It is dedicated to external data features such as weather data and temporal context, for all models. We use real observed weather data described in Section <ref>, to assess their potential in an optimistic setting in regards to their accuracy. At prediction time, we would then rely on weather forecasts.

All models contain two features as temporal context, that is the week number of the forecasted week and the month number of the Monday of the forecasted week. Weather features, for models that include them, consist in the average daily temperature, the accumulated snow (cm) and the accumulated precipitations (mm) of the forecasted week for the main 17 terminals in the network.

Input Layer 2 Features in Input Layer 2 are lagged observed or forecasted (when doing inferences) demand of commodities predicted by the model. The number of lags was found through hyperparameters optimization, and is given in Table <ref> below.

§.§ Hyper-parameters of neural networks

Model     Lags     c]@c@Weather
Features     c]@c@Hidden
layers     c]@c@Size
Layers      Dropout     c]@c@Learning
RNN     3     NO     1     700       0.25     0.1

RNN-W     4     YES     1     260       0.12     0.1

RNN-W-SPLIT1     3     YES     3     460       0.10     0.1

RNN-W-SPLIT2     8     YES     2     220       0.23     0.1

FFNN     3     NO     3     540     0.14     0.01

FFNN-W     4     YES     3     600       0.11     0.1

FFNN-W-SPLIT1     3     YES     2     700       0.06     0.01

FFNN-W-SPLIT2     5     YES     3     580       0.15     0.1

Table of hyperparameters of the neural networks