\documentclass[11pt,reqno]{amsart}
%\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}



\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]





%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{519.8} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{90B06} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Vol. 20, No. 2, pp. 144--145 (2023) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}




\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}
\renewcommand{\tablename}{Table}




\thispagestyle{empty}

\title[Column generation for petrol replenishment problem]{\large A column generation heuristic for a petrol station replenishment problem with complex unloading rules}

\author[B.G. GAVRISH]{{\bf B.G. Gavrish}}
\address{Boris Georgievich Gavrish  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Moscow State University,
\newline\hphantom{iii} Leninskiye Gory, 1-46
\newline\hphantom{iii} 119991, Moscow, Russia}%
\email{\textcolor{blue}{gavrish.boris@list.ru}}%




\author[Y.A. KOCHETOV]{{\bf Y.A. Kochetov}\protect\orcidicon{0000-0001-7123-8647}}
\address{Yury Andreevich Kochetov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Sobolev Institute of Mathematics,
\newline\hphantom{iii} pr. Koptyuga, 4,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}%
\email{\textcolor{blue}{jkochet@math.nsc.ru}}%

\thanks{\sc Gavrish, B.G., Kochetov, Y. A.
A column generation heuristic for a petrol station replenishment problem with complex unloading rules}
\thanks{\copyright \ 2025 Gavrish B.G., Kochetov Y. A.}
\thanks{\rm The research is carried out within the framework of the state contract of the Sobolev Institute of Mathematics (project FWNF-2022-0019).}
\thanks{\it Received January, 1, 2023, Published December, 31, 2023}%



\semrtop \vspace{1cm}
\maketitle {\small

\vspace{-12pt}
\centerline{{\it Communicated by} {\sc P.P. Petrov}
%({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:} The paper addresses a new single-depot multi-trip petrol station replenishment problem originating from a Russian petroleum company. We assume that each underground reservoir at each petrol station can be serviced in only one trip, and the service time of the station depends on the number of hoses on the truck. The truck's compartment closest to its cabin must be emptied last. The goal is to minimize the total size of the heterogeneous truck fleet required to service all the stations. We develop a linear 0-1 formulation for an exact approach and a heuristic based on the column generation method to obtain lower and upper bounds on the optimum. The approach is tested on a semi-synthetic dataset, and the bounds differ by at most two trucks in most cases. \medskip

\noindent{\bf Keywords:} vehicle routing, matheuristic, heterogeneous fleet, MIP formulation.
 \end{quote}
}

\bigskip

\section{Introduction}
In fuel distribution systems, petrol replenishment plays an important role in reducing expenses. These activities include procurement, transportation of products from depots to stations, and inventory management. Classification, models, and exact algorithms for multi-compartment multi-period delivery problems can be found in \cite{Coelho2015}. Typically, the objective of these models is to maximize a profit function or minimize the size of the truck fleet. The corresponding mixed-integer linear programs contain numerous binary variables and are computationally challenging. Therefore, researchers often develop various heuristic approaches. As a rule, the set of admissible truck routes is restricted by implementing an iterative procedure that includes more and more routes for each truck by solving a general formulation for each truck separately \cite{Cornillier2012}, or by splitting the set of stations into clusters based on proximity and solving each subproblem independently \cite{Cornillier2009}.

In this work, we consider a specific case of this problem inspired by the workflow of a Russian petroleum company. In our formulation, we assume that each truck must leave the depot fully loaded, and each of its compartments can be used to fill only one underground reservoir. As a result, we obtain a pure 0-1 formulation where the existence of a feasible solution is non-trivial. To address this challenge, we introduce lower and upper demand bounds for each underground reservoir and allow some reservoirs with a zero lower bound to be ignored. We focus on the single-shift case, but in practice, we deal with a multi-period delivery problem and plan to use the single-shift problem as a subproblem in future work.

We apply the well-known column generation approach to this problem. One of the first papers in which this approach is applied to the petrol replenishment problem is \cite{Avella2004}. They employ the branch-and-price method for this purpose. A recent example using the same technique is \cite{Bani2023}, where the authors apply column generation with a diving subroutine in the branch-and-price method for a real-world multi-depot, multi-period petrol replenishment problem. Other studies have used column generation itself in similar settings as a way to calculate the lower bounds on the optimal objective value. The authors of \cite{Wang2020} study a formulation similar to the one in this paper. An Adaptive Large Neighborhood Search heuristic is evaluated on data from a Chinese petroleum transportation company, compared against exact results from a MILP model and lower bounds from a column generation approach. The authors of \cite{Alyakoob2013} generate daily schedules for vessels (oil tankers) using column generation to construct integer solutions.

We present a novel formulation of the petrol replenishment problem and propose a 0-1 linear model for both exact and heuristic approaches. Each truck leaves the depot fully loaded and can pour fuel from each compartment into only one reservoir. Trucks may be equipped with multiple hoses that allow simultaneous filling of multiple reservoirs. For example, with two hoses, we need to partition the set of compartments into two subsets to minimize the station service time, assuming a constant pouring speed.

The compartment closest to the cabin must be emptied last. We allow split deliveries to reduce the truck fleet size: each station can be visited multiple times per shift, but each reservoir must be filled in a single trip. To solve the problem and minimize the fleet size, we design a column generation heuristic based on the exact solution to the restricted master problem. The pricing problem is the most computationally difficult step of the method. We solve it in parallel for each truck and impose limits on the maximum number of stations per trip and the set of neighboring stations. We study the influence of these limits on the upper and lower bounds on the optimum.    

Our computational results indicate that limiting trips to two stations is a very restrictive constraint. Allowing four or five stations leads to similar results, while three or four stations provide the best balance between running time and solution quality. The gap between the lower and upper bounds is not greater than two trucks in most cases. Therefore, we conclude that the method can generate reasonable daily plans with most trucks completing two trips per shift.

The paper is structured as follows. In Section \ref{genset}, we describe the model's assumptions and explain our motivation for multiple trips for servicing a single station. Section \ref{mathmod} presents a 0-1 linear programming formulation. In Section \ref{colgen}, the master and pricing problems are introduced for the column generation approach.  Section \ref{numsim} presents our computational results for a semi-synthetic dataset based on real-life instances. Section \ref{outro} concludes the work.

\section{Problem statement }\label{genset}

Let $K$ be the set of heterogeneous trucks and $V$ be the set of stations, including the depot. Let $c_{ij}$ be the travel time associated with the arc $(i, j), \ i,j \in V$. Let $\phi$ be the time to pour 1000 liters of petrol into an underground reservoir, and $\theta$ be the time needed for visit documentation per station. Let $\Delta_k$ be the time required to completely fill and prepare an empty truck $k \in K$ for a route. Each truck is characterized by the set of compartments $W_k$, each with capacity $Q_{wk}$ for compartment $w \in W_k$, and by its set of hoses $S_k$. We assume that the hose $0$ is used for the longest pouring time during each station visit. Let $v_{ik}$ be a Boolean parameter indicating whether truck $k$ is compatible with station $i$. Each station has a set of reservoirs, each containing one of the product lines: $(92, 95, 98, 100, \text{DT})$. For each reservoir, we know its product type $p$ and lower and upper bounds $l_p$ and $u_p$ on the quantity of this product to be delivered, where $0 \leq l_p \leq u_p$. Each truck  must complete all of its trips and return to the depot within the time window $[0, H_k]$. The objective is to service all stations using the minimal size truck fleet. 

\subsection{Assumptions}

\begin{enumerate}
    \item There are several types of trucks; the fleet is limited and heterogeneous.
    \item Each truck has between 3 and 6 compartments with known capacities ranging from 5,000 to 16,000 liters. Each truck is loaded at the depot. The quantity of product loaded into each compartment is equal to its capacity.
    \item The product from each compartment is transferred into one underground reservoir only.
    \item Each underground reservoir is filled in a single trip.
    \item Each truck leaving the depot must be completely filled.
    \item The compartment closest to the cabin must be emptied last.
    \item Each truck can simultaneously fill a number of reservoirs, limited by the number of its hoses.
    \item The capacity of each station's underground reservoirs is known (20,000 - 50,000 liters), and each reservoir contains only one product.
    \item The depot has sufficient products for an arbitrary replenishment plan.
    \item Only a single working shift is considered; all station orders must be delivered during this shift.
    \item Several trips can be performed by the same truck as long as they can be done within the allowed working hours per shift.
    \item Travel times are constant. Service time depends on the set of reservoirs being serviced simultaneously by a truck using its hoses.   
    \item Constraints regarding truck-station compatibility and station-station travel are known.
\end{enumerate}

\subsection{Split delivery motivation}
We assume that each station can be visited multiple times per shift, but each reservoir must be filled in only one trip. This assumption broadens the scope of our formulation compared to cases where all station orders must be completed in a single visit (see, for example, \cite{Cornillier2012}).
Savings from this approach were demonstrated in \cite{Dror1990} for the well-known capacitated vehicle routing problem. Now we illustrate how it can lead to a fleet reduction in the context of the new problem.

Suppose that each truck can perform only a single trip. There are three trucks, each with 5 compartments of 5,000 liters. We also consider three stations with the following orders: stations 1 and 3 each have three reservoirs, and station 2 has four reservoirs. We deliver between 5,000 and 10,000 liters to each reservoir.

\begin{figure}[h!]
\centering
\includegraphics[width=0.5\linewidth]{network.png}
\caption{(A): unsplit deliveries, (B): split deliveries} \label{fig:M1}
\end{figure}

We need three trucks to service the stations without split deliveries (one station per truck). Alternatively, the first truck can service all the orders at the first station and the first two orders at the second station, while the second truck fills the remaining reservoirs (see Fig. \ref{fig:M1}). This basic intuition remains valid in other settings and allows for much more flexible plans. Note that this still corresponds to the ‘unsplit-unsplit’ case according to the classification in \cite{Coelho2015}.  

\section{Integer linear programming formulation }\label{mathmod}

 For convenience and to simplify the model, we represent each station $f \in V$ by its set of reservoirs $R_f$ and denote the total set of reservoirs, including the depot, as $I$.  The travel time matrix is recalculated with $c_{ij} = 0$ for reservoirs $i$ and $j$ at the same station. As each truck can perform multiple trips,  we introduce the set $M_k$ of trips for truck $k$. We assume that each set $M_k$ is small, with at most $M$ elements.\\
 
Now we introduce the decision variables. 
\begin{itemize}
    \item $y_k \in \{0,1\}$: whether truck $k$ is used
    \item $g_{km} \in \{0,1\}$: whether truck $k$ conducts trip $m$
    \item $h_{fkm} \in \{0,1\}$: whether truck $k$ visits station $f$ in trip $m$
    \item $z_{ikm} \in \{0, 1\}$: whether truck $k$ fills reservoir $i$ in trip $m$
   \item $x_{ijkm} \in \{0,1\}$: whether truck $k$ travels from $i$ to $j$ in trip $m$
   \item $t_{ikwsm} \in \{0, 1\}$: whether truck $k$ uses  compartment $w$ to fill reservoir $i$ with hose $s$ in trip $m$.
\end{itemize}

According to the model's assumptions, we arrive at the following 0-1 linear formulation.
\begin{flalign} \label{obj1}
    \sum_{k} y_k \rightarrow \min &&
\end{flalign}
\begin{flalign} \label{1ineq1}
    t_{ikwsm} \leq v_{ik} && \forall (i, k, w,s, m) 
\end{flalign}
\begin{flalign}\label{1ineq2}
    \sum_{i,j,m}c_{ij} x_{ijkm} + \sum_{m} \Delta_k g_{km} + \frac{\phi }{1000} \sum_{i,w,m}t_{ikw0m}Q_{wk} + \theta \sum_{f,m}h_{fkm} \leq H_k && \forall k
\end{flalign}
\begin{flalign}\label{1ineq3}
    l_{i} \leq \sum_{w,k,m,s}t_{ikwsm} Q_{wk} \leq u_{i} && \forall i > 0 
\end{flalign}
\begin{flalign}\label{1ineq4}
    g_{km} \leq y_k && \forall (k, m)
\end{flalign}
\begin{flalign}\label{1ineq5}
    \sum_{i,s} t_{ikwsm} = g_{km} && \forall (w, k, m)
\end{flalign}
\begin{flalign}\label{1ineq6}
     \sum_{j} x_{jikm} \geq  \sum_{s}t_{ikwsm} && \forall ( i>0,k,m,w )
\end{flalign}
\begin{flalign}\label{1ineq88}
    t_{ikwsm} \leq z_{ikm} && \forall (i>0,k,w,s, m) \end{flalign}
\begin{flalign}\label{1ineq89}
    \sum_{k,m} z_{ikm} \leq 1 && \forall (i>0) \end{flalign}
\begin{flalign}\label{1ineq7}
    \sum_{j>0} x_{0jkm} = g_{km} && \forall (k, m) \end{flalign}
\begin{flalign}\label{1ineq8}
    \sum_{i} x_{ijkm} = \sum_{i} x_{jikm}  && \forall (j, k, m)
\end{flalign}
\begin{flalign}\label{1ineq9}
    \sum_{i, j \in S} x_{ijkm} \leq |S| - 1 && \forall (k, m, S): 2\leq |S| \leq |W_k|, S \subset I \setminus \{ 0 \} \end{flalign}
\begin{flalign}\label{1ineq10}
    t_{ik0sm} \leq x_{i0km} && \forall (i,k,s,m): i > 0 \end{flalign}
\begin{flalign}\label{1ineq11}
    \sum_s t_{ikwsm} \leq 1 && \forall (i,k,w,m): i > 0 \end{flalign}
\begin{flalign}\label{1ineq12}
    t_{ikwsm} \leq h_{fkm} && \forall (f, i,k,w,m): i \in R_{f} \end{flalign}
\begin{flalign}\label{1ineq13}
    \sum_{i,w:i\in R_f}t_{ikw0m}Q_{wk} \geq \sum_{i,w:i\in R_f}t_{ikwsm}Q_{wk}  && \forall (f,k,s>0,m)
\end{flalign}
In the objective function \eqref{obj1}, we minimize the truck fleet size. To do so, we must not service reservoirs with incompatible trucks \eqref{1ineq1}, must respect the time constraint for each truck \eqref{1ineq2}, and must adhere to the lower and upper bounds for each reservoir \eqref{1ineq3}. \eqref{1ineq4} is a technical condition linking the fleet by variables $y_k$ to the trip variables $g_{km}$. 
Constraint \eqref{1ineq5} ensures that each truck leaves the depot fully loaded. A truck can use its compartments to fill an underground reservoir only if it arrives there \eqref{1ineq6}. 
Constraints \eqref{1ineq88} and \eqref{1ineq89} ensure that each underground reservoir is filled in a single trip.

Additional trip constraints include: \eqref{1ineq7} ensures that a truck leaves the depot if and only if the trip is conducted. Each truck must leave a reservoir if it visits it \eqref{1ineq8}. Condition \eqref{1ineq9} is the well-known subtour elimination constraint, ensuring that the subset $S$ of reservoirs does not exceed the number of the truck’s compartments. Constraint \eqref{1ineq10} guarantees that the compartment closest to the cabin ($w = 0$) is emptied last. \eqref{1ineq11} ensures that one compartment cannot be emptied by two hoses simultaneously. Constraint \eqref{1ineq12} determines whether truck $k$ visits station $f$ in trip $m$. Finally, condition \eqref{1ineq13} ensures that hose $0$ is used for the largest amount of fuel at each station. This is needed for correct calculation of the service time for each trip and each station.

This formulation results in a rapid increase in the number of variables and constraints. To tackle the problem, we design a column generation heuristic based on the exact solution to the restricted master problem.  To eliminate non-promising solutions, we impose limits on the maximum number of stations per trip and the set of neighboring stations. As we show below, this can reduce running time, but may increase the required truck fleet size.

\section{Column generation approach} \label{colgen}

We reformulate the problem \eqref{obj1}--\eqref{1ineq13} as a special case of the set covering problem with a large set $R$ of admissible routes for the trucks. A route $r \in R$ is admissible for truck $k$ if the truck starts from the depot, visits some reservoirs, fills all of them, and returns back to the depot empty. The compartment closest to its cabin is emptied last. For each admissible route $r$, we can calculate its duration $\sigma_{kr}$, which must not exceed the shift length $H_k$, and a 0-1 indicator $b_{ikr}$ indicating whether the fuel was delivered to reservoir $i$. Note that the optimal strategy for using the truck's hoses is applied here to compute parameter $\sigma_{kr}$. We also introduce new variables $g_{kr} \in \{0, 1\}$: whether truck $k$ performs route $r$. It is easy to see that our problem is equivalent to the following one: 

\begin{flalign} \label{master-begin}
    \sum_{k} y_k \rightarrow \min && \end{flalign}
\begin{flalign} \label{}
    \sum_{k, r} b_{ikr}g_{kr} = 1 && \forall i: l_i > 0 \end{flalign}
\begin{flalign} \label{2ineq11}
    \sum_{k, r} b_{ikr} g_{kr} \leq 1 && \forall i: l_i = 0 \end{flalign}
\begin{flalign} \label{master-tricky-constraint}
     g_{kr} \leq y_{k} && \forall (k, r) \end{flalign}
\begin{flalign} \label{master-end}
    \sum_{r} \sigma_{kr} g_{kr}  \leq H_k && \forall k  \end{flalign}
    
Following the column generation framework, we need to find the optimal solution to the linear programming relaxation. To this end, we heuristically generate an initial feasible solution with a small subset of routes and iteratively enlarge it by including the routes with the most negative reduced cost identified through the pricing problem. At the final stage of our procedure, we obtain the lower bound for the optimal value and compute the upper bound as the exact 0-1 solution for the final subset of admissible routes. A similar approach has been successfully applied to many NP-hard combinatorial problems (see, for example,  \cite{Kochetov2023}, \cite{Kochetov2021}, \cite{Avella2012}).

\subsection{Initial solution}

We will use a simple fast heuristic to generate an initial feasible solution. This stage is less critical if a large fleet is available. Otherwise, additional effort is required to verify the problem’s feasibility. Assuming the problem has many solutions, we generate a subset of routes randomly, starting from the depot. The probability of each move depends on the distance to the neighboring station. For a sufficiently large set of admissible routes, we solve the master problem \eqref{master-begin} -- \eqref{master-end} and use the solution as the initial set of routes.


\subsection{Route generation}\label{tripgen}

We replace constraint \eqref{master-tricky-constraint} with a compact one $$     M y_k \geq \sum_r g_{kr},  \qquad \forall k $$
and consider the dual problem:
\begin{flalign} \label{1dual1}
    \sum_{i: l_i > 0} \alpha_i - \sum_{i: l_i = 0} \alpha_i  - \sum_k \delta_k \rightarrow \max &&
\end{flalign}
\begin{flalign}
    M \beta_k  + H_k \gamma_k -\delta_k \leq 1 && \forall k
\end{flalign}
\begin{flalign} \label{target-cg1}
    \sum_{i:l_i>0} b_{ikr} \alpha_i - \sum_{i:l_i=0} b_{ikr} \alpha_i  - \beta_{k} -  \sigma_{kr} \gamma_{k} \leq 0 && \forall (k, r)
\end{flalign}
\begin{flalign} \label{1dual60}
    \alpha_i \in \mathbb{R} && \forall i: l_i > 0
\end{flalign}
\begin{flalign}
     \alpha_i \geq 0, \ \beta_k \geq 0,\  \delta_k \geq 0, \ \gamma_k \geq 0 && \forall k, \ \forall i: l_i = 0
\end{flalign}

For a small subset of routes, we can find the optimal solution $\alpha_i^*, \ \beta_k^*,\  \delta_k^*,  \gamma_k^*$ to this linear program. We now need to verify whether constraint \eqref{target-cg1} holds for all other truck-route pairs. To this end, we create the following pricing problem for each truck with new decision variables: 
\begin{itemize}
    \item $x_{ij}\in \{0, 1\}$: whether the truck visits reservoir $i$ directly after reservoir $j$
    \item $h_f \in \{0, 1\}$: whether station $f$ is visited
    \item $z_{iws}\in \{0, 1\}$: whether compartment $w$ is devoted to filling reservoir $i$ with hose $s$
    \item $y_{i}\in \{0, 1\}$: whether reservoir $i$ is filled
    \item $u_{i} \geq 0$: auxiliary  variables for subtour elimination.
\end{itemize}

The pricing problem for each truck has the following form:

\begin{multline} \label{1building1}
    \sum_{i} \alpha_i^* y_i - \beta_k^* - \\  - \gamma_k^*
    (\sum_{i\neq j}c_{ij} x_{ij} + \Delta_k + \frac{\phi}{1000}  \sum_{w}z_{iw0}Q_w + \theta \sum_{f}h_{f})\rightarrow \max 
\end{multline}
\begin{flalign}
    z_{iws} \leq v_{ki} && \forall (i >0,w,s) 
\end{flalign}
\begin{flalign}
    \sum_{i\neq j}x_{ij} = \sum_{i\neq j}x_{ji} && \forall j \end{flalign}
\begin{flalign}
    \sum_{j>0}x_{0j} = 1 && \end{flalign}
\begin{flalign}
    z_{jws} \leq \sum_{i\neq j}x_{ij} && \forall (j>0, w, s)
\end{flalign}
\begin{flalign}
    \sum_{i\neq j}x_{ij} \leq \sum_{w,s}z_{jws} && \forall j>0
\end{flalign}

\begin{flalign} \label{cycle-ban1}
    u_i - u_j + |W_k| x_{ij} \leq |W_k| -1 && \forall (i>0,j>0)
\end{flalign}

\begin{flalign}
    l_i y_i \leq \sum_{w,s}z_{iws}Q_w \leq u_i y_i && \forall i>0
\end{flalign}

\begin{flalign}
    z_{iws} \leq y_i &&  \forall(i>0, w,s )
\end{flalign}

\begin{flalign}
    \sum_{i,s}z_{iws} = 1 && \forall w
\end{flalign}
\begin{flalign}
    z_{iws} \leq h_{f} && \forall (f, i,w,s): i \in R_{f}
\end{flalign}
\begin{flalign}
    \sum_{i\neq j}c_{ij} x_{ij} + \Delta_k + \frac{\phi}{1000}  \sum_{i,w}z_{iw0}Q_w + \theta \sum_{f}h_{f} \leq H_k && 
\end{flalign}
\begin{flalign}
    \sum_{i,w: i \in R_f}z_{iw0}Q_w \geq \sum_{i,w: i \in R_f}z_{iws}Q_w  && \forall (f,s>0)
\end{flalign}
\begin{flalign} \label{1building10}
    z_{i0s} \leq x_{i0} && \forall (i >0,s)
\end{flalign}

If the optimal value is positive,  we enlarge our subset of admissible routes by including the optimal solution of the pricing problem and repeat the process. Otherwise, we have obtained a lower bound for the minimal fleet size. Typically, this lower bound is extremely tight \cite{Kartak2015}, but the computational time can be significant. The most important question here is the following. How many iterations are needed to obtain the lower bound, and what can be achieved before the stopping criterion is met, especially for large-scale instances with long computation times?



\subsection{Pricing problem resolution}

It is worth noting that, even for instances of moderate size (e.g., $200$ reservoirs or more), the pricing problem becomes computationally hard to be solved exactly, even when using parallel processing with commercial solvers. Iterations of the pricing problem may take up to a few hours. To address this issue, we propose an alternative approach for cases with a small number of stations per route. In real-life applications, the set of available routes between stations is constrained by the company operating the truck fleet. This reduces the service time. According to company rules, a truck can only move to one of the few nearest stations. This implies that the set of feasible routes for each truck is relatively small and can be precomputed for daily use. In other words, this procedure can be considered a preprocessing stage. 

To generate all admissible routes, we use a modified Depth-First Search algorithm to traverse the graph of reservoirs, combined with an optimized Tank Truck Loading Procedure (TTLP) to check the feasibility of each route. The graph traversal is conducted in the following way. We select a truck and start from the depot. At each step, we iterate through all feasible next reservoirs, checking that the extended set of reservoirs can be serviced by the given truck. The proposed TTLP procedure returns a second Boolean value indicating whether no reservoir can be added to the current set so that the truck is able to fill all of them. Consequently, the recursion is stopped in this case. 

The TTLP procedure is presented below. First, corner cases are handled. When the number of reservoirs equals the number of truck compartments, the allocation problem reduces to a matching problem. If the number of reservoirs is smaller, the problem is solved using dynamic programming. $\text{DP}_i(S) \gets true$ if the subset of reservoirs $S$ can be serviced by truck compartments. We also denote by $i', w'$ the last compartment and the last reservoir on the route, respectively. After the procedure, hose allocation is performed by brute force and time feasibility is checked.


%\subsection{Tank Truck Loading Problem (TTLP) algorithm} %\label{ttlp}

\begin{center}
\begin{minipage}{.9\linewidth}   % keeps the code block together
\begin{tabbing}
\qquad\=\qquad\=\qquad\=\qquad\=\qquad\=\qquad\=\kill % ← 6 tab stops
\textbf{procedure} TTLP($Q_w,l_i,u_i$)\`// determines feasibility \\[2pt]

\> $n \gets |\{i\}|$ \\
\> $m \gets |\{w\}|$ \\[2pt]

\> \textbf{if} $n = 0$ \textbf{then} \\
\>\> \textbf{return} $(\text{false},\text{true})$ \\[2pt]

\> \textbf{else if} $n = 1$ \textbf{then} \\
\>\> \textbf{if} $\displaystyle\sum_{w} Q_w < l_1$ \textbf{then} \\
\>\>\> \textbf{return} $(\text{false},\text{false})$ \\
\>\> \textbf{else if} $\displaystyle\sum_{w} Q_w \le u_1$ \textbf{then} \\
\>\>\> \textbf{return} $(\text{true},\text{true})$ \\
\>\> \textbf{else return} $(\text{false},\text{true})$ \\[2pt]

\> \textbf{else if} $n > m$ \textbf{then return} $(\text{false},\text{false})$ \\[2pt]

\> \textbf{else if} $n = m$ \textbf{then} \\
\>\> $M \gets \text{match}(\{i\},\{w\}\mid (i',w')\!\in\! M)$\`// max.\ matching \cite{Kuhn1955}\\
\>\> \textbf{if} $|M| = n$ \textbf{then return} $(\text{true},\text{false})$ \\
\>\> \textbf{else return} $(\text{false},\text{false})$ \\[2pt]

\> \textbf{else}\`// case $1<n<m$ \\[2pt]

\>\> \textbf{if} $Q_{w'} > u_{i'}$ \textbf{then return} $(\text{false},\text{true})$ \\[2pt]

\>\> $F_i   \gets \{\,S\subseteq\{w\}\!\setminus\!w' : l_i \le \sum_{w\in S} Q_w \le u_i\,\}$ \\
\>\> $F_{i'}\gets \{\,S\subseteq\{w\}\!\setminus\!w' : l_i \le Q_{w'} + \sum_{w\in S} Q_w \le u_i\,\}$ \\[2pt]

\>\> $\mathrm{DP}_1(S)\gets\mathbf{1}_{F_1}(S)$ \\[2pt]

\>\> \textbf{for} $i = 2,\dots,n$ \textbf{do} $\mathrm{DP}_i(S)\gets\text{false}$ \\[2pt]

\>\> \textbf{for} $i = 2,\dots,n$ \textbf{do} \\
\>\>\> \textbf{for} $S \subseteq \{w\}\!\setminus\!w'$ \textbf{do} \\
\>\>\>\> \textbf{for} $S' \in F_i$ \textbf{do} \\
\>\>\>\>\> \textbf{if} $S' \subseteq S$ \textbf{then} \\
\>\>\>\>\>\> $\mathrm{DP}_i(S) \gets \mathrm{DP}_i(S)\,\lor\,\mathrm{DP}_{i-1}(S\setminus S')$ \\[2pt]

\>\> $\text{ANS} \gets \mathrm{DP}_n(\{w\}\!\setminus\!w')$ \\[2pt]

\>\> \textbf{return} $(\text{ANS},\text{true})$ \\[4pt]

\textbf{end procedure}
\end{tabbing}
\end{minipage}
\end{center}

The maximization in \eqref{1building1} is then solved by exhaustive search. After having generated all possible routes for each truck, the pricing problem becomes fairly straightforward. It is sufficient to iterate through all the routes at each step of the column generation procedure. As shown in the next section, most admissible routes are rejected by the pricing problem, resulting in a small final subset of routes. This results in a small restricted master problem \eqref{master-begin} -- \eqref{master-end}, which can be easily solved using a solver. This provides an upper bound for the fleet size in the original problem. To obtain the optimal fleet size, we would need to solve the restricted master problem for the entire set of admissible routes.


\section{Numerical simulations}\label{numsim}

The formulation that is studied here has direct practical applications. Petrol companies face vehicle routing problems on a daily basis. Several features of real-life problems studied in this work are worth highlighting.

Firstly, most reservoirs are nearly full, so the total number of orders is relatively low, even for a large network of stations (e.g., out of 120 stations with approximately 5 reservoirs each, only 60 reservoirs may require filling). Secondly, the truck-station compatibility constraints arise due to different truck types. There are two key features: the type of drain (left-sided, right-sided, or double) and the truck type (heavy or not). The type of station is determined by its ability to accommodate these types of truck. Finally, the depot is often located at a considerable distance from most stations. As a result, most trucks complete no more than two trips during a 12-hour shift.

To construct problem instances, we make assumptions based on statistical analysis of real-life data.

We set the number of hoses to 2 for each truck, and assume the same working day duration for all trucks, with $\phi=3$ minutes and $\theta = 15$ minutes.

The sizes of the truck compartment are modeled using a lognormal distribution with parameters $\mu = 9, \; \sigma = 0.3$. A truck has four compartments with probability 0.2, five compartments with probability 0.5, and six compartments with probability 0.3. A truck is considered heavy if its total capacity exceeds 40,000 liters. Forty percent of trucks are two-sided, while 30\% are left-sided and 30\% are right-sided. We assume that it takes three minutes to load and unload 1,000 liters, so $\Delta_k = \frac{\sum_{w} Q_{w}}{1000} \times 3$.

Stations are characterized by their reservoirs. The capacity of each reservoir follows a lognormal distribution with $\mu = 10$ and $\sigma = 0.4$, while the occupancy rate follows a beta distribution with $\alpha = 4.5$ and $\beta = 1.4$. The number of reservoirs for the simulations follows a Poisson distribution with $\mu=6$.
Regarding station types, we suppose that 30\% of the stations can accommodate heavy trucks and, independently, 40\% are two-sided and left and right-sided stations comprise 30\% each. 

The following strategy is used to choose locations. We place the depot at point $(0,0)$ and uniformly distribute the stations within the square of $(100,100) to (200,200)$. Distances are calculated in the Manhattan metric. The travel times are assumed to be proportional to distance and normalized, so it takes one hour for a truck to travel from the depot to $(100,100)$.

\subsection{Simulation strategy}

Our solution approach narrows the set of allowed routes for each truck by introducing two additional restrictions. Firstly, from each station, the truck can only proceed to a certain number of the nearest stations. Secondly, the total number of stations visited by a single truck during one trip from the depot is also usually limited (for example, in \cite{Cornillier2008}, the maximum number is set to $2$). 

We do not wish to set these restrictions arbitrarily. Instead, we perform a grid search to evaluate the performance of different solution strategies on $10$ test instances drawn from the same distribution described above with $N=30$ stations, $K=15$ trucks. 
We also impose a time limit of 30 minutes for each run to ensure practical applicability.

\subsection{Simulation results} 

The proposed algorithm was implemented in C++. The MILP problems were solved using the Google OR-Tools library with the CP-SAT solver backend. The experiments were conducted on an Apple M2 Max chip. In the following tables, $R_1$ denotes the maximum number of stations visited by a truck per trip, and $R_2$ represents the number of available neighboring stations from each station.

\begin{table}[!h]
    \centering
\begin{tabular}{|l|c|c|c|c|}\hline
$R_2 \; \backslash \; R_1$ &
        2 & 3 & 4 & 5 \\ \hline
  12 &   9 / 10 & 6 / 10 & 1 / 10 & 0 / 10 \\ \hline
  8 &   7 / 10 & 9 / 10 & 5 / 10 & 3 / 10 \\ \hline
  4 &  4 / 10 & 9 / 10 & 10 / 10 & 9 / 10 \\ \hline
\end{tabular}
    \caption{Number of successfully completed runs}
    \label{tab1}
\end{table}

The results presented in Table \ref{tab1} demonstrate that excessive constraints can lead to problem infeasibility, as evidenced by only 4 successful runs in the most constrained scenario. Conversely, relaxing these constraints results in increased computational complexity and longer solution times. From a practical implementation perspective, moderately constrained parameters (such as $R_1=3, \; R_2=4$) appear to offer the optimal balance between feasibility and computational efficiency. Restricting routes to 3 stations and 4 neighboring options reduces computational complexity while maintaining solution quality, as overly strict limits ($R_1=2$) lead to infeasibility, and relaxed limits ($R_1=5$) increase runtime without significant fleet reduction.


\begin{table}[h!]
\centering
\begin{tabular}{|l|c|c|c|c|}
  \hline
  $R_2 \; \backslash \; R_1$ &
2 & 3 & 4 & 5 \\ \hline & & & & \\
12 &
  \shortstack{1.16\\(0.78)} &
  \shortstack{13.13\\(6.10)} &
  \shortstack{12.17\\(0.00)} &
  -- \\ \hline & & & & \\
  8  &
  \shortstack{0.64\\(0.42)} &
  \shortstack{6.10\\(3.93)} &
  \shortstack{14.14\\(6.63)} &
  \shortstack{18.17\\(8.64)} \\ \hline & & & & \\

  4  &
  \shortstack{0.20\\(0.17)} &
  \shortstack{1.40\\(1.38)} &
  \shortstack{4.23\\(3.73)} &
  \shortstack{4.24\\(3.27)} \\ \hline
\end{tabular}
\caption{Route preprocessing time average and its standard deviation (minutes)}
\label{tab3}
\end{table}

Table \ref{tab3} provides deeper insight into the computational complexity under varying constraints by examining preprocessing times (the duration required to generate all feasible routes). The results show significant variation in the processing times for different combinations of parameters, with more relaxed constraints generally requiring longer preprocessing periods.

\begin{table}[h!]
\centering
\begin{tabular}{|l|c|c|c|c|}
\hline
$R_2 \; \backslash \; R_1$ &
2 & 3 & 4 & 5 \\ \hline & & & & \\
12 &
  \shortstack{6.22\\(0.74)} &
  \shortstack{5.00\\(0.33)} &
  \shortstack{4.00\\(0.00)} &
  -- \\ \hline & & & & \\
8 &
  \shortstack{6.43\\(0.65)} &
  \shortstack{5.33\\(0.67)} &
  \shortstack{5.00\\(0.40)} &
  \shortstack{4.67\\(0.44)} 
  \\ \hline & & & & \\
4 &
  \shortstack{6.75\\(0.38)} &
  \shortstack{5.56\\(0.72)} &
  \shortstack{5.30\\(0.82)} &
  \shortstack{5.22\\(0.35)} \\ \hline
\end{tabular}
\caption{Average number of trucks used and its standard deviation}
\label{tab2}
\end{table}

A notable finding revealed in Table \ref{tab2} is the substantial reduction in the number of trucks required when the maximum number of stations visited per trip is increased from 2 to 3. This difference suggests that limiting routes to only 2 stations per trip is overly restrictive and may lead to suboptimal resource utilization in the current formulation.

Further analysis focuses on the configurations that demonstrated the best performance (achieving at least 9 out of 10 successful solutions as shown in Table \ref{tab1}). Table \ref{longtab} presents detailed information about the number of columns utilized in each case. While the total number of generated routes varies considerably across instances, the column generation algorithm consistently achieves convergence after adding a relatively small number of additional routes. It is important to mention that the exact solution to the linear problem does not serve as an effective initial solution for the binary integer programming problem, thus validating the choice of a column generation approach for route selection.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.5\linewidth]{barchart.png}
    \caption{Distribution of instances based on the number of trucks with a single route}
    \label{fig:histogram}
\end{figure}

\begin{table}
\centering
\renewcommand{\arraystretch}{1}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Instance &
$R_1=3\,/\,R_2=4$ & $4 / 4$ & $5 / 4$ & $3 / 8$ & $2 / 12$ \\ \hline & & & & & \\
1 &
  \shortstack{115357\\10\\146} &
  \shortstack{349088\\9\\159} &
  \shortstack{524621\\10\\169} &
  \shortstack{490099\\9\\165} &
  \shortstack{66343\\10\\68} \\ \hline & & & & & \\
2 &
  \shortstack{24252\\9\\124} &
  \shortstack{75924\\10\\116} &
  \shortstack{128158\\9\\129} &
  \shortstack{119096\\9\\154} &
  \shortstack{12805\\12\\67} \\ \hline & & & & & \\
3 &
  \shortstack{113074\\11\\154} &
  \shortstack{541929\\11\\92} &
  \shortstack{924221\\12\\74} &
  \shortstack{520067\\11\\202} &
  \shortstack{29945\\12\\59} \\ \hline & & & & & \\
4 &
  \shortstack{20764\\8\\78} &
  \shortstack{61589\\8\\146} &
  \shortstack{66399\\8\\160} &
  \shortstack{84660\\8\\76} &
  \shortstack{5839\\10\\35} \\ \hline & & & & & \\
5 &
  \shortstack{168129\\10\\112} &
  \shortstack{596907\\11\\120} &
  \shortstack{877185\\11\\128} &
  \shortstack{731613\\13\\93} &
  \shortstack{43629\\11\\42} \\ \hline & & & & & \\
6 &
  \shortstack{14437\\12\\50} &
  \shortstack{143420\\9\\118} &
  \shortstack{484035\\10\\173} &
  \shortstack{414452\\11\\91} &
  \shortstack{27445\\11\\63} \\ \hline & & & & & \\
7 &
  \shortstack{330906\\14\\137} &
  \shortstack{1342994\\13\\110} &
  -- & -- &
  \shortstack{98280\\12\\96} \\ \hline & & & & & \\
8 &
  \shortstack{77515\\10\\136} &
  \shortstack{266839\\10\\150} &
  \shortstack{401444\\10\\106} &
  \shortstack{261738\\9\\139} &
  \shortstack{21360\\11\\61} \\ \hline & & & & & \\
9 &
  -- &
  \shortstack{204061\\9\\152} &
  \shortstack{296980\\9\\177} &
  \shortstack{162125\\9\\62} &
  -- \\ \hline & & & & & \\
10 &
  \shortstack{80379\\11\\103} &
  \shortstack{298804\\9\\136} &
  \shortstack{875230\\9\\111} &
  \shortstack{298362\\9\\73} &
  \shortstack{19563\\12\\39} \\ \hline
\end{tabular}
\caption{Number of columns across instances. \textsuperscript{1}}
\small\textsuperscript{1} In each cell, top–to–bottom: total generated routes, initial solution routes, added routes via CG.
\label{longtab}
\end{table}


A crucial practical consideration is the distribution of routes among trucks during their 12-hour shifts. The solution results indicate that the majority of trucks complete two routes, efficiently utilizing their designated time. Only a small number of vehicles, as visualized in Fig. \ref{fig:histogram}, are assigned to single routes per shift, demonstrating that the solution algorithm successfully handles complex unloading constraints while producing efficient schedules.


\section{Conclusion}\label{outro}

We studied a new formulation of the petrol replenishment problem for a single shift. Based on the specific characteristics of the petrol replenishment in Russia, we introduced a set of novel assumptions that significantly affected the nature of the problem. We presented a new 0-1 linear formulation for the problem, which includes complex unloading rules, and developed a column generation heuristic. Numerical simulations on a semisynthetic dataset demon\-strate that the proposed procedure can produce near-optimal solutions.

In future work, it would be of interest to explore generalizations of the problem discussed here. Extending the planning horizon to a week would make the setting more realistic, as current orders are interconnected with past and future ones. In addition, travel times may differ between night and day shifts. Consequently, some stations may need to be visited only during night shifts. Determining the optimal partition of stations between night and day shifts is also a promising research direction. \\
% \\
%\textbf{Acknowledgement:} 
%The research is carried out within the framework of the %state contract of the Sobolev Institute of Mathematics %(project FWNF-2022-0019). 

\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).

\bibitem{Coelho2015} L.C.~Coelho, G.~Laporte,
\href{https://doi.org/10.1016/j.ejor.2014.10.059}{\it Classification, models and exact algorithms for multi-compartment delivery problems},
European Journal of Operational Research, {\bf 242}:3 (2015), 854--864.


\bibitem{Cornillier2012} F.~Cornillier, F.~Boctor, J.~Renaud,
\href{http://dx.doi.org/10.1016/j.ejor.2012.02.007}{\it Heuristics for the multi-depot petrol station replenishment problem with time windows},
European Journal of Operational Research, {\bf 220}:2 (2012), 361--369.

\bibitem{Cornillier2009} F.~Cornillier, G.~Laporte, F.F.~Boctor, J.~Renaud,
\href{http://dx.doi.org/10.1016/j.cor.2007.11.007}{\it The petrol station replenishment problem with time windows},
Computers \&\ Operations Research, {\bf 36}:3 (2009), 919--935.

\bibitem{Avella2004} P.~Avella, M.~Boccia, A.~Sforza,
\href{http://dx.doi.org/10.1016/S0377-2217(02)00676-8}{\it Solving a fuel delivery problem by heuristic and exact approaches},
European Journal of Operational Research, {\bf 152}:1 (2004), 170--179.


\bibitem{Bani2023} A.~Bani, I.~El\,Hallaoui, A.I.~Corréa, A.~Tahir,
\href{http://dx.doi.org/10.1016/j.ejor.2023.04.016}{\it Solving a real-world multi-depot multi-period petrol replenishment problem with complex loading constraints},
European Journal of Operational Research, {\bf 311}:1 (2023), 154--172.


\bibitem{Wang2020} L.~Wang, J.~Kinable, T.~van~Woensel,
\href{https://doi.org/10.1016/j.cor.2020.104904}{\it The fuel replenishment problem: a split-delivery multi-compartment vehicle routing problem with multiple trips},
Computers \&\ Operations Research, {\bf 118} (2020), 104904.


\bibitem{Alyakoob2013} S.M.~Al-Yakoob, H.D.~Sherali,
\href{https://doi.org/10.1016/j.apm.2012.05.028}{\it A column generation approach for determining optimal fleet mix, schedules, and transshipment facility locations for a vessel transportation problem},
Applied Mathematical Modelling, {\bf 37}:4 (2013), 2374--2387.

\bibitem{Dror1990} M.~Dror, P.~Trudeau,
\href{http://dx.doi.org/10.1002/nav.3800370304}{\it Split delivery routing},
Naval Research Logistics, {\bf 37}:3 (1990), 383--402.


\bibitem{Kochetov2023} Yu.A.~Kochetov, A.V.~Ratushnyi,
\href{http://dx.doi.org/10.21538/0134-4889-2024-30-1-109-127}{\it Upper and lower bounds for the optimum in a temporal bin packing problem},
Trudy Instituta Matematiki i Mekhaniki UrO RAN, {\bf 30}:1 (2023), 109--127.

\bibitem{Kochetov2021} Y.~Kochetov, A.~Kondakov,
\href{http://dx.doi.org/10.2298/YJOR200117009K}{\it A hybrid VNS matheuristic for a bin packing problem with a color constraint},
YUJOR, {\bf 31}:3 (2021), 285--298.

\bibitem{Avella2012} P.~Avella, M.~Boccia, S.~Salerno, I.~Vasilyev,
\href{http://dx.doi.org/10.1016/j.cor.2011.09.016}{\it An aggregation heuristic for large scale $p$-median problem},
Computers \&\ Operations Research, {\bf 39}:7 (2012), 1625--1632.


\bibitem{Kartak2015} V.M.~Kartak, A.V.~Ripatti, G.~Scheithauer, S.~Kurz,
\href{http://dx.doi.org/10.1016/j.dam.2015.02.020}{\it Minimal proper non-IRUP instances of the one-dimensional cutting stock problem},
Discrete Applied Mathematics, {\bf 187} (2015), 120--129.

\bibitem{Kuhn1955} H.W.~Kuhn,
\href{https://doi.org/10.1002/nav.3800020109}{\it The Hungarian method for the assignment problem},
Naval Research Logistics Quarterly, {\bf 2}:1–2 (1955), 83--97.


\bibitem{Cornillier2008} F.~Cornillier, F.F.~Boctor, G.~Laporte, J.~Renaud,
\href{http://dx.doi.org/10.1057/palgrave.jors.2602374}{\it An exact algorithm for the petrol station replenishment problem},
Journal of the Operational Research Society, {\bf 59}:5 (2008), 607--615.



\end{thebibliography}


\end{document}
