In this Section, the dominated convergence theorem is used implicitly to
justify any differentiation under the integral sign. When dealing with
polyhedral functions such as
$\theta\longmapsto\max_{v\in\mathcal{V}}\theta^{\top}v$, we often write
$\nabla_{\theta}$ for simplicity even though they are only subdifferentiable,
because the set of non-differentiability has measure zero. We denote the
standard Gaussian density by $\nu$.

###  A.1 Additive perturbation

These proofs are already given by [6], but we include them for comparison
purposes.

####  A.1.1 Derivatives

Proof of Proposition 3.1.

###### Proof.

The following change of variable is a diffeomorphism:

$u=\theta+\varepsilon z\quad\iff\quad\frac{u-\theta}{\varepsilon}=z.$

We apply it to the definition of $\widehat{p}_{\varepsilon}^{+}(v|\theta)$.

$\displaystyle\widehat{p}_{\varepsilon}^{+}(v|\theta)$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathbbm{1}{\\{f(\theta+\varepsilon
z)=v\\}}\leavevmode\nobreak\ \nu(z)\leavevmode\nobreak\ \mathrm{d}z$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathbbm{1}{\\{f(u)=v\\}}\leavevmode\nobreak\
\nu\left(\frac{u-\theta}{\varepsilon}\right)\frac{\mathrm{d}u}{\varepsilon^{d}}.$

We now differentiate with respect to $\theta$ before applying the reverse
change of variable:

$\displaystyle\nabla_{\theta}\widehat{p}_{\varepsilon}^{+}(v|\theta)$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathbbm{1}{\\{f(u)=v\\}}\leavevmode\nobreak\
\left(\frac{-1}{\varepsilon}\nabla\nu\left(\frac{u-\theta}{\varepsilon}\right)\right)\frac{\mathrm{d}u}{\varepsilon^{d}}$
$\displaystyle=\frac{-1}{\varepsilon}\int_{\mathbb{R}^{d}}\mathbbm{1}{\\{f(\theta+\varepsilon
z)=v\\}}\leavevmode\nobreak\ \nabla\nu(z)\leavevmode\nobreak\ \mathrm{d}z$
$\displaystyle=\frac{1}{\varepsilon}\int_{\mathbb{R}^{d}}\mathbbm{1}{\\{f(\theta+\varepsilon
z)=v\\}}\leavevmode\nobreak\ z\nu(z)\leavevmode\nobreak\ \mathrm{d}z.$

The last equality holds because the standard Gaussian density satisfies
$\nabla\nu(z)=-z\nu(z)$. From there, we deduce the Jacobian of
$\widehat{f}_{\varepsilon}^{+}(\theta)$:

$\displaystyle J_{\theta}\widehat{f}_{\varepsilon}^{+}(\theta)$
$\displaystyle=\sum_{v\in\mathcal{V}}v\nabla_{\theta}\widehat{p}_{\varepsilon}^{+}(\theta,y)^{\top}$
$\displaystyle=\frac{1}{\varepsilon}\int_{\mathbb{R}^{d}}\underbrace{\Big{(}\sum_{v\in\mathcal{V}}v\leavevmode\nobreak\
\mathbbm{1}{\\{f(\theta+\varepsilon z)=v\\}}\Big{)}}_{f(\theta+\varepsilon
z)}\leavevmode\nobreak\ z^{\top}\nu(z)\leavevmode\nobreak\ \mathrm{d}z$

We arrive at the following simple expression, which was already given by [6]:

$J_{\theta}\widehat{f}_{\varepsilon}^{+}(\theta)=\frac{1}{\varepsilon}\mathbb{E}\left[f(\theta+\varepsilon
Z)Z^{\top}\right]$

∎

####  A.1.2 Regularization

Proof of Proposition 3.2.

###### Proof.

Because $\Omega_{\varepsilon}^{+}=(F_{\varepsilon}^{+})^{*}$ is a Fenchel
conjugate, it is automatically convex. Furthermore,

$\displaystyle\Omega_{\varepsilon}^{+}(\mu)$
$\displaystyle=\sup_{\theta\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-
F_{\varepsilon}^{+}(\theta)\right\\}$
$\displaystyle=\sup_{\theta\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-\mathbb{E}\left[\max_{v\in\mathcal{V}}(\theta+\varepsilon
Z)^{\top}v\right]\right\\}.$

We consider $\mu\notin\operatorname{conv}(\mathcal{V})$. By convex separation,
there exists $\tilde{\theta}\in\mathbb{R}^{d}$ and $\alpha>0$ such that
$\tilde{\theta}^{\top}\mu\geq\alpha+\tilde{\theta}^{\top}v$ for all
$v\in\mathcal{V}$. This implies that, for all $t>0$,

$\displaystyle\Omega_{\varepsilon}^{+}(\mu)$ $\displaystyle\geq
t\tilde{\theta}^{\top}\mu-\mathbb{E}\left[\max_{v\in\mathcal{V}}\left(t\tilde{\theta}+\varepsilon
Z\right)^{\top}v\right]$ $\displaystyle\geq
t\tilde{\theta}^{\top}\mu-t\mathbb{E}\left[\max_{v\in\mathcal{V}}\tilde{\theta}^{\top}v\right]-\varepsilon\mathbb{E}\left[\max_{v\in\mathcal{V}}Z^{\top}v\right]$
$\displaystyle\geq
t\alpha-\varepsilon\mathbb{E}\left[\max_{v\in\mathcal{V}}Z^{\top}v\right]\xrightarrow[t\to+\infty]{}+\infty$

We have shown that $\mu\notin\operatorname{dom}(\Omega_{\varepsilon}^{+})$,
and therefore
$\operatorname{dom}(\Omega_{\varepsilon}^{+})\subset\operatorname{conv}(\mathcal{V})$.
We define

$f_{Z}(\theta,v)=(\theta+\varepsilon Z)^{\top}v\qquad\text{so that}\qquad
F_{\varepsilon}^{+}(\theta)=\mathbb{E}\left[\max_{v\in\mathcal{V}}f_{Z}(\theta,v)\right].$

Danskin’s theorem helps us compute the gradient of $F_{\varepsilon}^{+}$:

$\displaystyle\nabla_{\theta}F_{\varepsilon}^{+}(\theta)$
$\displaystyle=\mathbb{E}\left[\nabla_{\theta}\left(\max_{v\in\mathcal{V}}f_{Z}(\theta,v)\right)\right]$
$\displaystyle=\mathbb{E}\left[\nabla_{1}f_{Z}\left(\theta,\operatorname*{argmax}_{v\in\mathcal{V}}f_{Z}(\theta,v)\right)\right]$
$\displaystyle=\mathbb{E}\left[\operatorname*{argmax}_{v\in\mathcal{V}}f_{Z}(\theta,v)\right]$
$\displaystyle=\mathbb{E}\left[\operatorname*{argmax}_{v\in\mathcal{V}}(\theta+\varepsilon
Z)^{\top}v\right]=\widehat{f}_{\varepsilon}^{+}(\theta).$

As shown by [6, Proposition 2.2], the function $\Omega_{\varepsilon}^{+}$ is a
Legendre type function, which means that

$\nabla_{\theta}F_{\varepsilon}^{+}=\nabla_{\theta}(\Omega_{\varepsilon}^{+})^{*}=(\nabla_{\theta}\Omega_{\varepsilon}^{+})^{-1}.$

From this, we deduce

$\displaystyle\nabla_{\theta}F_{\varepsilon}^{+}(\theta)$
$\displaystyle=\operatorname*{argmax}_{\mu\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-\Omega_{\varepsilon}^{+}(\mu)\right\\}$
$\displaystyle=\operatorname*{argmax}_{\mu\in\operatorname{dom}(\Omega_{\varepsilon}^{+})}\left\\{\theta^{\top}\mu-\Omega_{\varepsilon}^{+}(\mu)\right\\}=\widehat{f}_{\Omega_{\varepsilon}^{+}}(\theta).$

Hence, we can conclude:

$\widehat{f}_{\varepsilon}^{+}(\theta)=\nabla_{\theta}F_{\varepsilon}^{+}(\theta)=\operatorname*{argmax}_{\mu\in\operatorname{dom}(\Omega_{\varepsilon}^{+})}\\{\theta^{\top}\mu-\Omega_{\varepsilon}^{+}(\mu)\\}=\widehat{f}_{\Omega_{\varepsilon}^{+}}(\theta).$

To recover the formula given in Proposition 3.2, we simply remember that
$\operatorname{dom}(\Omega_{\varepsilon}^{+})\subseteq\operatorname{conv}(\mathcal{V})$.
∎

###  A.2 Multiplicative perturbation

####  A.2.1 Derivatives

Proof of Proposition 3.3.

###### Proof.

Suppose $\theta\in\mathbb{R}^{d}$ only has positive components. Then the
following change of variable is a diffeomorphism:

$u=\theta\odot e^{\varepsilon
z-\varepsilon^{2}\bm{1}/2}\quad\iff\quad\frac{\log(u)-\log(\theta)}{\varepsilon}+\frac{\varepsilon\bm{1}}{2}=z$

We apply it to the definition of
$\widehat{p}_{\varepsilon}^{\odot}(v|\theta)$.

$\displaystyle\widehat{p}_{\varepsilon}^{\odot}(\theta,y)$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathbbm{1}{\left\\{f\left(\theta\odot
e^{\varepsilon
z-\varepsilon^{2}\bm{1}/2}\right)=v\right\\}}\leavevmode\nobreak\
\nu(z)\leavevmode\nobreak\ \mathrm{d}z$
$\displaystyle=\int_{(0,+\infty)^{d}}\mathbbm{1}{\\{f(u)=v\\}}\leavevmode\nobreak\
\nu\left(\frac{\log(u)-\log(\theta)}{\varepsilon}+\frac{\varepsilon\bm{1}}{2}\right)\frac{\mathrm{d}u}{\varepsilon^{d}\prod_{i}u_{i}}.$

We now differentiate with respect to $\theta$ before applying the reverse
change of variable:

$\displaystyle\nabla_{\theta}\widehat{p}_{\varepsilon}^{\odot}(\theta,y)$
$\displaystyle=\int_{(0,+\infty)^{d}}\mathbbm{1}{\\{f(u)=v\\}}\leavevmode\nobreak\
\left(\frac{-1}{\varepsilon\theta}\odot\nabla\nu\left(\frac{\log(u)-\log(\theta)}{\varepsilon}+\frac{\varepsilon\bm{1}}{2}\right)\right)\frac{\mathrm{d}u}{\varepsilon^{d}\prod_{i}u_{i}}$
$\displaystyle=\frac{-1}{\varepsilon\theta}\odot\int_{\mathbb{R}^{d}}\mathbbm{1}{\left\\{f\left(\theta\odot
e^{\varepsilon
z-\varepsilon^{2}\bm{1}/2}\right)=v\right\\}}\leavevmode\nobreak\
\nabla\nu(z)\leavevmode\nobreak\ \mathrm{d}z$
$\displaystyle=\frac{1}{\varepsilon\theta}\odot\int_{\mathbb{R}^{d}}\mathbbm{1}{\left\\{f\left(\theta\odot
e^{\varepsilon
z-\varepsilon^{2}\bm{1}/2}\right)=v\right\\}}\leavevmode\nobreak\
z\nu(z)\leavevmode\nobreak\ \mathrm{d}z.$

From there, we deduce the Jacobian of
$\widehat{f}_{\varepsilon}^{\odot}(\theta)$:

$\displaystyle J_{\theta}\widehat{f}_{\varepsilon}^{\odot}(\theta)$
$\displaystyle=\sum_{v\in\mathcal{V}}v\nabla_{\theta}\widehat{p}_{\varepsilon}^{\odot}(\theta,y)^{\top}$
$\displaystyle=\frac{1}{\varepsilon\theta}\odot\int_{\mathbb{R}^{d}}\underbrace{\Big{(}\sum_{v\in\mathcal{V}}v\leavevmode\nobreak\
\mathbbm{1}{\left\\{f\left(\theta\odot e^{\varepsilon
z-\varepsilon^{2}\bm{1}/2}\right)=v\right\\}}\Big{)}}_{f\left(\theta\odot
e^{\varepsilon z-\varepsilon^{2}\bm{1}/2}\right)}\leavevmode\nobreak\
z^{\top}\nu(z)\leavevmode\nobreak\ \mathrm{d}z$

We arrive at a simple variant of the previous expression:

$J_{\theta}\widehat{f}_{\varepsilon}^{\odot}(\theta)=\frac{1}{\varepsilon\theta}\odot\mathbb{E}\left[f\left(\theta\odot
e^{\varepsilon z-\varepsilon^{2}\bm{1}/2}\right)Z^{\top}\right]$

∎

####  A.2.2 Regularization

Proof of Proposition 3.4.

###### Proof.

Because $\Omega_{\varepsilon}^{\odot}=(F_{\varepsilon}^{\odot})^{*}$ is a
Fenchel conjugate, it is automatically convex. Furthermore,

$\displaystyle\Omega_{\varepsilon}^{\odot}(\mu)$
$\displaystyle=\sup_{\theta\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-
F_{\varepsilon}^{\odot}(\theta)\right\\}$
$\displaystyle=\sup_{\theta\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-\mathbb{E}\left[\max_{v\in\mathcal{V}}\left(\theta\odot
e^{\varepsilon Z-\varepsilon^{2}\bm{1}/2}\right)^{\top}v\right]\right\\}$
$\displaystyle=\sup_{\theta\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-\mathbb{E}\left[\max_{v\in\mathcal{V}}\theta^{\top}\left(v\odot
e^{\varepsilon Z-\varepsilon^{2}\bm{1}/2}\right)\right]\right\\}.$

This last expression shows why we do not have
$\operatorname{dom}(\Omega_{\varepsilon}^{\odot})\subset\operatorname{conv}(\mathcal{V})$
(unlike in the additive case). Indeed, even when
$\mu\notin\operatorname{conv}(\mathcal{V})$, the multiplicative scaling of $v$
might allow it to compensate the inner product $\theta^{\top}\mu$ and stop
$\Omega_{\varepsilon}^{\odot}(\mu)$ from going to $+\infty$. We define

$f_{Z}(\theta,v)=\left(\theta\odot e^{\varepsilon
Z-\varepsilon^{2}\bm{1}/2}\right)^{\top}v\qquad\text{so that}\qquad
F_{\varepsilon}^{\odot}(\theta)=\mathbb{E}\left[\max_{v\in\mathcal{V}}f_{Z}(\theta,v)\right].$

Danskin’s theorem helps us compute the gradient of $F_{\varepsilon}^{\odot}$:

$\displaystyle\nabla_{\theta}F_{\varepsilon}^{\odot}(\theta)$
$\displaystyle=\mathbb{E}\left[\nabla_{\theta}\left(\max_{v\in\mathcal{V}}f_{Z}(\theta,v)\right)\right]$
$\displaystyle=\mathbb{E}\left[\nabla_{1}f_{Z}\left(\theta,\operatorname*{argmax}_{v\in\mathcal{V}}f_{Z}(\theta,v)\right)\right]$
$\displaystyle=\mathbb{E}\left[e^{\varepsilon
Z-\varepsilon^{2}\bm{1}/2}\odot\operatorname*{argmax}_{v\in\mathcal{V}}f_{Z}(\theta,v)\right]$
$\displaystyle=\mathbb{E}\left[e^{\varepsilon
Z-\varepsilon^{2}\bm{1}/2}\odot\operatorname*{argmax}_{v\in\mathcal{V}}\left(\theta\odot
e^{\varepsilon Z-\varepsilon^{2}\bm{1}/2}\right)^{\top}v\right]$
$\displaystyle=\widehat{f}_{\varepsilon}^{\odot\mathrm{scaled}}(\theta)\neq\widehat{f}_{\varepsilon}^{\odot}(\theta)$

We could prove in a way similar to [6, Proposition 2.2] that
$\Omega_{\varepsilon}^{\odot}$ is a Legendre type function, which means that

$\nabla_{\theta}F_{\varepsilon}^{\odot}=\nabla_{\theta}(\Omega_{\varepsilon}^{\odot})^{*}=(\nabla_{\theta}\Omega_{\varepsilon}^{\odot})^{-1}.$

From this, we deduce

$\displaystyle\nabla_{\theta}F_{\varepsilon}^{\odot}(\theta)$
$\displaystyle=\operatorname*{argmax}_{\mu\in\mathbb{R}^{d}}\left\\{\theta^{\top}\mu-\Omega_{\varepsilon}^{\odot}(\mu)\right\\}$
$\displaystyle=\operatorname*{argmax}_{\mu\in\operatorname{dom}(\Omega_{\varepsilon}^{\odot})}\left\\{\theta^{\top}\mu-\Omega_{\varepsilon}^{\odot}(\mu)\right\\}=\widehat{f}_{\Omega_{\varepsilon}^{+}}(\theta).$

This time we cannot replace $\operatorname{dom}(\Omega_{\varepsilon}^{\odot})$
by $\operatorname{conv}(\mathcal{V})$, but we still obtain a similar
conclusion:

$\widehat{f}_{\varepsilon}^{\odot\mathrm{scaled}}(\theta)=\nabla_{\theta}F_{\varepsilon}^{\odot}(\theta)=\operatorname*{argmax}_{\mu\in\operatorname{dom}(\Omega_{\varepsilon}^{\odot})}\\{\theta^{\top}\mu-\Omega_{\varepsilon}^{\odot}(\mu)\\}=\widehat{f}_{\Omega_{\varepsilon}^{\odot}}(\theta).$

∎

###  A.3 Inexact oracles

Proof of Proposition 3.5.

###### Proof.

We start with additive perturbation. By Proposition 3.1, we have:

$\displaystyle
J_{\theta}\widehat{g}_{\varepsilon}^{+}(\theta)-J_{\theta}\widehat{f}_{\varepsilon}^{+}(\theta)$
$\displaystyle=\frac{1}{\varepsilon}\mathbb{E}\left[g(\theta+\varepsilon
Z)Z^{\top}\right]-\frac{1}{\varepsilon}\mathbb{E}\left[f(\theta+\varepsilon
Z)Z^{\top}\right]$
$\displaystyle=\frac{1}{\varepsilon}\mathbb{E}\left[\left(g(\theta+\varepsilon
Z)-f(\theta+\varepsilon Z)\right)Z^{\top}\right].$

We bound the spectral norm of the error using Jensen’s inequality:

$\displaystyle\left\lVert
J_{\theta}\widehat{g}_{\varepsilon}^{+}(\theta)-J_{\theta}\widehat{f}_{\varepsilon}^{+}(\theta)\right\lVert^{2}$
$\displaystyle\leq\frac{1}{\varepsilon^{2}}\mathbb{E}\left\lVert\left(g(\theta+\varepsilon
Z)-f(\theta+\varepsilon Z)\right)Z^{\top}\right\rVert^{2}$
$\displaystyle=\frac{1}{\varepsilon^{2}}\mathbb{E}\left[\lVert
g(\theta+\varepsilon Z)-f(\theta+\varepsilon Z)\rVert^{2}\lVert
Z\rVert^{2}\right]$ $\displaystyle\leq\frac{1}{\varepsilon^{2}}\lVert
g-f\rVert_{\infty}^{2}\mathbb{E}\left[\lVert Z\rVert^{2}\right]$
$\displaystyle=\frac{d}{\varepsilon^{2}}\lVert g-f\rVert_{\infty}^{2}.$

We now move on to multiplicative perturbation. For multiplicative
perturbation, following the same proof starting from Proposition 3.3 yields a
similar inequality:

$\displaystyle\left\lVert
J_{\theta}\widehat{g}_{\varepsilon}^{\odot}(\theta)-J_{\theta}\widehat{f}_{\varepsilon}^{\odot}(\theta)\right\lVert^{2}$
$\displaystyle\leq\frac{d}{\varepsilon^{2}\min_{i}\lvert\theta_{i}\rvert^{2}}\lVert
g-f\rVert_{\infty}^{2}.$

∎

##  Appendix B More details on applications

###  B.1 Stochastic vehicle scheduling problem

####  B.1.1 Linear program formulation of (VSP):

The deterministic Vehicle Scheduling Problem (VSP), can be formulated as the
following linear program:

$\displaystyle\min$ $\displaystyle\,\sum_{a\in A}\theta_{a}y_{a}$ (29)
$\displaystyle s.t.$
$\displaystyle\sum_{a\in\delta^{-}(v)}y_{a}=\sum_{a\in\delta^{+}(v)}y_{a},$
$\displaystyle\forall v\in V\backslash\\{o,d\\}$
$\displaystyle\sum_{a\in\delta^{-}(v)}y_{a}=1,$ $\displaystyle\forall v\in
V\backslash\\{o,d\\}$ $\displaystyle y_{a}\in\\{0,1\\},$ $\displaystyle\forall
a\in A$

Since the constraint matrix is totally unimodular, the constraint polytope
gives a perfect formulation, therefore binary constraints can be relaxed, and
this problem can be solved using either a flow algorithm or a linear
programming solver.

####  B.1.2 Two MILP formulations for the StoVSP

##### MILP formulation

We denote respectively $c_{\text{vehicle}}$ and $c_{\text{delay}}$ the cost of
one vehicle, and the cost of one minute of delay. The StoVSP problem can be
formulated as the following MIP with quadratic constraints:

$\displaystyle\min_{d,y}$
$\displaystyle\,c_{\text{delay}}\dfrac{1}{|S|}\sum\limits_{s\in
S}\sum\limits_{v\in
V\backslash\\{o,d\\}}d_{v}^{s}+c_{\text{vehicle}}\sum\limits_{a\in\delta^{+}(o)}y_{a}$
(30a) s.t.
$\displaystyle\sum_{a\in\delta^{-}(v)}y_{a}=\sum_{a\in\delta^{+}(v)}y_{a}$
$\displaystyle\forall v\in V\backslash\\{o,d\\}$ (30b)
$\displaystyle\sum_{a\in\delta^{-}(v)}y_{a}=1$ $\displaystyle\forall v\in
V\backslash\\{o,d\\}$ (30c) $\displaystyle
d_{v}^{s}\geq\gamma_{v}^{s}+\sum_{\begin{subarray}{c}a\in\delta^{-}(v)\\\
a=(u,v)\end{subarray}}(d_{u}^{s}-s_{u,v}^{s})y_{a}$ $\displaystyle\forall v\in
V\backslash\\{o,d\\},\forall s\in S$ (30d) $\displaystyle
d_{v}^{s}\geq\gamma_{v}^{s}$ $\displaystyle\forall v\in
V\backslash\\{o,d\\},\forall s\in S$ (30e) $\displaystyle y_{a}\in\\{0,1\\}$
$\displaystyle\forall a\in A$ (30f)

These quadratic delay constraints (30b) can be linearized using McCormick
inequalities, and we can then solve the resulting MILP with industrial
solvers. This approach does not scale well on instances with large number of
tasks and scenarios.

##### Column generation formulation

Another possible formulation for StoVSP is the column generation approach. Let
$\mathcal{P}$ the set of feasible vehicle routes. Let us define the cost of a
vehicle route $P\in\mathcal{P}$ by
$c_{P}^{s}=c_{\text{vehicle}}+c_{\text{delay}}\times\sum_{v\in P}d_{v}^{s}$.
We introduce decision variables $y_{P}\in\\{0,1\\}$ that equals $1$ if the
vehicle route $P$ is chosen. We can then formulate StoVSP as follows:

$\displaystyle\min$ $\displaystyle\frac{1}{|S|}\sum_{s\in
S}\sum_{P\in\mathcal{P}}c_{P}^{s}y_{P}$ (31) s.t. $\displaystyle\sum_{p\ni
v}y_{P}=1$ $\displaystyle\forall v\in V\backslash\\{o,d\\}$
$\displaystyle\quad(\lambda_{v}\in\mathbb{R})$ $\displaystyle
y_{P}\in\\{0,1\\}$ $\displaystyle\forall p\in\mathcal{P}$

Since there is an exponential number of variables, this formulation can be
solved using a column generation algorithm. We denote $\lambda_{v}$ the dual
variables of (31). The associated sub-problem is a constrained shortest path
problem :

$\min_{P\in\mathcal{P}}(c_{P}-\sum_{v\in P}\lambda_{v})$ (32)

The sub-problem (32) can be solved using a stochastic constrained shortest
path algorithm (see [41] for in-depth details about these algorithms, and
ConstrainedShortestPaths.jl202020https://github.com/BatyLeo/ConstrainedShortestPaths.jl
for a Julia implementation).

Solving (31) using column generation gives good bounds on instances with up to
50 tasks. A Branch-and-price could be used to try to find the optimal
solution.

####  B.1.3 Data generation

##### Instance generator

To generate our dataset, we use a similar approach as by [43]. We define a
squared map of 50 minutes width, divided in 25 squared ($10\times 10$)
districts. Each task $v$ has a (uniformly drawn) start point, a start time
$t_{v}^{b}$, end point, and end time $t_{v}^{e}$.

Figure 9: Example of a map with 10 tasks. Task 0 represents the depot, in the
middle of the map. Square with task indices are starting points, and tasks
without are end points.

We draw each scenario $s\in S$ separately and independently. For each district
$d$ and hour $h$ of the day, we draw
$\beta_{d,h}\sim\log\mathcal{N}(\mu_{d,h},\sigma_{d,h})$, with
$(\mu_{d,h},\sigma_{d,h})\sim\mathcal{U}([1,3])\times\mathcal{U}([0.4,0.6])$.
For each district $d$ and hour $h$ of the day, $\zeta_{d,h}^{dis}$ models the
congestion in the district at this time, and is computed as follows:

$\begin{cases}\forall d,\,\zeta^{dis}_{d,0}=\beta_{d,0}\\\ \forall
d,h,\,\zeta^{dis}_{d,h+1}=\frac{1}{2}\zeta^{dis}_{d,h}+\beta_{d,h}\end{cases}$
(33)

For each hour $h$ of the day, $\zeta^{inter}_{h}$ models the congestion on
roads between districts, and is computed similarly. With
$I\sim\log\mathcal{N}(\mu=0.02,\sigma=0.05)$:

$\begin{cases}\zeta^{inter}_{0}=I\\\
\zeta^{inter}_{h+1}=(\zeta^{inter}_{h}+0.1)I\end{cases}$ (34)

Let $v$ be a task corresponding to a trip between district $d_{1}$ and
$d_{2}$. We compute the perturbed start and end times are:

$\begin{cases}\xi_{v}^{b}=t_{v}^{b}+\beta_{v}\\\
\xi_{v}^{e}=\xi_{v}^{b}+t_{v}^{e}-t_{v}^{b}+\zeta^{dis}_{d_{1},h(\xi_{1})}+\zeta^{inter}_{h(\xi_{2})}+\zeta^{dis}_{d_{2},h(\xi_{3})}\end{cases}$
(35)

with $\xi_{1}=\xi_{v}^{b}$, $\xi_{2}=\xi_{1}+\zeta^{dis}_{d_{1},h(\xi_{1})}$,
and $\xi_{3}=\xi_{2}+t_{v}^{e}-t_{v}^{b}+\xi^{inter}_{h(\xi_{2})}$. And the
perturbed travel times along arcs $a=(u,v)$:

$\xi_{a}^{tr}=\xi_{v}^{b}+t_{a}^{tr}+\zeta^{dis}_{d_{1},h(\xi_{1})}+\zeta^{inter}_{h(\xi_{2})}+\zeta^{dis}_{d_{2},h(\xi_{3})}$
(36)

with $\xi_{1}=\xi_{u}^{e}$, $\xi_{2}=\xi_{1}+\zeta^{dis}_{d_{1},h(\xi_{1})}$,
and $\xi_{3}=\xi_{2}+t_{a}^{tr}+\xi^{inter}_{h(\xi_{2})}$.

##### Encoding features

For an instance, we compute a matrix of 20 features for every arc of the
corresponding graph. Let $a=(u,v)\in A$. The first feature is the length of
$a$ in minutes, i.e. the deterministic travel time between $u$ and $v$. The
second feature is equal to $c_{\text{vehicle}}$ if $a$ is connected to the
source (i;.e if $u=o$), else 0. Then, we add the 9 deciles of
$\xi_{v}^{b}-(\xi_{u}^{e}+\xi_{a}^{tr})$ which represents the slack between
$u$ and $v$. Finally, we add the values of the cumulative probability
distribution of the slack, evaluated in $-100$, $-50$, $-20$, $-10$, $0$,
$50$, $200$, and $500$.

##### Label solutions

We label each instance with a local search heuristic initialized with the
solution of the deterministic problem. Its implementation can be found on the
GitHub
repository212121https://github.com/BatyLeo/StochasticVehicleScheduling.jl.

####  B.1.4 Training plots

Figure 10: Learning by imitation metric evolution during learning Figure 11:
Learning by experience metric evolution during learning

###  B.2 Two-stage minimum weight spanning tree

This subsection contains additional material about the two-stage stochastic
spanning tree problem.

####  B.2.1 Lagrangian heuristic

We present the heuristic based on Lagrangian relaxation used to label our
datasets.

##### Lagrangian relaxation

An option to compute a lower bound on the optimal value of the _two-stage
minimum weight spanning tree_ is to apply a Lagrangian relaxation. For this,
we duplicate complicating variables $y$ for each scenario and introduce
associated $\theta_{es}$ dual variable for each new constraint:

$\begin{array}[]{rll}\min&\displaystyle\sum_{e\in
E}c_{e}y_{e}+\frac{1}{|S|}\sum_{e\in E}\sum_{s\in S}d_{es}z_{es}\\\
\mathrm{s.t.}&\displaystyle\sum_{e\in E}y_{e}+z_{es}=|V|-1&\forall s\in S\\\
&\displaystyle\sum_{e\in E(Y)}y_{e}+z_{es}\leq|Y|-1,&\forall
Y,\,\emptyset\subsetneq Y\subsetneq E,\,\forall s\in S\\\
&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}y_{es}=y_{e},}&{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\forall
e\in E,\,\forall s\in S\text{ ($\theta_{es}$ dual variables)}}\\\
&y,z\in\\{0,1\\}\end{array}$ (37)

Lagrangian dual problem can be formulated as follows:

$\max_{\theta}G(\theta)=\left\\{\begin{aligned} \min_{y}\quad&\sum_{e\in
E}(c_{e}+\dfrac{1}{|S|}\sum_{s\in S}\theta_{es})y_{e}\\\
\mathrm{s.t.}\quad&0\leq y\leq M\end{aligned}\right.+\frac{1}{|S|}\sum_{s\in
S}\left\\{\begin{aligned} \min_{y_{s},z_{s}}\quad&\sum_{e\in
E}d_{es}z_{es}-\theta_{es}y_{es}\\\
\mathrm{s.t.}\quad&y_{s}+z_{s}\in\text{spanning tree
polytope}\end{aligned}\right.,$ (38)

where $M$ is a large constant. In theory $M$ is not needed, we can take
$M=+\infty$, but taking a finite $M$ leads to more informative gradients. $G$
can be maximized by gradient ascent:

$(\nabla_{\theta}G(\theta))_{es}=\frac{1}{|S|}(y_{e}-y_{es})$ (39)

##### Lagrangian heuristic

Once we find a good $\theta$ with the Lagrangian relaxation, we can retrieve
the corresponding $y_{es}$ solution, and use it as input of the Lagrangian
heuristic 1. The general idea is to rebuild a forest from the edges that are
selected in the first stage for most scenarios.

Data: $y_{es}\in\\{0,1\\},\,\forall e\in E,\,\forall s\in S$

forall  _$e\in E$_ do

if  _$\frac{1}{|S|}\sum_{s\in S}y_{es} >0.5$_ then

$w_{e}\leftarrow 1$;

else

$w_{e}\leftarrow-1$;

end if

end forall

Find maximum weight forest for weights $w$;

Result: $y_{e}\in\\{0,1\\},\,\forall e\in E$

Algorithm 1 Lagrangian heuristic

In the numerical experiments, we use the Lagrangian heuristic to build the
training set. Practically, we stop the Lagrangian relaxation after 50.000
iterations or if the gap between the lower bound provided by the $G(\theta)$
and the upper bound obtained from the heuristic is non-grater than 0.1%.

####  B.2.2 Instance features

We detail here all the features used for training.

##### Basic feature set

Basic features only contain twelve features. Let $e$ be an edge. The first one
is the first stage cost $c_{e}$, and the other 11 are the quantiles of the
second stage costs $(d_{es})_{s\in S}$ for values
$\\{0,\,0.1,\,0.2,\,0.3,\,0.4,\,0.5,\,0.6,\,0.7,\,0.8,\,0.9,\,1\\}$. All
quantiles in the advanced feature set use the same values.

##### Advanced feature set

Advanced features contain 79 features in total. Let $e$ be an edge. The first
12 features are the same as the basic features. To these, we add quantiles of
best stage costs (i.e for each scenario $s$ the value $\min(c_{e},d_{es})$),
quantiles of neighbors first stage costs, and quantiles of neighbors second
stage costs. Then, we compute several single stage minimum spanning trees with
Kruskal’s algorithm: one with first stage costs, one for each scenario with
corresponding second stage costs, and one for each scenario with corresponding
best stage costs $\min(c_{e},d_{es})$. We use these spanning trees to compute
new useful features: a binary features that tells if $e$ is in the first stage
spanning tree, and quantile features on its presence in second stage trees.
Finally, we use best stage spanning trees to compute quantiles for the cost of
$e$ in the best stage costs spanning tree: one set of quantiles for first edge
costs, and another set for second stage costs.