Minimizing the waiting time for a one-way shuttle service

TL;DR

This paper develops efficient algorithms to optimize shuttle schedules for minimizing user waiting times in a terminal with known, variable arrival rates, combining theoretical guarantees with practical performance evaluations.

Contribution

It introduces novel algorithms for shuttle scheduling that incorporate various optimization techniques and provide proven performance guarantees.

Findings

Algorithms achieve near-optimal waiting time reductions.

Methods perform well in practical numerical experiments.

Techniques include convex optimization and shortest path approaches.

Abstract

Consider a terminal in which users arrive continuously over a finite period of time at a variable rate known in advance. A fleet of shuttles has to carry the users over a fixed trip. What is the shuttle schedule that minimizes their waiting time? This is the question addressed in the present paper. We propose efficient algorithms for several variations of this question with proven performance guarantees. The techniques used are of various types (convex optimization, shortest paths,...). The paper ends with numerical experiments showing that most of our algorithms behave also well in practice.

Tables

Table 1. Table 1. Numerical results for problem P no return max subscript superscript absent no return {}^{\max}_{\mbox{{\tiny no return}}}

$D$	$S$	$\epsilon$	LB	UB	gap	CPU
					(%)	( s )
1P	100	${10}^{-4}$	27.2	27.2	0.0	0
1P	150	${10}^{-4}$	18.1	18.1	0.0	0
1P	200	${10}^{-4}$	13.6	13.6	0.0	0
1P	250	${10}^{-4}$	11.0	11.0	0.0	0
2P	100	${10}^{-4}$	27.0	27.0	0.0	0
2P	150	${10}^{-4}$	18.0	18.0	0.0	0
2P	200	${10}^{-4}$	13.5	13.5	0.0	0
2P	250	${10}^{-4}$	10.8	10.8	0.0	0

Table 2. Table 2. Numerical results for problem P no return ave subscript superscript absent ave no return {}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}

$D$	$S$	$M$	LB	UB	gap	CPU
					(%)	( s )
1P	100	32	17.3	19.2	10.0	34
1P	100	128	18.7	19.2	2.5	1930
1P	200	32	7.7	9.6	19.4	70
1P	200	128	9.1	9.6	5.0	4035
2P	100	32	17.5	19.4	9.9	38
2P	100	128	18.9	19.4	2.5	2387
2P	200	32	7.9	9.7	19.2	76
2P	200	128	9.2	9.7	5.0	4463

Table 3. Table 3. Numerical results for problem P return max subscript superscript absent return {}^{\max}_{\mbox{{\tiny return}}}

$D$	$S$	$M$	LB	UB	gap	CPU
					(%)	( s )
1P∗	1	16	168.6	214.2	21.3	104
1P∗	1	32	184.5	210.8	12.5	1654
2P∗	1	16	101.0	131.0	22.9	106
2P∗	1	32	106.9	126.3	15.4	1848

Equations

$$\tau\colon y\in[0,D(T)]\longmapsto\left\{\begin{array}[]{ll}\inf\left\{t\in[0,T]\colon D(t)>y\right\}\in[0,T]&\mbox{if}\;y<D(T)\\ T&\mbox{otherwise}\end{array}\right.$$

$$\mathrm{bar}\tau :y\in [0,D(T)]⟼\inf \left\{t\in [0,T]:D(t)\ge y\right\}\in [0,T],$$

$$g^{\max }(\mathbf{d},\mathbf{y})$$

$$g^{\mathrm{ave}}(\mathbf{d},\mathbf{y})$$

$${}_{\text{mbox}noreturn}$$

$$g^{\max}(\boldsymbol{d},\boldsymbol{y})\geq g^{\max}(\boldsymbol{d}^{\prime},\boldsymbol{y})=\max_{j}\big{(}d^{\prime}_{j}-\tau(y_{j-1})\big{)}.$$

$$d_{j}=\max\big{(}d_{j-1},\mathrm{bar}\tau(y_{j})+\nu(y_{j}-y_{j-1})\big{)}$$

$$\mathrm{bar}y_{j}=\left\{\begin{array}[]{ll}y_{j}-\varepsilon&\mbox{for $j=j_{0},\ldots,j_{1}-1$}\\ y_{j}&\mbox{otherwise}\end{array}\right.$$

$$\mathrm{bar}d_{j}=\left\{\begin{array}[]{ll}\max\big{(}\mathrm{bar}d_{j-1},\mathrm{bar}\tau(\mathrm{bar}y_{j})+\nu(\mathrm{bar}y_{j}-\mathrm{bar}y_{j-1})\big{)}&\mbox{for $j=j_{0},\ldots,j_{1}$}\\ d_{j}&\mbox{otherwise,}\end{array}\right.$$

$$\sum_{j=j_{0}}^{j_{1}}\int_{\mathrm{bar}y_{j-1}}^{\mathrm{bar}y_{j}}\big{(}\mathrm{bar}d_{j}-\mathrm{bar}\tau(u)\big{)}du\leq\int_{\mathrm{bar}y_{j_{0}-1}}^{\mathrm{bar}y_{j_{1}}}\big{(}\mathrm{bar}d_{j_{1}}-\mathrm{bar}\tau(u)\big{)}du<\int_{y_{j_{0}-1}}^{y_{j_{1}}}\big{(}d_{j_{1}}-\mathrm{bar}\tau(u)\big{)}du=\sum_{j=j_{0}}^{j_{1}}\int_{y_{j-1}}^{y_{j}}\big{(}d_{j}-\mathrm{bar}\tau(u)\big{)}du,$$

$${}_{\text{mbox}return}$$

$$d_s,d_{s+S},d_{s+2S},\dots$$

$$\frac{\nu D(T)}{S}+\left(\lceil \frac{D(T)}{CS}\rceil -1\right)\pi .$$

$$g^{\mathrm{ave}}(\mathbf{d},\mathbf{y})\le OPT+O\left(\frac{S^2}{M}\right)$$

$$\left(2\lceil \frac{T}{\pi }\rceil +1\right)S+\left(\frac{\nu }{\pi }+\frac{1}{C}\right)D(T).$$

$$\left(\frac{k}{S}-1\right)\pi \le T+\frac{\nu D(T)}{S}+\frac{\pi D(T)}{CS},$$

$$k\le \frac{TS}{\pi }+\frac{\nu D(T)}{\pi }+\frac{D(T)}{C}+S.$$

$$g_{\text{mbox}return}^{\max }(\mathbf{d},\mathbf{y})\le OPT+O\left(\frac{S^2}{M}\right)$$

$$\begin{array}[]{rll}\operatorname{Min}&\displaystyle{\max_{s\in\left\{1,\ldots,S\right\}}\left(\nu D_{s}+\left(\left\lceil\frac{D_{s}}{C}\right\rceil-1\right)\pi\right)}\\ \mbox{s.c.}&\displaystyle{\sum_{s=1}^{S}D_{s}=D}&\\ &D_{s}\geq 0&s=1,\ldots,S.\end{array}$$

$$\frac{1}{D}\sum _{j=1}^{+\infty }\left((j-1)\pi +\nu \sum _{i=1}^j{x}_i\right)x_j=\frac{1}{D}\left(\sum _{j=1}^{+\infty }(j-1)\pi x_j+\frac{1}{2}\nu \sum _{j=1}^{+\infty }{x}_j^2\right)+\frac{\nu D}{2},$$

$$P(D)$$

$$\begin{array}[]{l}x_{0}^{*}=0\\ x_{j}^{*}=\left\{\begin{array}[]{ll}C&\quad\mbox{if $j\leq a$,}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \displaystyle{\frac{D-aC}{\theta(a)-a}+\frac{\pi}{\nu}\left(\frac{a+\theta(a)+1}{2}-j\right)}&\quad\mbox{if $a+1\leq j\leq\theta(a)$,}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 0&\quad\mbox{otherwise,}\end{array}\right.\end{array}$$

$$\theta (a)=a+\lceil \frac{-1+\sqrt{1+\frac{8\nu C}{\pi }(D-aC)}}{2}\rceil .$$

$$P_0^n(D)$$

$$\nu x_j^{\ast }+(j-1)\pi +\lambda +{\mu }_j-{\omega }_j=0\text{mbox}and{\omega }_j{x}_j^{\ast }={\mu }_j(x_j^{\ast }-C)=0.$$

$$x_{j}^{*}=\left\{\begin{array}[]{ll}C&\quad\mbox{if $j\leq a$,}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ \displaystyle{\frac{D-aC}{b-a}+\frac{\pi}{\nu}\left(\frac{a+b+1}{2}-j\right)}&\quad\mbox{if $a+1\leq j\leq b$,}\vskip 3.0pt plus 1.0pt minus 1.0pt\\ 0&\quad\mbox{otherwise.}\end{array}\right.$$

$$(b-a)(b-a-1)<\frac{2\nu }{\pi }(D-aC).$$

$$(b-a)(b-a+1)\ge \frac{2\nu }{\pi }(D-aC).$$

$$-\frac{1}{2}+\frac{\sqrt{1+\frac{8\nu }{\pi }(D-aC)}}{2}\le b-a<\frac{1}{2}+\frac{\sqrt{1+\frac{8\nu }{\pi }(D-aC)}}{2}$$

$$v_0(D-\epsilon )\le v(D)$$

$$x^{\prime}_{j}=\left\{\begin{array}[]{ll}\min(x_{j},D-\varepsilon-\sum_{i=1}^{j-1}x_{i}^{\prime})&\mbox{for $j\leq N_{\varepsilon}$}\\ 0&\mbox{for $j\geq N_{\varepsilon}+1$.}\end{array}\right.$$

Full text

1. Introduction

The original motivation of this paper comes from a partnership of the authors with Eurotunnel, the company operating the tunnel under the Channel. Eurotunnel is currently facing an increasing congestion due to the trucks waiting in the terminal before being loaded in the shuttles. A way to address this issue consists in scheduling the shuttles so that the trucks do not wait too long in the terminal.

In railway transportation, a traditional point of view considers that the demand can be smoothed by offering sufficiently enough departures over a day. Timetabling is then guided by other considerations, such as robustness, maintainability, or rolling stock. For instance, Swiss, Dutch and German companies usually design periodic timetables, which present many advantages [6, 9]. The way to optimize this kind of timetables has been the topic of many researches, initiated by Serafini and Ukovich [15] and by Voorhoeve [16] explicitly in the railway context, see [10, 12, 13, 14] for further works. In the context of periodic timetables, a way to adapt the schedules to a demand with strong variations consists in inserting new departures at peak-hours and deleting departures when the demand is low.

Since the trip of the trucks in the tunnel is a small part of their whole journey, it is a reasonable approximation to assume that they cannot choose their arrival time in the terminal. Moreover, increasing the size of the fleet is not always doable in practice (the shuttles are expensive and the tunnel is used by other vehicles, which limits the maximal number of shuttle trips over a day). We face thus a different problem than the one addressed in the aforementioned literature: the demand is assumed to be fixed and nonelastic to the departures times, and the number of shuttles cannot be adjusted to the demand. Given a fleet of shuttles and a demand of transportation known in advance, the problem consists in designing a schedule for the shuttles that minimizes the waiting time of the users. There are timetabling problems with similar features, see [1, 2, 3, 5, 4, 8] for instance, but these articles are at a more macroscopic level than what we require to solve our problem. Moreover, in the present work, the schedules have to be designed in an offline manner. In a transportation context, and especially for Eurotunnel, computing the schedule in advance is mandatory.

We study several versions of the problem, mainly according to two features. The first feature is whether the shuttles are allowed to come back at the terminal after having realized the trip. The second feature is the objective function. We consider in turn the following two quantities to be minimized: the maximum waiting time and the average waiting time. The first objective is perhaps a fairer one regarding the users, while the second one is relevant for the global efficiency.

It seems that the question we address in the present paper is new. Moreover, it may be relevant for any situation where a demand, known in advance, has to be processed by batches and for which we want to minimize the processing time. This kind of situation is often met in chemical industry. An example whose motivation is very close to ours considers a test to be performed on samples, which arrive continuously [11]. The test takes a certain amount of time and can be performed on several samples simultaneously. The question is then to determine the test schedule that minimizes the processing time.

We propose efficient algorithms for the different versions. “Efficient” here means “theoretically efficient”, according to their time complexity and the performance guarantee. It also means “practically efficient” for almost all cases, as shown by numerical experiments conducted on real-world instances. It might also be worth noting that, depending on the version considered, the proof techniques rely on various fields of optimization (convexity, Karush-Kuhn-Tucker conditions, binary search, domination results, shortest paths in finite graphs,…).

2. Model

2.1. The problems

We are given a fleet of $S$ shuttles, for which departure times have to be determined. All shuttles have a capacity $C\geq 0$ and are situated in the same loading terminal at the beginning of the day. The users are infinitesimal and arrive continuously in the terminal over a finite period of time, modeled as the interval $[0,T]$, following the cumulative nondecreasing function $D:[0,T]\rightarrow\mathbb{R}_{+}$, where $D(t)$ is the total number of users arrived in the terminal during the interval $[0,t]$. We assume throughout the paper that $D(T)>0$. The shuttles have to carry the users over a fixed trip.

When the users arrive in the terminal, they enter a queue. This queue closes when all the users who will leave with the next shuttle have arrived in the queue and users can enter a new queue only if the previous one is closed (this is how it works at Eurotunnel). When a queue is closed, the users in that queue can start boarding the shuttle. The process in illustrated on Figure 1. Loading a shuttle with a total of $x$ users takes a time $\nu x$. Note that setting $\nu$ to zero allows to model the case where the users do not have to wait for the last user before boarding. Even if no users arrive strictly after time $T$, loading and departures are allowed after that instant.

Two possibilities are considered regarding the shuttles. Either return is not allowed: once the shuttles leave, they never come back to the terminal; or it is allowed: once the shuttles leave, they come back to the terminal after a time equal to $\pi\geq 0$. Two objective functions to be minimized are considered: the maximum waiting time and the average waiting time.

We have thus four problems:

• •

P${}^{\max}_{\mbox{{\tiny no return}}}$, which consists of not allowing return and minimizing the maximum waiting time.

• •

P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, which consists of not allowing return and minimizing the average waiting time.

• •

P${}^{\max}_{\mbox{{\tiny return}}}$, which consists of allowing return and minimizing the maximum waiting time.

• •

P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$, which consists of allowing return and minimizing the average waiting time.

Practical constraints impose that overtake is not possible and thus, when return is allowed, the departure orders of the shuttles remain the same over the whole period. It is nevertheless possible to have simultaneous trips. This is an approximation in the case of Eurotunnel (we neglect the security distance and the length of the shuttles). For other situations, as the chemical application mentioned in the introduction, it may match what is met in practice.

2.2. The demand

Throughout the paper, we assume that $D(\makebox[4.30554pt]{\cdot})$ is upper semicontinuous. It allows to model discontinuity in the arrival process (batch of users arriving simultaneously). Yet, a weaker requirement could lead to mathematical difficulties, e.g., nonexistence of optimal solutions even for very simple cases.

The pseudo-inverses of $D(\makebox[4.30554pt]{\cdot})$, defined by

[TABLE]

and

[TABLE]

play an important role in the paper. Note that they are nondecreasing functions and that $\tau(y)\geq\mathrm{bar}\tau(y)$ for all $y\in[0,D(T)]$. Times $\tau(y)$ and $\mathrm{bar}\tau(y)$ are respectively interpreted as the arrival times of the first user after the $y$ firsts and of the last user of these $y$ firsts. Since $D(\makebox[4.30554pt]{\cdot})$ is upper semicontinuous, we have the following properties, with proofs for sake of completeness.

Lemma 1**.**

We have $D(\tau(y))\geq D(\mathrm{bar}\tau(y))\geq y$ for every $y\in[0,D(T)]$.

Proof.

Since $D(\makebox[4.30554pt]{\cdot})$ is nondecreasing, the first inequality is a direct consequence of the inequality $\tau(y)\geq\mathrm{bar}\tau(y)$, which is obvious from the definition. To prove the second inequality, consider $(t_{n})$ a nonincreasing sequence converging toward $\mathrm{bar}\tau(y)$ such that $D(t_{n})\geq y$ for all $n$. By the upper semicontinuity of $D(\makebox[4.30554pt]{\cdot})$, we get then $D(\mathrm{bar}\tau(y))\geq y$. ∎

Lemma 2**.**

$\tau(\makebox[4.30554pt]{\cdot})$* is upper semicontinuous and $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is lower semicontinuous.*

Proof.

We first prove that $\tau(\makebox[4.30554pt]{\cdot})$ is upper semicontinuous. Let $\alpha$ be some real number such that $\{y\colon\tau(y)<\alpha\}$ is nonempty and take from it an arbitrary element $y_{0}$. We want to prove that $\{y\colon\tau(y)<\alpha\}$ is open for the induced topology on $[0,D(T)]$. If $y_{0}=D(T)$, then this set is $[0,D(T)]$ and thus open. Otherwise, by the definition of $\tau(\makebox[4.30554pt]{\cdot})$, we know that there exists $t_{0}<\alpha$ such that $D(t_{0})>y_{0}$. For any element $y$ in $[0,D(t_{0}))$, we have $\tau(y)\leq t_{0}$, and thus $[0,D(t_{0}))$ is an open set containing $y_{0}$ and fully contained in $\{y\colon\tau(y)<\alpha\}$. The set $\{y\colon\tau(y)<\alpha\}$ is thus an open set of $[0,D(T)]$ for every real number $\alpha$, which precisely means that $\tau(\makebox[4.30554pt]{\cdot})$ is upper semicontinuous.

We prove now that $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is lower semicontinuous. Let $\alpha$ be some real number. Consider a converging sequence $(y_{n})$ such that $\mathrm{bar}\tau(y_{n})\leq\alpha$ for all $n$. For every fixed $n$, there exists thus a sequence $(t_{n,k})$ indexed by $k$ such that $D(t_{n,k})\geq y_{n}$ and $t_{n,k}\leq\alpha+\frac{1}{k}$ for all $k$. Now, consider the sequence $(t_{n,n})$. It is such that $D(t_{n,n})\geq y_{n}$ and $t_{n,n}\leq\alpha+\frac{1}{n}$ for all $n$. Since $[0,T]$ is compact, we can extract a nonincreasing converging subsequence $(t_{n})$ from the sequence $(t_{n,n})$ such that $D(t_{n})$ converges towards some real number nonsmaller than $\lim_{n\to\infty}y_{n}$ with $t_{n}\leq\alpha+\frac{1}{n}$ for all $n$. It implies that $\mathrm{bar}\tau(\lim_{n\to\infty}y_{n})\leq\alpha$, which means that $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is lower semicontinuous. ∎

Lemma 3**.**

If $D(\makebox[4.30554pt]{\cdot})$ is increasing, then $\tau(y)=\mathrm{bar}\tau(y)$ for every $y\in[0,D(T)]$.

Proof.

If $\mathrm{bar}\tau(y)=T$, then the equality is obvious. We can thus assume that $\mathrm{bar}\tau(y)<T$. For every $t>\mathrm{bar}\tau(y)$, we have $D(t)>D(\mathrm{bar}\tau(y))$ since $D(\makebox[4.30554pt]{\cdot})$ is increasing, and Lemma 1 implies that $D(t)>y$. By definition of $\tau(\makebox[4.30554pt]{\cdot})$, we have $\tau(y)\leq\mathrm{bar}\tau(y)$. The reverse inequality being clear from the definitions, we get the result. ∎

2.3. Mathematical model

For the four problems P${}^{\max}_{\mbox{{\tiny no return}}}$, P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, P${}^{\max}_{\mbox{{\tiny return}}}$, and P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$, a feasible solution is characterized by two nondecreasing sequences of nonnegative real numbers $\boldsymbol{d}=d_{1},d_{2},\ldots$ and $\boldsymbol{y}=y_{1},y_{2},\ldots$. The $d_{j}$’s are the successive departure times of the shuttles, and the $y_{j}$’s are their successive cumulative loads: the $j$th departure occurs at time $d_{j}$ with a load of $y_{j}-y_{j-1}$ users, where we set $y_{0}=0$.

Denote by $g^{\max}(\boldsymbol{d},\boldsymbol{y})$ the value of the maximum waiting time and by $g^{\operatorname{ave}}(\boldsymbol{d},\boldsymbol{y})$ the value of the average waiting time. There are explicit expressions of these objective functions. Note that $\tau(y_{j})$ can be interpreted as the first arrival time of a user leaving with the “$(j+1)$th shuttle”.

[TABLE]

where the indices $j$ range over all departures.

Problems P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ can be written under the following form,

[TABLE]

where $g(\makebox[4.30554pt]{\cdot})$ is either $g^{\max}(\makebox[4.30554pt]{\cdot})$ or $g^{\operatorname{ave}}(\makebox[4.30554pt]{\cdot})$. Constraint (i) ensures that the total amount of users in any shuttle does not exceed the shuttle capacity. Constraint (ii) ensures that the indices of the $y_{j}$ variables are consistent. Constraint (iii) ensures that the shuttles do not overtake. Constraint (iv) ensures that every user eventually leaves the terminal in a shuttle. Constraint (v) ensures that the departure time of a shuttle occurs once the last user of this shuttle has arrived and the loading is over.

Problems P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ always admit optimal solutions when they are feasible, i.e., when $CS\geq D(T)$. Indeed, $\mathrm{bar}\tau(y_{j})+\nu(y_{j}-y_{j-1})$ is upper-bounded by $T+\nu C$ and adding a constraint $d_{j}\leq T+\nu C$ for all $j$ does not change the optimal value; since $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is lower semicontinuous (Lemma 2), the set of feasible solutions of the optimization problem obtained with this new contraint is compact; its objective function is lower semicontinuous (and even continuous in the case of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$).

The following properties for P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ will be useful in some proofs.

Claim 1**.**

Replacing $g^{\max}(\makebox[4.30554pt]{\cdot})$ by \max_{j}\big{(}d_{j}-\tau(y_{j-1})\big{)} does not change the optimal value of P${}^{\max}_{\mbox{{\tiny no return}}}$.

Proof.

Let $(\boldsymbol{d},\boldsymbol{y})$ be a feasible solution of P${}^{\max}_{\mbox{{\tiny no return}}}$. We are going to build a feasible solution $(\boldsymbol{d}^{\prime},\boldsymbol{y})$ (with the same $\boldsymbol{y}$) such that

[TABLE]

We set $d_{1}^{\prime}=\mathrm{bar}\tau(y_{1})+\nu y_{1}$ and define inductively $d^{\prime}_{j}=\max(d_{j-1}^{\prime},\mathrm{bar}\tau(y_{j})+\nu(y_{j}-y_{j-1}))$. We have $d_{j}^{\prime}\leq d_{j}$ for all $j$ and it implies the inequality in (1). Let us prove the equality of (1): if \max_{j}\big{(}d^{\prime}_{j}-\tau(y_{j-1})\big{)} is attained for a $\mathrm{bar}\jmath$ such that $y_{\mathrm{bar}\jmath-1}<D(T)$, then there exists $k\geq\mathrm{bar}\jmath$ such that $y_{k}>y_{k-1}=y_{\mathrm{bar}\jmath-1}$ and $d_{k}^{\prime}\geq d_{\mathrm{bar}\jmath}^{\prime}$, which means that the maximum is also attained for a $k$ such that $y_{k}>y_{k-1}$; and if \max_{j}\big{(}d^{\prime}_{j}-\tau(y_{j-1})\big{)} is attained for a $\mathrm{bar}\jmath$ such that $y_{\mathrm{bar}\jmath-1}=D(T)$, then there exists $\ell\leq\mathrm{bar}\jmath-1$ such that $y_{\ell-1}<y_{\ell}=y_{\mathrm{bar}\jmath-1}$ and by construction $d^{\prime}_{\ell}=d^{\prime}_{\ell+1}=\cdots=d^{\prime}_{S}$ (since $y_{\ell}=y_{\ell+1}=\cdots=y_{S}=D(T)$), which means that the maximum is also attained for an $\ell$ such that $y_{\ell}>y_{\ell-1}$. ∎

Claim 2**.**

If $D(\makebox[4.30554pt]{\cdot})$ is increasing, for any of P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, there is an optimal solution such that $d_{j}=\mathrm{bar}\tau(y_{j})+\nu(y_{j}-y_{j-1})$ for all $j\in\{1,\ldots,S\}.$

Proof.

Let $(\boldsymbol{d},\boldsymbol{y})$ be an optimal solution (we do not care which objective function is used yet). Without loss of generality, we can assume that $d_{1}=\mathrm{bar}\tau(y_{1})+\nu y_{1}$ and that for all $j\in\{2,\ldots,S\}$ we have

[TABLE]

(just redefine $d_{j}$ according to these equalities if necessary). When $\nu=0$, a straightforward induction on $j$ shows that we have then always $d_{j}=\mathrm{bar}\tau(y_{j})$. We can thus assume that $\nu>0$.

Suppose for a contradiction that there is a $j$ such that $d_{j}>\mathrm{bar}\tau(y_{j})+\nu(y_{j}-y_{j-1})$. Denote by $j_{1}$ the smallest index for which this inequality holds. We necessarily have $d_{j_{1}}=d_{j_{1}-1}$ (because of the equality (2)). Denote by $j_{0}$ the smallest index $j<j_{1}$ such that $d_{j}=d_{j_{1}}$. Note that since $D(\makebox[4.30554pt]{\cdot})$ is increasing, we have that $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is continuous (it is upper and lower semi-continuous with Lemma 3).

For some small $\varepsilon>0$, we define $(\mathrm{bar}\boldsymbol{d},\mathrm{bar}\boldsymbol{y})$ as follows:

[TABLE]

and

[TABLE]

where $\mathrm{bar}d_{0}=0$. We first check that $(\mathrm{bar}\boldsymbol{d},\mathrm{bar}\boldsymbol{y})$ is a feasible solution of (P${}_{\mbox{{\tiny no return}}}$).

The definition of $j_{1}$ implies that $d_{j_{0}}>0$. Thus if $j_{0}=1$, we have $y_{1}>0$ and for a small enough $\varepsilon$, the vector $\mathrm{bar}\boldsymbol{y}$ satisfies constraint (ii). Otherwise, we have $\mathrm{bar}\tau(y_{j_{0}-1})+\nu(y_{j_{0}-1}-y_{j_{0}-2})=d_{j_{0}-1}<d_{j_{0}}=\mathrm{bar}\tau(y_{j_{0}})+\nu(y_{j_{0}}-y_{j_{0}-1})$. It implies that $y_{j_{0}-1}<y_{j_{0}}$ (as otherwise the equality would imply that $y_{j_{0}-1}<y_{j_{0}-2}$). Thus, for a small enough $\varepsilon$, we have $\mathrm{bar}\boldsymbol{y}$ satisfies constraint (ii). It also satisfies obviously constraint (iv).

For $j\in\{2,\ldots,j_{1}\}\cup\{j_{1}+2,\ldots,S\}$, checking $\mathrm{bar}d_{j-1}\leq\mathrm{bar}d_{j}$ is straightforward. The remaining case is $j=j_{1}+1$. A direct induction shows that $\mathrm{bar}d_{j}\leq d_{j}$ for $j\leq j_{1}-1$. Since $\mathrm{bar}\tau(y_{j_{1}})+\nu(y_{j_{1}}-y_{j_{1}-1})<\mathrm{bar}\tau(y_{j_{1}-1})+\nu(y_{j_{1}-1}-y_{j_{1}-2})$ (because $d_{j_{1}-1}=d_{j_{1}}$), for $\varepsilon$ small enough, we have $\mathrm{bar}d_{j_{1}-1}\geq\mathrm{bar}\tau(\mathrm{bar}y_{j_{1}})+\nu(\mathrm{bar}y_{j_{1}}-\mathrm{bar}y_{j_{1}-1})$. Here, we use the fact that $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is continuous. Thus $\mathrm{bar}d_{j_{1}}=\mathrm{bar}d_{j_{1}-1}$. Since we have $\mathrm{bar}d_{j_{1}-1}\leq d_{j_{1}-1}$ by the above induction, we finally obtain $\mathrm{bar}d_{j_{1}}\leq d_{j_{1}}\leq d_{j_{1}+1}=\mathrm{bar}d_{j_{1}+1}$. Therefore, $\mathrm{bar}\boldsymbol{d}$ satisfies constraint (iii).

Constraint (i) is satisfied for all $j$, except maybe for $j=j_{1}$. We have proved that $\mathrm{bar}d_{j_{1}}=\mathrm{bar}d_{j_{1}-1}$. Since $\mathrm{bar}d_{j_{1}-1}=\mathrm{bar}\tau(\mathrm{bar}y_{j})+\nu(\mathrm{bar}y_{j}-\mathrm{bar}y_{j-1})$ for some $j^{\prime}\leq j_{1}-1$, we have $\mathrm{bar}\tau(\mathrm{bar}y_{j_{1}})+\nu(\mathrm{bar}y_{j_{1}}-\mathrm{bar}y_{j_{1}-1})\leq\mathrm{bar}d_{j_{1}}=\mathrm{bar}\tau(\mathrm{bar}y_{j^{\prime}})+\nu(\mathrm{bar}y_{j^{\prime}}-\mathrm{bar}y_{j^{\prime}-1})$, and thus $\nu(\mathrm{bar}y_{j_{1}}-\mathrm{bar}y_{j_{1}-1})\leq\nu(\mathrm{bar}y_{j^{\prime}}-\mathrm{bar}y_{j^{\prime}-1})\leq\nu C$. Therefore constraint (i) is also satisfied for $j=j_{1}$.

Since the constraint (v) is clearly satisfied, $(\mathrm{bar}\boldsymbol{d},\mathrm{bar}\boldsymbol{y})$ is a feasible solution of (P${}_{\mbox{{\tiny no return}}}$).

A careful examination of the arguments used when we checked constraint (ii) shows that actually $\mathrm{bar}d_{j_{0}}<d_{j_{0}}$. The same induction as the one used we checked constraint (iii) shows that $\mathrm{bar}d_{j_{1}-1}<d_{j_{1}-1}$. We have proved that $\mathrm{bar}d_{j_{1}-1}\geq\mathrm{bar}\tau(\mathrm{bar}y_{j_{1}})+\nu(\mathrm{bar}y_{j_{1}}-\mathrm{bar}y_{j_{1}-1})$. Thus $\mathrm{bar}d_{j_{1}}=\mathrm{bar}d_{j_{1}-1}$, and $\mathrm{bar}d_{j_{1}}<d_{j_{1}}$. We have

[TABLE]

which in contradiction with the optimality assumption. This settles the case of $g^{\operatorname{ave}}(\cdot)$. The other case is dealt with similarly. ∎

Problems P${}^{\max}_{\mbox{{\tiny return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ can be written almost identically under the following form. We use infinitely many variables since there is no a priori reason to have a bounded number of departures, and there are indeed special cases for which there is no optimal solution with a finite number of departures. However, if $\pi>0$, we prove that any optimal solution of P${}^{\max}_{\mbox{{\tiny return}}}$ requires a finite number of departures, see Proposition 3. The case of P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ remains open.

[TABLE]

where $g(\makebox[4.30554pt]{\cdot})$ is either $g^{\max}(\makebox[4.30554pt]{\cdot})$ or $g^{\operatorname{ave}}(\makebox[4.30554pt]{\cdot})$. Constraints (i), (ii), (iii), (iv), and (v) have the same meaning as for the previous problems. Constraint (vi) ensures that the time between two consecutive departures of a same shuttle is not smaller than the time required for a full trip plus the time needed to load the users.

In the model (P${}_{\mbox{{\tiny return}}}$), the shuttles are not identified. Note however that their schedules can be easily be recovered: the departure times of a shuttle $s$ is of the form

[TABLE]

and the time at which the loading starts for a shuttle with departure time $d_{j}$ can be chosen to be $d_{j}-\nu(y_{j}-y_{j-1})$ (the loading starts as late as possible).

While it can be shown that problem P${}^{\max}_{\mbox{{\tiny return}}}$ always admits an optimal solution when it is feasible (see Proposition 3), we were not able to settle the case of problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$.

2.4. Computational model

We assume that the following operations take constant time:

• •

Evaluation of $D(t)$ for any $t\in[0,T]$.

• •

Integration of $D(\makebox[4.30554pt]{\cdot})$ between two values.

• •

Evaluation of $\tau(y)$ and $\mathrm{bar}\tau(y)$ for any $y\in\mathbb{R}_{+}$.

• •

Evaluation of $\sup\{y\colon\mathrm{bar}\tau(y)+\nu y\leq\alpha\}$ for any $\alpha\in\mathbb{R}_{+}$.

Note that if $D(\makebox[4.30554pt]{\cdot})$ is piecewise affine with a natural description, as it is usually the case in practice, these assumptions are easily matched. Moreover, we set as constants of the computational model the capacity $C$, the length of the period $T$, the cumulative demand $D(\makebox[4.30554pt]{\cdot})$, the loading rate $\nu$, and the return time $\pi$. The complexity functions will be expressed in terms of $S$ and the accuracy of the computed solution.

3. Main results

In the present section, we present our main findings. Many results state the existence of algorithms with a guarantee that the returned solution has a value close to the optimal value $OPT$ of the considered problem. Except for two easy results – Corollary 1 and Proposition 3 – all proofs are postponed to other sections.

We organize the results presented in that section in three subsections. The first subsection – Section 3.1 – deals with the special case where $D(\makebox[4.30554pt]{\cdot})$ is a constant function, i.e., when all users are in the loading terminal from the beginning of the period, and with returns allowed. It seems to us that these results are also interesting in themselves, because they form a natural situation for which there is a very efficient algorithm. The second subsection – Section 3.2 – deals with the general case where the shuttles are not allowed to come back, i.e., with the case covered by the problems P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$. The case where the shuttles are allowed to come back, i.e., when we deal with the problems P${}^{\max}_{\mbox{{\tiny return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$, is discussed in Section 3.3.

3.1. All users in the terminal from the beginning

In this subsection, we present results regarding the four problems when $D(t)=D(T)$ for all $t\in[0,T]$ (all users are from the beginning in the terminal). For the problems for which return is not allowed (P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$), an obvious optimal solution is given by $y_{j}^{*}=jD(T)/S$ and $d_{j}^{*}=\nu D(T)/S$ for $j\in\left\{1,\ldots,S\right\}$ and the optimal value is $\nu D(T)/S$ for both problems, provided that $D(T)\leq CS$ (otherwise, there is no feasible solution at all): the shuttles take all the same amount of users, start immediately the loading process, and have the same departure time.

The rest of the section is devoted to the results regarding the problems P${}^{\max}_{\mbox{{\tiny return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$. For the first one, there are closed-from expressions for the optimal value and an optimal solution.

Proposition 1**.**

When $D(t)=D(T)$ for all $t\in[0,T]$, the optimal value of P${}^{\max}_{\mbox{{\tiny return}}}$ is

[TABLE]

In the proof, we actually provide a closed-form expression for an optimal solution. For P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ however, there does not seem to be a closed-form expression for an optimal solution, and not even for the optimal value. There is nevertheless an efficient algorithm.

Proposition 2**.**

Suppose $\pi>0$. When $D(t)=D(T)$ for all $t\in[0,T]$, the optimal value of P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ can be computed in constant time and an optimal solution can be computed in $O(S)$.

If $\pi=0$, the optimal value is $\frac{\nu D(T)}{2S}$, and it is not too difficult to see that there is no optimal solution. In a transportation context, $\pi=0$ looks unrealistic. However, the chemical application mentioned in the introduction could be a situation where this equality could be met: as soon as the test is over for a batch, we can start the test for a new one.

3.2. When return is not allowed

We have the existence of an efficient approximation algorithm for P${}^{\max}_{\mbox{{\tiny no return}}}$. The algorithm is actually an easy binary search (Section 5.1.1).

Theorem 1**.**

Let $\rho>0$. A feasible solution $(\boldsymbol{d},\boldsymbol{y})$ of P${}^{\max}_{\mbox{{\tiny no return}}}$ – if the problem is feasible – satisfying $g^{\max}(\boldsymbol{d},\boldsymbol{y})\leq OPT+\rho$ can be computed in $O\left(S\log{\frac{1}{\rho}}\right)$.

With an additional assumption on $D(\makebox[4.30554pt]{\cdot})$, this theorem provides actually an approximation scheme.

Corollary 1**.**

If $D(\makebox[4.30554pt]{\cdot})$ is increasing, the algorithm of Theorem 1 computes in $O(S\log\frac{S}{\varepsilon})$ a $(1+\varepsilon)$-approximation for P${}^{\max}_{\mbox{{\tiny no return}}}$.

A schedule for the shuttles requires to specify $S$ real numbers. Taking an output sensitive point of view, this corollary states thus the existence of a polynomial approximation scheme in this particular case.

Proof of Corollary 1.

Suppose $D(\makebox[4.30554pt]{\cdot})$ increasing. Let $(\boldsymbol{d},\boldsymbol{y})$ be a feasible solution. According to Lemma 3, we have then $\tau(y_{j-1})=\mathrm{bar}\tau(y_{j-1})$ for every $j$ and the maximum waiting time for shuttle $j$ is at least $\tau(y_{j})+\nu(y_{j}-y_{j-1})-\tau(y_{j-1})$. Note that if $y_{j}=y_{j-1}$, this quantity is zero. Hence, the sum of the maximum waiting times over all nonempty shuttles is at least $T+\nu D(T)$ and the optimal value $OPT$ of P${}^{\max}_{\mbox{{\tiny no return}}}$ is at least $(T+\nu D(T))/S$. Setting $\rho$ to $\varepsilon(T+\nu D(T))/S$ leads to the result. ∎

For P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, there exists an efficient approximation algorithm too. The algorithm is also described later (Section 5.2.1). We already outline that this algorithm is not a binary search as in the former case, but consists in building a quite simple weighted “approximative” graph, in which a shortest path is computed.

Theorem 2**.**

Suppose that $D(\makebox[4.30554pt]{\cdot})$ admits right derivatives everywhere (denoted $D^{\prime}_{+}(t)$) and that $\inf_{t\in[0,T)}D^{\prime}_{+}(t)$ is positive. Then, for any positive integer $M$, a feasible solution $(\boldsymbol{d},\boldsymbol{y})$ of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ – if the problem is feasible – satisfying

[TABLE]

can be computed in $O\left(SM^{3}\right)$.

As for Corollary 1 above, this theorem could be interpreted as a polynomial approximation scheme by using the fact that $D(\makebox[4.30554pt]{\cdot})$ is increasing.

3.3. When return is allowed

The following proposition implies that when $\pi$ is larger than [math], any optimal solution of P${}^{\max}_{\mbox{{\tiny return}}}$ requires a finite number of nonempty departures.

Proposition 3**.**

If $\pi>0$, there exists an optimal solution of P${}^{\max}_{\mbox{{\tiny return}}}$ and the number of nonempty departures in any optimal solution is at most

[TABLE]

Proof.

A feasible solution with an infinite number of nonempty departures has an infinite objective value and it thus strictly dominated by any solution by a finite number of departures. Thus, the set of feasible solutions can be reduced to the solutions where the number of nonempty departures is finite. Since $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ is lower semicontinuous (Lemma 2), the set of feasible solutions is compact and the objective function is lower semicontinuous which leads then to existence of an optimal solution.

The schedule consisting in making the shuttles wait until time $T$, loading them at full capacity (except maybe for the last departure), and making them leave as soon as the loading is completed provides a feasible solution of P${}^{\max}_{\mbox{{\tiny return}}}$ with a value $T+\nu D(T)/S+(\lceil D(T)/(CS)\rceil-1)\pi$.

Consider an optimal solution of P${}^{\max}_{\mbox{{\tiny return}}}$ and denote by $k$ the number of departures after time $T$. The users in the last shuttle to leave have waited at least $(k/S-1)\pi$. We have thus

[TABLE]

which implies that

[TABLE]

Before time $T$, the number of departures is at most $\lceil T/\pi\rceil S$. ∎

The next theorem states that there exists an algorithm computing arbitrarily good feasible solutions for P${}^{\max}_{\mbox{{\tiny return}}}$ within reasonable computational times when $S$ is small. As for Section 3.2, this algorithm is described later in the paper (Section 6). It is based on the computation of a shortest path in an “approximative” graph, as for Theorem 2. It also uses Proposition 3 in a crucial way (actually a slight variation of it: Lemma 10).

Theorem 3**.**

Suppose that $D(\makebox[4.30554pt]{\cdot})$ admits right derivatives everywhere, $\pi$ is positive, and $\inf_{t\in[0,T)}D^{\prime}_{+}(t)$ is positive. Then, for any positive integer $M$, a feasible solution $(\boldsymbol{d},\boldsymbol{y})$ of P${}^{\max}_{\mbox{{\tiny return}}}$ satisfying

[TABLE]

can be computed in $O\left(\beta^{3S}M^{3S+2}\right)$, where $\beta$ depends only on the constants of the computational model.

As above, the theorem actually ensures that the algorithm is an approximation scheme since we can bound from below $OPT$ using only the input values. If $S$ is considered as constant, this becomes even a polynomial approximation scheme.

We do not know whether there is a counterpart to Proposition 3 for problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$. If such a counterpart would exist, then almost the same technique as the one used in Section 6 would lead to a theorem similar to Theorem 3 for P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$. The existence of such a theorem remains thus open.

4. All users in the terminal from the beginning

Consider the case where all users are in the loading terminal from the beginning. To ease the reading, and for the present section only, we use $D$ to denote the quantity $D(T)$.

We treat first the case of problem P${}^{\max}_{\mbox{{\tiny return}}}$.

Proof of Proposition 1.

For P${}^{\max}_{\mbox{{\tiny return}}}$, when $S=1$, an optimal solution is obtained by loading at full capacity the shuttle for each departure (except maybe for the last departure for which the shuttle load is $D-C\lfloor D/C\rfloor$) and by making the shuttle leave immediately after each loading process. The optimal value is then $\nu D+(\lceil D/C\rceil-1)\pi$. When $S>1$, consider the problem Q defined as the problem P${}^{\max}_{\mbox{{\tiny return}}}$ without the constraint that the shuttles do not overtake (constraint (iii) in (P${}_{\mbox{{\tiny return}}}$)). The optimal value of Q provides a lower bound of the optimal value of P${}^{\max}_{\mbox{{\tiny return}}}$. Since there is no constraint linking the different shuttles, problem Q can be solved separately for each shuttle $s$ with a demand $D_{s}$ to carry, such that $\sum_{s}D_{s}=D$. The optimal solutions of Q are thus obtained from the optimal solutions of

[TABLE]

The solution given by $D_{s}=D/S$ for all $s$ is clearly optimal (and it is actually the unique optimal solution when $\nu>0$). Hence, there is an optimal solution for Q in which all shuttles have the same departure times and, for each travel, carry the same amount of users. Its value is $\nu D/S+(\lceil D/(CS)\rceil-1)\pi$. Since the shuttles do not overtake in this optimal solution of Q, it is actually a feasible solution for the problem P${}^{\max}_{\mbox{{\tiny return}}}$, and thus an optimal solution for this latter problem (its value being equal to a lower bound). ∎

The rest of this section is devoted to the proof of Proposition 2, which ensures the existence of an efficient algorithm solving problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ when $D(t)=D$ for all $t\in[0,T]$. We start by considering the special case of problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ when $S=1$. In such a case, it is always beneficial to define $d_{j}=(j-1)\pi+\nu y_{j}$. Assuming that the $y_{j}$’s are given, it provides a feasible solution since $\mathrm{bar}\tau(y)=0$ for all $y\in[0,D]$. The objective function of P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ becomes thus

[TABLE]

where $x_{j}=y_{j}-y_{j-1}$. Solving P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ when $S=1$ reduces thus to solving

[TABLE]

which is a convex program (with infinitely many variables). We will show that there is always an optimal solution of $P(D)$ with a finite support. Then, we will solve $P_{0}(D)$, defined as the program $P(D)$ with the additional constraint $|\{j\colon x_{j}\neq 0\}|<+\infty$, with the help of the Karush-Kuhn-Tucker conditions (that do not apply otherwise).

Lemma 4**.**

Suppose that $\pi>0$. Then $P_{0}(D)$ has an optimal solution and it is necessarily of the form

[TABLE]

with $a\in\mathbb{Z}_{+}$ such that $a\leq\frac{D}{C}$ and where

[TABLE]

Proof.

Consider the following program

[TABLE]

Note that $P_{0}^{n}(D)$ is actually $P_{0}(D)$ with the additional constraint that $\sup\{j\colon x_{j}\neq 0\}\leq n$.

For $n<D/C$, $P_{0}^{n}(D)$ has no feasible solutions, and for $n\geq D/C$, the set of feasible solutions is nonempty. In this case, by compactness and continuity of the objective function, $P_{0}^{n}(D)$ has an optimal solution $\boldsymbol{x}^{*}$. We necessarily have $x_{j}^{*}\geq x_{j+1}^{*}$ for every $j\in\left\{1,\ldots,n-1\right\}$, otherwise, exchanging the two values would strictly decrease the objective function. Let $a$ be the largest index $j$ such that $x_{j}^{*}=C$, with the convention that $a=0$ if there are no such index $j$, and let $b+1$ be the smallest index $j$ such that $x_{j}^{*}=0$, with the convention that $b=n$ if there is no such index $j$.

The constraints being all affine, the Karush-Kuhn-Tucker conditions apply. There is thus a real number $\lambda\in\mathbb{R}$ and two collections $\boldsymbol{\mu},\boldsymbol{\omega}\in\mathbb{R}_{+}^{n}$ such that for every $j\in\{1,\ldots,n\}$ we have

[TABLE]

Summing this equality from $j=a+1$ to $j=b$ and noting that $\mu_{j}=\omega_{j}=0$ and $\sum_{j=a+1}^{b}x_{j}^{*}=D-aC$ by definition of $a$ and $b$ provide an expression of $\lambda$ in terms of $a$ and $b$. Replacing $\lambda$ by this expression in the same equality leads to

[TABLE]

Using this equality for $j=b$ gives the following equation.

[TABLE]

Equation (3) specialized for $j=b+1$ gives

[TABLE]

These two inequalities together – treated as conditions on a second order polynomial in $b-a$ – imply the necessary condition

[TABLE]

which imposes a unique integer value for $b-a$ and $b$ takes a unique value $\theta(a)$ for each value of $a$. We have proved that any optimal solution of $P_{0}^{n}(D)$ is of this form. Now, note that by definition of $a$, we necessarily have $a\leq\lfloor D/C\rfloor$. It means that there are only finitely many optimal solutions of the $P_{0}^{n}(D)$’s when $n$ goes to infinity. Since the set of feasible solutions of the $P_{0}^{n}(D)$’s is nondecreasing when $n$ goes to infinity, it means actually that there exists an $n_{0}$ such that any optimal solution of $P_{0}^{n}(D)$ for $n\geq n_{0}$ is an optimal solution of $P_{0}^{n_{0}}(D)$. Moreover, any feasible solution of $P_{0}(D)$ is a feasible solution of $P_{0}^{n}(D)$ for some $n\geq n_{0}$, and thus is dominated by the optimal solutions of $P_{0}^{n_{0}}(D)$. These latter are thus the optimal solutions of $P_{0}(D)$. ∎

Let $v(D)$ and $v_{0}(D)$ be the optimal values of respectively $P(D)$ and $P_{0}(D)$. Note that $v(D)\leq v_{0}(D)$.

Lemma 5**.**

If $\pi>0$, we have

[TABLE]

for every $\varepsilon\in(0,D]$.

Proof.

Let $\varepsilon\in(0,D]$. Consider a feasible solution $\boldsymbol{x}$ of $P(D)$. Let $N_{\varepsilon}\in\mathbb{Z}_{+}$ be such that $\sum_{j=N_{\varepsilon}+1}^{+\infty}x_{j}<\varepsilon$. Define inductively

[TABLE]

This $\boldsymbol{x}^{\prime}$ is a feasible solution of $P_{0}(D-\varepsilon)$. Since $x_{j}^{\prime}\leq x_{j}$ for all $j$, the value given by $\boldsymbol{x}^{\prime}$ to the objective value of $P_{0}(D-\varepsilon)$ is nonlarger that the value obtained by $\boldsymbol{x}$ for $P(D)$. The inequality follows. ∎

Proof of Proposition 2.

Let us deal with the case $S=1$. Using the fact that $P_{0}(D)$ is a convex program, we easily get that $v_{0}(\makebox[4.30554pt]{\cdot})$ is a convex function. It is thus continuous on $(0,+\infty)$, and since $v_{0}(0)=0$, we have that $v_{0}(\makebox[4.30554pt]{\cdot})$ is continuous everywhere on $[0,+\infty)$. Making $\varepsilon$ tend toward [math] in Lemma 5 and the inequality $v(D)\leq v_{0}(D)$ imply that $v_{0}(D)=v(D)$. Since any feasible solution of $P_{0}(D)$ is a feasible solution of $P(D)$ with the same value for the objective function, every optimal solution of $P_{0}(D)$ is an optimal solution of $P(D)$. An algorithm computing an optimal solution of $P(D)$ can then be derived from Lemma 4: we just have to try all the finitely many possible values for $a$. The proof for any value of $S$ will be obtained by showing that an optimal solution in this case consists just in replicating optimal solutions for the one-shuttle case.

When $S>1$, consider the problem Q defined as the problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$ without the constraint that the shuttles do not overtake (constraint (iii) in (P${}_{\mbox{{\tiny return}}}$)). The optimal value of Q provides a lower bound of the optimal value of P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$. Since there is no constraint linking the different shuttles, problem Q can be solved separately for each shuttle $s$ with a demand $D_{s}$ to carry, such that $\sum_{s}D_{s}=D$. The optimal solutions of Q are thus obtained from the optimal solutions of

[TABLE]

The fact that $P(D)$ is a convex program implies that the map $v(\makebox[4.30554pt]{\cdot})$ is convex. As a consequence, the solution $D_{s}=D/S$ for all $s$ is an optimal solution of the previous program. Hence, there is an optimal solution for Q in which all shuttles have the same departure times and, for each travel, carry the same amount of users. Since the shuttles do not overtake in this optimal solution of Q, it is actually a feasible solution for the problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$, and thus an optimal solution for this latter problem (its value being equal to a lower bound). ∎

5. When return is not allowed

5.1. Minimizing the maximum waiting time

5.1.1. The algorithm

If $CS<D(T)$, there is no feasible solution. We can thus assume that $CS\geq D(T)$. The algorithm is a binary search starting with the values $h^{+}=T+\nu D(T)$ and $h^{-}=0$ which are respectively upper and lower bounds of the optimal value. While the gap $h^{+}-h^{-}$ is larger than $\rho$, we consider the tentative value $h=\frac{h^{+}+h^{-}}{2}$ and the system

[TABLE]

where $\mathcal{S}_{j}^{h}=\left\{y\in\mathbb{R}_{+}\colon y\leq C+y_{j-1},\mathrm{bar}{\tau}(y)+\nu(y-y_{j-1})-\tau(y_{j-1})\leq h,y\leq D(T)\right\}$. Each iteration of the binary search consists in deciding whether (Sh) has a feasible solution or not, and it can be done in $O(S)$ by computing the values of the $y_{j}$’s and the $d_{j}$’s iteratively (here we use in particular the computational assumptions on $D$). As we are going to prove, (Sh) has a feasible solution if and only if the problem has a feasible solution with a value of the objective function at most $h$. If (Sh) has a feasible solution, we update thus the value of $h^{+}$ with the current value of $h$, otherwise, we update $h^{-}$ with $h$. When $h^{+}-h^{-}\leq\rho$, the solution of program $(\text{S}_{h^{+}})$ is feasible for P${}^{\max}_{\mbox{{\tiny no return}}}$ and its value $h^{+}$ is at most at $\rho$ from the optimal value.

5.1.2. Proof of Theorem 1

For any fixed $h$, P${}^{\max}_{\mbox{{\tiny no return}}}$ has a feasible solution with a value of the objective function at most $h$ if and only if the following system has a feasible solution.

[TABLE]

We claim that (Qh) has a feasible solution if and only if (Sh) has one. Once this equivalence is established, the correctness of the binary search described above is almost immediate using Claim 1.

Let $(\boldsymbol{d},\boldsymbol{y})$ be a feasible solution of (Sh). We use without further mention that $\mathcal{S}_{j}^{h}$ is closed. It satisfies the constraints $\textup{(Qi)},\textup{(Qii)}$, and (Qv). We have $y_{j-1}\leq C+y_{j-1}$ and $y_{j-1}\leq D(T)$. Since $\mathrm{bar}{\tau}(y)\leq\tau(y)$ for all $y$, we also have $\mathrm{bar}{\tau}(y_{j-1})+\nu(y_{j-1}-y_{j-1})-\tau(y_{j-1})\leq h$. It means that $y_{j-1}$ belongs to $\mathcal{S}_{j}^{h}$, and thus $y_{j-1}\leq y_{j}$. Hence, $(\boldsymbol{d},\boldsymbol{y})$ satisfies also constraint (Qiii). Since $\mathrm{bar}{\tau}(y_{j})+\nu(y_{j}-y_{j-1})-\tau(y_{j-1})\leq h$, the solution also satisfies constraint (Qvi) and since $\tau(\makebox[4.30554pt]{\cdot})$ is nondecreasing, it satisfies constraint (Qiv). Therefore, it is a feasible solution of (Qh) and the existence of a feasible solution of (Sh) implies the existence of a feasible solution of (Qh).

For the converse implication, suppose that (Qh) admits a feasible solution, and consider the optimization problem consisting in maximizing $\sum_{j=1}^{S}y_{j}$ over its feasible solutions. These feasible solutions form a compact set of $\mathbb{R}_{+}^{S}$ since it is obviously bounded and since the semicontinuities of $\tau(\makebox[4.30554pt]{\cdot})$ and $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ imply that it is closed. There is thus an optimal solution $(\boldsymbol{d}^{*},\boldsymbol{y}^{*})$ to that optimization problem. Suppose for a contradiction that there is a $j$ such that $y_{j}^{*}<\sup\mathcal{S}_{j}^{h}$. Denote $j_{0}$ the largest such index. Let us slightly increase $y_{j_{0}}^{*}$, while letting the other $y_{j}^{*}$ untouched. Redefine $d_{j}^{*}$ to be $h+\tau(y_{j-1}^{*})$ for all $j\geq j_{0}$. The pair $(\boldsymbol{d}^{*},\boldsymbol{y}^{*})$ remains feasible for (Qh) (we use here the fact that $\sup\mathcal{S}_{j}^{h}$ is nondecreasing with $j$), while increasing the quantity $\sum_{j=1}^{S}y_{j}^{*}$, which is a contradiction with the optimality assumption. Thus, we have $y_{j}^{*}=\sup\mathcal{S}_{j}^{h}$ for all $j$ and $d_{j}^{*}:=h+\tau(y_{j-1}^{*})$ for all $j$ provides a feasible solution for (Sh). ∎

5.2. Minimizing the average waiting time

5.2.1. The algorithm

The following map will be useful in the description of the algorithm.

[TABLE]

Define the directed graph $\mathcal{G}=(\mathcal{V},\mathcal{A})$ by

[TABLE]

where we use the following notations:

[TABLE]

Set for each arc a=\big{(}(z,r),(z^{\prime},r^{\prime})\big{)} a weight w(a)=f^{\operatorname{ave}}\big{(}\mathrm{bar}\tau(r^{\prime})+\nu(z^{\prime}-\eta),r+\eta,r^{\prime}\big{)}.

If $CS<D(T)$, there is no feasible solution. We can thus assume that $CS\geq D(T)$. The algorithm consists first in computing a path $\tilde{p}$ minimizing $\sum_{a\in\mathcal{A}(p)}w(a)$, among all paths $p$ with at most $S$ arcs starting at $(0,0)\in\mathcal{V}$ and ending at a vertex $(z,r)$ with $r=R\eta$. Such paths exist, see Lemma 6 below. The computation of $\tilde{p}$ can be done in $O(S|\mathcal{A}|)$ via dynamic programming. Let the vertex sequence of $\tilde{p}$ be \big{(}(z_{0},r_{0}),(z_{1},r_{1}),\ldots,(z_{n},r_{n})\big{)}. The algorithm consists then in defining recursively

[TABLE]

and

[TABLE]

and outputting the pair $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$. The construction of the graph is sketched on Figure 2.

As it will be shown below, this $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is a feasible solution of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ providing a value to the objective function within a $O\left(\frac{S^{2}}{M}\right)$ gap to the optimal value.

5.2.2. Proof of Theorem 2

We provide three lemmas, which are proved in a separate section at the very end of the paper to ease the reading. Theorem 2 results immediately from their combination.

In the proofs of the lemmas, we assume that $M$ is large enough so that $\eta<D(T)$. Since in Theorem 2, $M$ appears in ‘big O’ formulas, it is a valid assumption. Anyway it is what is sought in practice: the larger $M$, the larger the accuracy of the solution. An $\eta$ of same order of magnitude of $D(T)$ would be useless.

Lemma 6**.**

For every optimal solution $(\boldsymbol{d}^{*},\boldsymbol{y}^{*})$, there is a path $p$ with at most $S$ arcs starting at $(0,0)\in\mathcal{V}$ and ending at a vertex $(z,r)$ with $r=R\eta$ and such that

[TABLE]

Lemma 7**.**

The pair $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is a feasible solution of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$.

Lemma 8**.**

The following inequality holds:

[TABLE]

5.3. When the demand function is a step function

Better complexity results can be obtained when the demand is a step function. A step function is a function that can be written as a finite linear combination of indicator functions of intervals. The assumption on $D(\makebox[4.30554pt]{\cdot})$ being a step function means that the users arrive only on a finite number of instants. As it has already been noted, the assumption $\nu=0$ is equivalent to the assumption that every user boards a shuttle as soon as he arrives in the terminal.

Proposition 4**.**

Assume that $D(\makebox[4.30554pt]{\cdot})$ is a step function defined with $K$ discontinuities, supposed to be part of the input. Suppose moreover that $\nu=0$. Then for each of P${}^{\max}_{\mbox{{\tiny no return}}}$ and of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, there is an algorithm computing an optimal solution in $O(K^{2}S)$.

It turns out that when $C$ and the values taken by $D(\makebox[4.30554pt]{\cdot})$ are integer, the loads of the shuttles in the optimal solution returned by the algorithm are also integer. We cover thus the case where the users are atoms.

The algorithm

We provide only the algorithm for P${}^{\max}_{\mbox{{\tiny no return}}}$, the other case can be dealt similarly. Let $t_{1}<\cdots<t_{K}$ be the $K$ discontinuities. Define the directed graph $\mathcal{G}=(\mathcal{V},\mathcal{A})$ by

[TABLE]

where $Q=\lfloor D(T)/C\rfloor$. Note that the vertex set is a finite subset of $\mathbb{R}_{+}$. Set for each arc $a=(y,y^{\prime})$ a weight $w(a)=\mathrm{bar}\tau(y^{\prime})-\tau(y)$. We consider the two vertices [math] and $D(T)$ (obtained with $k=K$ and $q=0$).

If $CS<D(T)$, there is no feasible solution. We can thus assume that $CS\geq D(T)$. The algorithm consists first in computing a [math]-$D(T)$ path $\tilde{p}$ with $S$ arcs minimizing $\max_{a\in\mathcal{A}(\tilde{p})}w(a)$. Within the proof of Proposition 4 below, we show that from any feasible solution we can build a [math]-$D(T)$ path with $S$ arcs in $\mathcal{G}$. Thus, when the problem is feasible, such paths exist in $\mathcal{G}$. The computation of $\tilde{p}$ can be done in $O(S|\mathcal{A}|)$ via dynamic programming. Let the vertex sequence of $\tilde{p}$ be $(\tilde{y}_{0},\tilde{y}_{1},\ldots,\tilde{y}_{S})$. The end of the algorithm consists in defining $\tilde{d}_{j}=\mathrm{bar}\tau(\tilde{y}_{j})$ for all $j\in\{1,\ldots,S\}$ and outputting the pair $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$.

As it will be shown below, this $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is an optimal solution of P${}^{\max}_{\mbox{{\tiny no return}}}$.

Proof of Proposition 4

According to Claim 1, we replace the objective function of P${}^{\max}_{\mbox{{\tiny no return}}}$ by $\max_{j\in\left\{1,\ldots,S\right\}}(d_{j}-\tau(y_{j-1}))$. It can easily be checked that $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is feasible for P${}^{\max}_{\mbox{{\tiny no return}}}$. It provides a value $\max_{j\in\left\{1,\ldots,S\right\}}(\mathrm{bar}\tau(\tilde{y}_{j})-\tau(\tilde{y}_{j-1}))$ for P${}^{\max}_{\mbox{{\tiny no return}}}$ (with the alternative objective function), and this value coincides with $\max_{a\in\mathcal{A}(\tilde{p})}w(a)$. The path $\tilde{p}$ describes therefore a solution of P${}^{\max}_{\mbox{{\tiny no return}}}$ with a value equal to $\max_{a\in\mathcal{A}(\tilde{p})}w(a)$.

Conversely, let $(\boldsymbol{d},\boldsymbol{y})$ be any feasible solution of P${}^{\max}_{\mbox{{\tiny no return}}}$. Let $\mathrm{bar}\boldsymbol{y}$ be the sequence defined by $\mathrm{bar}y_{j}=\min\{y\in\mathcal{V}\colon y\geq y_{j}\}$. On the one hand, we have $\mathrm{bar}y_{j-1}\leq\mathrm{bar}y_{j}$ because $y_{j-1}\leq y_{j}$. On the other hand, we have $\mathrm{bar}y_{j-1}+C\geq y_{j-1}+C\geq y_{j}$. If $\mathrm{bar}y_{j-1}+C\in\mathcal{V}$, we have $\mathrm{bar}y_{j-1}+C\geq\mathrm{bar}y_{j}$ by definition of $\mathrm{bar}y_{j}$. If $\mathrm{bar}y_{j-1}+C\notin\mathcal{V}$, then $\mathrm{bar}y_{j-1}+C>D(T)\geq\mathrm{bar}y_{j}$ since $D(T)\in\mathcal{V}$. Thus, $(\mathrm{bar}y_{j-1},\mathrm{bar}y_{j})\in\mathcal{A}$ for all $j\in\{1,\ldots,S\}$. We have $\mathrm{bar}y_{0}=0$ and $\mathrm{bar}y_{S}=D(T)$ and the sequence $\mathrm{bar}\boldsymbol{y}$ is a [math]-$D(T)$ path $p$ with $S$ arcs.

Second, we prove that $\mathrm{bar}\tau(\mathrm{bar}y_{j})-\tau(\mathrm{bar}y_{j-1})\leq d_{j}-\tau(y_{j-1})$ as follows. There exists a unique $k$ such that $D(t_{k})<y_{j}\leq D(t_{k+1})$. By definition of $D(\makebox[4.30554pt]{\cdot})$, we have $D(t)=D(t_{k})$ for all $t\in[t_{k},t_{k+1})$, and thus $\mathrm{bar}\tau(y_{j})=t_{k+1}$. Since $D(t_{k+1})\in\mathcal{V}$, we have $\mathrm{bar}y_{j}\leq D(t_{k+1})$ by definition of $\mathrm{bar}y_{j}$, and hence $\mathrm{bar}\tau(\mathrm{bar}y_{j})\leq t_{k+1}$ (directly by definition of $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$). Therefore, $\mathrm{bar}\tau(\mathrm{bar}y_{j})-\tau(\mathrm{bar}y_{j-1})\leq\mathrm{bar}\tau(y_{j})-\tau(y_{j-1})\leq d_{j}-\tau(y_{j-1})$ (where we use the fact that $\tau(\makebox[4.30554pt]{\cdot})$ is nondecreasing).

Finally, we have \max_{a\in\mathcal{A}(p)}w(a)\leq\max_{j\in\{1,\ldots,S\}}\big{(}d_{j}-\tau(y_{j-1})\big{)}. As the path $\tilde{p}$ is optimal, $\max_{a\in\mathcal{A}(\tilde{p})}w(a)$ is a lower bound on the value taken by the (alternative) objective function on $(\boldsymbol{d},\boldsymbol{y})$. ∎

6. When return is allowed

The algorithm

The following map will be useful:

[TABLE]

We introduce the following two sets

[TABLE]

where

[TABLE]

Define the directed graph $\mathcal{G}=(\mathcal{V},\mathcal{A})$ by

[TABLE]

where

[TABLE]

We adopt the convention $D(t)=D(T)$ when $t\geq T$ and we define $r_{0}=r_{1}-z$. Set for each arc a=\big{(}(z,\boldsymbol{q},\boldsymbol{r}),(z^{\prime},\boldsymbol{q}^{\prime},\boldsymbol{r}^{\prime})\big{)} a weight $w(a)=\max_{k\in\left\{1,\ldots,S\right\}}f^{\max}(q_{k}^{\prime}-\eta,r_{k-1}^{\prime}+\eta,r^{\prime}_{k})$ where $r^{\prime}_{0}=r^{\prime}_{1}-z^{\prime}$.

The algorithm consists first in computing a path $\tilde{p}$ minimizing $\max_{a\in A(p)}w(a)$ among all paths $p$ starting at $(0,\boldsymbol{0},\boldsymbol{0})\in\mathcal{V}$ (the ‘all zero’ vector) and ending at a vertex $(z,\boldsymbol{q},\boldsymbol{r})$ with $r_{S}=R\eta$. Such paths exist, see Lemma 9 below. It can be done in $O(|\mathcal{V}||\mathcal{A}|)$ via dynamic programming. Let the vertex sequence of $\tilde{p}$ be \big{(}(0,\boldsymbol{0},\boldsymbol{0}),(z_{0},\boldsymbol{q}^{0},\boldsymbol{r}^{0}),\ldots,(z_{n},\boldsymbol{q}^{n},\boldsymbol{r}^{n})\big{)}. The vector $\boldsymbol{r}^{i}$ models the cumulative loads of the $S$ shuttles when they perform their $i$th departure. The vector $\boldsymbol{q}^{i}$ models the times at which the loading of the $S$ shuttles starts when they perform their $i$th departure. These quantities are computed only approximatively (with an accuracy $\eta$).

The algorithm consists then in defining recursively for all $j=iS+k$ with $i=0,\ldots,n$ and $k=1,\ldots,S$

[TABLE]

where $\tilde{y}_{0}=0$, $r^{i}_{0}=r^{i-1}_{S}$, $r^{0}_{0}=0$, and $\tilde{\gamma}=(1+2\nu+1/\alpha)$. For $j=(n+1)S+1,\ldots,N$

[TABLE]

and outputting the pair $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$.

As it will be stated below, this $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is a feasible solution of P${}^{\max}_{\mbox{{\tiny return}}}$ providing a value to the objective function $g^{\max}(\makebox[4.30554pt]{\cdot})$ within a $O\left(\frac{S^{2}}{M}\right)$ gap to the optimal value.

Proof of Theorem 3

We provide four lemmas. The proof of Lemma 10 is almost the one of Proposition 3 and the proofs of the three others follow the same scheme as the ones of Lemmas  6, 7, and 8. They are thus omitted. Theorem 3 results immediately from their combination.

Lemma 9**.**

For every optimal solution $(\boldsymbol{d}^{*},\boldsymbol{y}^{*})$, there is a path $p$ starting at $(0,\boldsymbol{0},\boldsymbol{0})\in\mathcal{V}$ and ending at a vertex $(z,\boldsymbol{q},\boldsymbol{r})$ with $r_{S}=R\eta$ and such that

[TABLE]

Lemma 10**.**

There is an optimal solution of P${}^{\max}_{\mbox{{\tiny return}}}$ for which $T^{+}$ is an upper bound on the loading time of the last departure.

Lemma 11**.**

The pair $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is a feasible solution of P${}^{\max}_{\mbox{{\tiny return}}}$.

Lemma 12**.**

The following inequality holds:

[TABLE]

7. Experimental results

In this section, we test the performance of the algorithms described in previous sections for problems P${}^{\max}_{\mbox{{\tiny no return}}}$, P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, and P${}^{\max}_{\mbox{{\tiny return}}}$. As explained in Section 3, we do not have such an algorithm for problem P${}^{\operatorname{ave}}_{\mbox{{\tiny return}}}$.

7.1. Data

Our experiments are based on real data provided by our partner Eurotunnel. They are related to the transportation of freight trucks between France and Great Britain. Some parameters are fixed as constants of the problems and do not vary from an instance to another. For the constants $C,T,\nu,\pi$ of our problem, we take the real values used in practice by the company:

[TABLE]

(The value taken for $\pi$ is actually the duration of a trip going from France to Great-Britain, and not of the round trip, which lasts approximatively twice this quantity.)

Two functions $D(\makebox[4.30554pt]{\cdot})$ are used. The first one (“1P”) is a piecewise affine map obtained by averaging the real demand over several days. It turns out that this function has a peak period in the morning. The second function (“2P”), also piecewise affine, is obtained from the first by artificially adding a second peak period in the evening. In both cases, $D(\makebox[4.30554pt]{\cdot})$ is increasing and $D(T)=2016$. For problems P${}^{\max}_{\mbox{{\tiny no return}}}$ and P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, we consider $S\in[100,250]$ since the number of shuttle trips in every direction is within this range for a typical day.

The numerical experiments are performed on a Macbook Pro of 2014 with four 2.2 Ghz processors and 16 Gb of ram.

7.2. Results

The problems P${}^{\max}_{\mbox{{\tiny no return}}}$, P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, and P${}^{\max}_{\mbox{{\tiny return}}}$ are solved with algorithms described in this article. The results are summarized in the following tables.

Table 1 gives the numerical results for problem P${}^{\max}_{\mbox{{\tiny no return}}}$. The next column is the number of shuttles $S$ in the fleet. The third column provides the parameter $\varepsilon$ of the algorithm, which is an a priori upper bound on the optimality gap (Corollary 1). The two next columns give respectively the lower bound and the upper bound (value of the feasible solution returned by the algorithm), both expressed in minutes. The next column is the optimality gap. The last column provides the CPU time spent solving the problem.

Table 2 gives the numerical results for problem P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$. The columns are the same as for Table 1 except the third one which provides here the parameter $M$ of the algorithm. We know from Theorem 2 that the gap between the upper bound and the lower bound converges asymptotically to 0 when $M$ goes to infinity. We tried $M=32$ and $M=128$.

Table 3 gives the numerical results for problem P${}^{\max}_{\mbox{{\tiny return}}}$. The columns are the same as for Table 2. Since the computation time was prohibitively long as soon as $S\geq 2$, we made experiments for $S=1$. To get realistic waiting times for the users, we divided the demand functions by 3.5 leading to (“1P∗”) and (“2P∗”). Again, we know from Theorem 3 that for large $M$, we will be close to the optimal solution and we tried $M=16$ and $M=32$.

7.3. Comments

In Table 1, the results for problem P${}^{\max}_{\mbox{{\tiny no return}}}$ are extremely convincing, the optimal solutions were found almost immediately. In Table 2, the algorithm for problem P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ was able to find provable good solutions within reasonable computation times. We may note that increasing $M$ after some threshold does not seem to improve the quality of the return solution. This was confirmed by other experiments not shown here. It may indicate that the algorithm could be used efficiently in practice. In Table 3, the same holds for P${}^{\max}_{\mbox{{\tiny return}}}$ once we have accepted to work with one shuttle. Finding an efficient algorithm with at least two shuttles seems to remain a challenging task.

8. Proofs of Lemmas of Section 5.2

Claim 3**.**

We have $r_{j}\leq\tilde{y}_{j}\leq r_{j}+\eta$ for $j=0,\ldots,n$.

Proof.

We have $\tilde{y}_{j}\leq r_{j}+\eta$ by definition. Using $r_{j}-r_{j-1}\leq M\eta$ in a feasible path, a direct induction shows that $\tilde{y}_{j}\geq r_{j}$ for $j=0,\ldots,n$. ∎

Claim 4**.**

Suppose that $\alpha>0$. Then for all $y\in[0,D(T)]$ and $\delta\in[0,D(T)-y]$, we have $\mathrm{bar}\tau(y+\delta)\leq\mathrm{bar}\tau(y)+\delta/\alpha$ and $\tau(y+\delta)\leq\tau(y)+\delta/\alpha$.

Proof.

Diewert [7] extended the Mean Value Theorem to semicontinuous functions. According to his result, for any $0\leq a\leq b\leq T$, there exists $c\in[a,b)$ such that

[TABLE]

Since

[TABLE]

we have $D(a)+\alpha(b-a)\leq D(b)$ for any $0\leq a\leq b\leq T$. With $a=\mathrm{bar}\tau(y)$ and $b=\mathrm{bar}\tau(y)+\delta/\alpha$, we get $y+\delta\leq D(\mathrm{bar}\tau(y))+\delta\leq D(\mathrm{bar}\tau(y)+\delta/\alpha)$ (the first inequality is given by Lemma 1). By definition of $\mathrm{bar}\tau$, we have $\mathrm{bar}\tau(y+\delta)\leq\mathrm{bar}\tau(y)+\delta/\alpha$. The second inequality is proved along the same lines. ∎

Proof of Lemma 6.

Let $(\boldsymbol{d}^{*},\boldsymbol{y}^{*})$ be an optimal solution of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$ such that $d_{j}^{*}=\mathrm{bar}\tau(y_{j}^{*})+\nu(y_{j}^{*}-y_{j-1}^{*})$ for all $j\in\{1,\ldots,S\}$ (Claim 2). Consider the sequence $\lfloor y^{*}_{1}/\eta\rfloor\eta,\ldots,\lfloor y^{*}_{S}/\eta\rfloor\eta$ and remove the repetitions. Since the sequence is nondecreasing, we obtain an increasing sequence $\boldsymbol{r}=r_{1},\ldots,r_{n}$. We introduce $\sigma\colon\{1,\ldots,n\}\rightarrow\{1,\ldots,S\}$ with $\sigma(j)$ being the smallest index such that $r_{j}=\lfloor y_{\sigma(j)}^{*}/\eta\rfloor\eta$. We then define $z_{j}=r_{j}-r_{j-1}$ for $j\in\left\{1,\ldots,n\right\}$, with $r_{0}=0$. We prove that the sequence $(z_{j},r_{j})_{j\in\left\{1,\ldots,n\right\}}$ provides a feasible path from the vertex $(0,0)$ to $(z_{n},r_{n})$ in $\mathcal{G}$. First note that $r_{n}=R\eta$ since $y_{S}^{*}=D(T)$ and that $z_{j}>0$. For all $j\in\left\{1,\ldots,n\right\}$, we have $z_{j}=r_{j}-r_{j-1}=\left(\lfloor y_{\sigma(j)}^{*}/\eta\rfloor-\lfloor y_{\sigma(j)-1}^{*}/\eta\rfloor+\lfloor y_{\sigma(j)-1}^{*}/\eta\rfloor-\lfloor y_{\sigma(j-1)}^{*}/\eta\rfloor\right)\eta<M\eta+\eta$, since $\lfloor y_{\sigma(j)-1}^{*}/\eta\rfloor=\lfloor y_{\sigma(j-1)}^{*}/\eta\rfloor$ and $y^{*}_{\sigma(j)}-y^{*}_{\sigma(j)-1}\leq C$. Thus $z_{j}\leq M\eta$. Moreover by definition, $r_{j}\leq R\eta$. Therefore $(z_{j},r_{j})\in\mathcal{V}$ for all $j\in\left\{1,\ldots,n\right\}$. Let us now prove that $((z_{j-1},r_{j-1}),(z_{j},r_{j}))\in\mathcal{A}$ for all $j\in\left\{2,\ldots,n\right\}$. By definition, $z_{j}+r_{j-1}=r_{j}$. Note that because of the definition of $r_{j}$, we have $r_{j}\leq y_{\sigma(j)}^{*}\leq y_{\sigma_{(j+1)}-1}^{*}<r_{j}+\eta$. Combining these inequalities for all $j$ with Claim 4 leads to

[TABLE]

The sequence $(z_{j},r_{j})_{j\in\left\{1,\ldots,n\right\}}$ is then a feasible path $p$ from the vertex $(0,0)$ to $(z_{n},r_{n})$ in $\mathcal{G}$, with at most $S$ arcs. The only thing that remains to be checked in that the claimed inequality holds.

We have f^{\operatorname{ave}}(d_{\sigma(j)}^{*},y_{\sigma(j)-1}^{*},y_{\sigma(j)}^{*})\geq f\big{(}\mathrm{bar}\tau(r_{j})+\nu(z_{j}-\eta),r_{j-1}+\eta,r_{j}\big{)} for all $j\in\left\{1,\ldots,n\right\}$ since $f^{\operatorname{ave}}(\makebox[4.30554pt]{\cdot})$ is nonincreasing in the second term and nondecreasing in the first and third terms. Thus,

[TABLE]

Since this inequality holds for any optimal solution of P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$, we get the conclusion. ∎

Proof of Lemma 7.

We are going to check that $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is feasible for P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$.

For $j=1,\ldots,n$, we have $\tilde{y}_{j}-\tilde{y}_{j-1}\leq C$ by definition of $\tilde{\boldsymbol{y}}$. For $j=n+2,\ldots,S$, we have $\tilde{y}_{j}-\tilde{y}_{j-1}=0$. Finally, we have $\tilde{y}_{n+1}-\tilde{y}_{n}\leq D(T)-r_{n}<\eta\leq C$ (where we use Claim 3 to bound $\tilde{y}_{n}$). Thus, $\tilde{\boldsymbol{y}}$ satisfies constraint (i).

For $j=1,\ldots,n$, if $r_{j}>r_{j-1}$, we have $\tilde{y}_{j-1}\leq r_{j-1}+\eta\leq r_{j}\leq\tilde{y}_{j}$ (the last inequality being Claim 3) and if $r_{j}=r_{j-1}$, necessarily $r_{j}=r_{j-1}=0$ and $\tilde{y}_{j-1}=\tilde{y}_{j}=\eta$. Thus, $\tilde{\boldsymbol{y}}$ satisfies constraint (ii).

Consider $j\in\{2,\ldots,n\}$. We have

[TABLE]

The first inequality is obtained with the help of Claim 3. For the second one, we use Claim 4 and also that $z_{j}=r_{j}-r_{j-1}$ and $z_{j-1}=r_{j-1}-r_{j-2}$ which hold because \tilde{p}=\big{(}(z_{0},r_{0}),(z_{1},r_{1}),\ldots,(z_{n},r_{n})\big{)} is a path in $\mathcal{G}$. For the last inequality, we use $\mathrm{bar}\tau(r_{j})-\mathrm{bar}\tau(r_{j-1})+\nu(z_{j}-z_{j-1})+\frac{1}{2}\gamma\eta\geq 0$, which holds again because $\tilde{p}$ is a path, and the definition of $\gamma$. For $j\geq n+1$, we have $\tilde{d}_{j}\geq\tilde{d}_{j-1}$ by definition. Constraint (iii) is thus satisfied for all $j$.

If $n<S$, then $\tilde{y}_{S}=D(T)$ by definition. From now on, we suppose thus that $n=S$. We also suppose that $S\geq 2$. The case $S=1$ being easy to check (and anyway, for a complexity point of view, this case does not matter). If $\tilde{y}_{S-1}=r_{S-1}+\eta$, then $\tilde{y}_{S-1}+C=r_{S-1}+\eta+C\geq r_{S}+\eta>D(T)$ (here we use that $z_{S}\leq C$ and that $r_{S}=R\eta$) and thus $\tilde{y}_{S}=D(T)$. If $\tilde{y}_{S-1}=D(T)$, then $\tilde{y}_{S}=D(T)$ since $\tilde{y}_{S-1}\leq\tilde{y}_{S}\leq D(T)$. Hence, in all these cases, $\tilde{\boldsymbol{y}}$ satisfies constraint (iv). The only remaining case is when $\tilde{y}_{S-1}=\tilde{y}_{S-2}+C$. If $j$ is an index in $\left\{1,\ldots,S-2\right\}$ such that $\tilde{y}_{j}=r_{j}+\eta$, then we have $r_{j+1}+\eta\leq r_{j}+C+\eta=\tilde{y}_{j}+C$ and $r_{j+1}+\eta\leq D(T)$, and thus $\tilde{y}_{j+1}=r_{j+1}+\eta$. It implies that as soon as some $j_{0}\in\left\{1,\ldots,S-1\right\}$ is such that $\tilde{y}_{j_{0}}=r_{j_{0}}+\eta$, we have $\tilde{y}_{S-1}=r_{S-1}+\eta$, which is a case we have already dealt with. Since $r_{j}+\eta\leq r_{S}\leq D(T)$ for $j\in\{1,\ldots,S-1\}$, we are left with the case where $\tilde{y}_{j}=\tilde{y}_{j-1}+C$ for every $j\in\left\{1,\ldots,S-1\right\}$. In this situation, we have $\tilde{y}_{S-1}=(S-1)C$ and hence $\tilde{y}_{S-1}+C=CS\geq D(T)$. Since $r_{S}+\eta>D(T)$, we get that $\tilde{y}_{S}=D(T)$, and $\tilde{\boldsymbol{y}}$ satisfies constraint (iv) in every case.

For $j=1,\ldots,n$, we have $\tilde{d}_{j}\geq\mathrm{bar}\tau(\tilde{y}_{j})+\nu(\tilde{y}_{j}-\tilde{y}_{j-1})$ by definition, and for $j\geq n+1$, we have $\tilde{d}_{j}\geq T+\nu(\tilde{y}_{n+1}-\tilde{y}_{n})\geq\mathrm{bar}\tau(\tilde{y}_{j})+\nu(\tilde{y}_{j}-\tilde{y}_{j-1})$. Thus $\tilde{\boldsymbol{d}}$ satisfies constraint (v) and $(\tilde{\boldsymbol{d}},\tilde{\boldsymbol{y}})$ is feasible for P${}^{\operatorname{ave}}_{\mbox{{\tiny no return}}}$. ∎

Proof of Lemma 8.

Our goal is to bound from above the following quantity

[TABLE]

We proceed by splitting the expression into two parts: the sum from $j=1$ to $j=n$, and the sum from $j=n+1$ to $j=S$.

Using Claims 3 and 4, we have $\mathrm{bar}\tau(\tilde{y}_{j})+\nu(\tilde{y}_{j}-\tilde{y}_{j-1})\leq q_{j}+\eta/\alpha+\nu\eta$, where $q_{j}=\mathrm{bar}\tau(r_{j})+\nu(r_{j}-r_{j-1})$. Thus we have for all $j\leq n$,

[TABLE]

since $f^{\operatorname{ave}}(\makebox[4.30554pt]{\cdot})$ is nonincreasing in the second term and nondecreasing in the first and third terms and where we extend the definition of $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ by letting $\mathrm{bar}\tau(y)=T$ for all $y>D(T)$.

For the second part, we proceed as follows. Since $r_{n}+\eta=(R+1)\eta>D(T)$, Claim 3 immediately implies $D(T)-\tilde{y}_{n}\leq\eta$. With Claim 4, we get thus $T\leq\mathrm{bar}\tau(\tilde{y}_{n})+\eta/\alpha$, where we used T=\mathrm{bar}\tau\big{(}D(T)-\tilde{y}_{n}+\tilde{y}_{n}\big{)}. This provides

[TABLE]

Using again the fact that $f^{\operatorname{ave}}(\makebox[4.30554pt]{\cdot})$ is nonincreasing in the second term and nondecreasing in the first and third terms and with the help of Claim 3, we get

[TABLE]

since the terms indexed by $j=n+2,\ldots,S$ are all zero and since $D(T)<r_{n}+\eta$.

We aim at comparing the upper bounds in Equations (5) and (6) with

[TABLE]

We first compare the $j$th term of the bound in (5) with the $j$th term of the sum in (7).

[TABLE]

with

[TABLE]

Since $\mathrm{bar}\tau(\makebox[4.30554pt]{\cdot})$ in nondecreasing, we get

[TABLE]

Using $j\gamma\leq n\gamma$ and $\gamma=2(1/\alpha+2\nu)$, we obtain

[TABLE]

We now bound the term in Equation (6). Let $I=f^{\operatorname{ave}}(q_{n}+\eta/\alpha+\nu\eta+n\gamma\eta,r_{n},r_{n}+\eta)$. We have

[TABLE]

We have thus

[TABLE]

Using $r_{n}\leq D(T)$ and $\mathrm{bar}\tau(r_{n})\leq T$ leads to

[TABLE]

∎