# Enhanced transport of flexible fibers by pole vaulting in turbulent wall-
bounded flow

Jérémie Bec  Université Côte d’Azur, Inria, CNRS, Calisto team, Sophia
Antipolis, France  Mines Paris, PSL University, CNRS, Cemef, Sophia Antipolis,
France  Christophe Brouzet  Université Côte d’Azur, CNRS, Institut de Physique
de Nice, France  Christophe Henry  Université Côte d’Azur, Inria, CNRS,
Calisto team, Sophia Antipolis, France

###### Abstract

Long, flexible fibers transported by a turbulent channel flow sample non-
linear variations of the fluid velocity along their length. As the fibers
tumble and collide with the boundaries, they bounce off with an impulse that
propels them toward the center of the flow, similar to pole vaulting. As a
result, the fibers migrate away from the walls, leading to depleted regions
near the boundaries and more concentrated regions in the bulk. These higher
concentrations in the center of the channel result in a greater net flux of
fibers than what was initially imposed by the fluid. This effect becomes more
pronounced as fiber length increases, especially when it approaches the
channel height.

The concentration and orientation of fibers in wall-bounded turbulent flows
have been extensively studied in the past decades [1, 2]. Earliest work has
focused on numerical simulations of rigid fibers, either shorter [2, 3, 4] or
longer [5] than the smallest length scale of the flow. These simulations have
shown that the fibers tend to adopt a random orientation in the center of the
channel, while aligning with the streamwise direction close to the walls.
Near-wall coherent structures were found to play an important role, as fibers
could accumulate in vortices, depending on their length and inertia, leading
to large inhomogeneities in concentration and velocity distribution. Recent
experiments [6, 7, 8, 9, 10, 11, 12, 13] confirmed these findings.
Accumulation in near-wall structures can result from various effects, such as
interactions with high-speed streaks [11], turbophoresis (that causes inertial
particles to migrate toward the wall), as well as interactions with
boundaries, especially for rod-like objects [14, 8, 5, 12]. Still, these
studies focused on rigid fibers much shorter than the channel height. It is
expected that longer and more flexible fibers will be deformed by the flow,
leading to new behaviors not yet fully explored.

Recent simulations have demonstrated that flexible fibers with inertia also
tend to accumulate in the near-wall region [15, 16]. The fibers considered
there are longer than the smallest scale of the flow, but still much shorter
than the channel height. To date, there are no studies on the behavior of
flexible fibers with sizes comparable to the channel height and, to the best
of our knowledge, they have only been investigated in homogeneous isotropic
turbulence. In this context, laboratory experiments have characterized the
rigid-flexible transition [17] and curvature statistics [18] of long flexible
fibers. Long elastic chains have been simulated [19, 20], suggesting that they
tend to align with vortex filaments and remain trapped for extended periods.
Other numerical studies on long flexible fibers in homogeneous isotropic
turbulence have demonstrated their potential for measuring two-point
statistics of turbulence [21], or exploring modulation mechanisms of
turbulence [22]. The aim of this study is to investigate the near-wall
dynamics of elongated flexible fibers induced by collisions with boundaries,
filling the gap in our understanding of their behavior in wall-bounded
turbulent flows.

In this study, our focus is on long, inextensible, and inertialess fibers
immersed in a channel flow with moderate Reynolds number. We consider a one-
way coupling between the fibers and the flow [23], which means that the fibers
do not produce any feedback on the flow nor interact with one another. The
over-damped slender-body theory is used to describe the fibers, whereby each
fiber, with length $\ell$ and thickness $a$, is parametrized by its arc-length
$s\in[-\tfrac{\ell}{2},\tfrac{\ell}{2}]$, and its position
${\boldsymbol{X}}(s,t)$ is given by:

$\displaystyle\partial_{t}{\boldsymbol{X}}$ $\displaystyle=$
$\displaystyle{\boldsymbol{u}}({\boldsymbol{X}},t)+\frac{c}{8\pi\rho_{\mathrm{f}}\nu}\mathbb{D}\left[\partial_{s}(T\partial_{s}{\boldsymbol{X}})-E\partial_{s}^{4}{\boldsymbol{X}}\right],$
with
$\displaystyle\mathbb{D}=\mathbb{1}+\partial_{s}{\boldsymbol{X}}\partial_{s}{\boldsymbol{X}}^{\mathsf{T}}\quad\mbox{and}\quad|\partial_{s}{\boldsymbol{X}}|=1.$
(1)

The tension $T(s,t)$ is here a Lagrange multiplier associated to the
inextensibility constraint, $E=YI$ is the fiber’s bending modulus (with $Y$
the Young modulus and $I=\pi a^{4}/64$ the fiber’s moment of inertia),
$\rho_{\rm f}$ the fluid density, $\nu$ its kinematic viscosity and
$c=-[1+2\,\text{log}(a/\ell)]$ a parameter related to the shape [24]. Equation
(1) involves the fluid velocity ${\boldsymbol{u}}({\boldsymbol{X}},t)$ along
the fiber. It is obtained in our numerical simulations from the open-source
spectral code Channelflow 2.0 [25], used to solve the incompressible
Navier–Stokes equation in a rectangular channel of size $L_{x}=4\pi h$,
$L_{y}=2h$, $L_{z}=4\pi h/3$ (see Fig. 1), periodic in the $x$ and $z$
directions and with no-slip boundary conditions at $y=\pm h$. The fluid flow
is integrated with a resolution of $128\times 129\times 128$ collocation
points. It is maintained in motion by imposing a constant bulk velocity, such
that the friction Reynolds number is equal to
$\textrm{Re}_{\tau}=u_{\tau}\,h/\nu\approx 180$ (with $u_{\tau}$ the friction
velocity). From now on, we will be using classical wall units (denoted by a
superscript $+$), where time and length scales are non-dimensionalized by
$\tau_{\nu}=\nu/u_{\tau}^{2}$ and $\delta_{\nu}=\nu/u_{\tau}$.

Figure 1: Instantaneous distributions of fibers. (a) $\ell^{+}=30$ and (b)
$\ell^{+}=270$. The colored panels show the modulus of fluid vorticity in the
planes $x=L_{x}$, $y=-L_{y}/2=h$ and $z=0$. Note that near-wall structures or
streaks are visible close to the plane $y=-h$. The mean flow is from left to
right.

The dynamics of fibers in Eq. (1) is integrated numerically by a second-order
finite-difference scheme [26, 24]. To prevent them from crossing the walls at
$y=\pm h$, we adopt a soft-boundary approach where the fluid velocity outside
the domain is replaced by ${\boldsymbol{u}}=(0,v,0)$, with $v=-\gamma\,(y+h)$
for $y<-h$ and $v=-\gamma\,(y-h)$ for $y>h$. The parameter $\gamma$ is
adjusted to the time step $\Delta t$ as $\gamma=0.5/\Delta t$. In our
simulations, fiber lengths take values in the range
$\ell^{+}=\ell/\delta_{\nu}\in[15,270]$, i.e. $\ell/h\in[0.08,1.5]$. Each
fiber is discretized with a given number of grid points $N_{s}$ along its arc-
length, chosen to keep a constant step size $\delta
s^{+}=\ell^{+}/N_{s}\approx 0.9$. The number of simulated fibers $N_{\rm
fib}\in[100,1250]$ decreases with increasing length $\ell^{+}$ for the sake of
computational time.

In the model we are using, inextensible and inertialess fibers are
characterized by two length scales. The first scale is their extension length
$\ell$, which is expressed in wall units as $\ell^{+}=\ell/\delta_{\nu}$. This
value indicates how far into the turbulent boundary layer the fibers can
extend. The second scale is the elastic length $\ell_{\rm
E}=\left[c\,E\tau_{\nu}/(8\pi\,\rho_{\mathrm{f}}\,\nu)\right]^{1/4}$, which is
maintained constant to $\ell_{\rm E}^{+}\approx 7.1$ in our simulations. It is
determined by balancing in Eq. (1) the bending rate of the fibers, given by
$cE/(8\pi\rho_{\mathrm{f}}\nu\ell^{4})$, with the typical shear rate
$\tau_{\nu}^{-1}$ of the channel flow. Note that the definition of elastic
length introduced in [17] for homogeneous isotropic turbulence is recovered
when using the Kolmogorov time scale $\tau_{\eta}$ instead of $\tau_{\nu}$.
The dimensionless flexibility of the fibers, defined as
$\mathcal{F}=\ell/\ell_{\rm E}$, measures their ability to bend near the walls
of the channel, with smaller values of $\mathcal{F}$ indicating more rigid
fibers. In order to examine how flexible fibers behave in comparison to rigid
rods (which are formally defined as having $\ell^{+}=0$ and $\mathcal{F}=0$),
we also integrated the Jeffery equations for threadlike ellispoids in the same
flow along tracer trajectories. Figure 1 shows two snapshots of the fiber
distribution with different lengths. Upon qualitative inspection, it is
apparent that short, rigid fibers remain very straight, whereas long, flexible
fibers bend into complex shapes.

We first report results on the statistics of the fiber orientation and
curvature as a function of their distance to the wall. Figure 2(a) represents
the average orientation of the local tangent vector to the fiber conditioned
on the distance $y^{+}$ to the nearest boundary. Only the streamwise ($x$) and
wall-normal ($y$) components are represented, the spanwise component ($z$)
being deduced from the inextensibility condition
$|\partial_{s}{\boldsymbol{X}}|^{2}=1$. Our findings confirm the previously-
observed behavior of stiff rods ($\ell^{+}=0$, black lines) that tend to align
with the mean flow near the boundaries and have a random orientation in the
bulk [2]. Finite-size, flexible fibers show a qualitatively similar trend.
However, their preferential orientation with the mean flow persists at larger
distances, in particular for the longer fibers. Notably, longer fibers show a
non-zero alignment with the wall-normal direction in the viscous sublayer, as
observed for flexible fibers with inertia [15].

Figure 2: Near-wall alignment and tumbling. (a) Mean-square components of the
tangent $\partial_{s}{\boldsymbol{X}}$ to the fiber in the streamwise (solid
lines) and wall-normal (dashed lines) directions, as a function of the
distance $y^{+}$ to the wall and for various fibers lengths $\ell^{+}$. The
horizontal dashed line corresponds to an isotropic orientation
$\langle(\partial_{s}X_{i})^{2}\rangle=1/3$. (b) Root-mean square amplitude of
the tumbling rate $\partial_{t}\partial_{s}{\boldsymbol{X}}$, again
conditioned on $y^{+}$. The dashed line shows here the value
$\langle|\partial_{t}\partial_{s}{\boldsymbol{X}}|^{2}\rangle\approx
0.1/\tau^{2}_{\eta}$ expected for rods in homogeneous isotropic turbulence
[27, 2], where the Kolmogorov time $\tau_{\eta}=(\nu/\varepsilon)^{1/2}$ is
computed using the bulk-flow turbulent dissipation rate $\varepsilon$.

Regarding angular velocities, there are several ways to define a tumbling rate
for flexible particles. We measure here the amplitude
$|\partial_{t}\partial_{s}{\boldsymbol{X}}|$ of the time derivative of the
tangent vector along the fiber arc-length. This definition retrieves when
$\ell\to 0$ the well-known tumbling rate of rigid fibers
$|\mathrm{d}\boldsymbol{p}/\mathrm{d}t|$, where $\boldsymbol{p}$ is the rod’s
orientation [2, 28]. Figure 2(b) displays the root-mean-square value of the
fibers tumbling rate conditioned on the wall-normal distance $y^{+}$. In the
channel center, they are all comparable and approach the value expected for
rods in homogeneous isotropic turbulence [27, 2]. Interestingly, we observe
that the smallest fibers with $\ell^{+}=15$ tumble similarly to rigid rods,
except in close proximity to the wall, at distances $y^{+}\lesssim\ell^{+}/2$,
where tumbling is slightly enhanced compared to rods. For intermediate sizes,
tumbling is depleted far from the boundary but enhanced close to the wall.
This effect becomes significantly more pronounced for the longest fibers
($\ell^{+}=180,270$), where tumbling is drastically amplified in close
proximity to the wall. These strong angular velocities may result from wall
contacts [12] and, as we will see, give rise to important motions away from
the boundary.

Insights into the near-boundary dynamics can be gained by considering a simple
two-dimensional model where the fiber is confined in the $(x,y)$ plane and the
flow over the wall is a pure linear shear flow
${\boldsymbol{u}}=\sigma\,y\,\hat{\boldsymbol{x}}$. In this scenario the
elastic length directly reads $\ell_{\rm
E}=\left[c\,E/(8\pi\,\rho_{\mathrm{f}}\,\nu\sigma)\right]^{1/4}$. If a fiber
is initially straight but slightly slanted downward, it will undergo a
tumbling motion following a Jeffery orbit. If the fiber’s center of mass is
initially within a distance $\ell/2$ from the wall, the fiber will touch the
boundary during this rotating motion. Upon collision, it will fold and the
resulting displacement will be strongly influenced by its elasticity. Figure
3(a) illustrates that when the fiber is rather rigid it bends and is propelled
further from the wall. This effect, comparable to pole vaulting, is
responsible for the fiber finishing its course at a higher height than
initially. Similar pole-vaulting effects have been observed experimentally for
rod-like particles [14, 8, 11, 12]. On the other hand, when the fiber is more
flexible, as shown in Fig. 3(b), it buckles upon contact with the wall and
remains folded, ending its course in a hook shape. As a result of this
folding, the final height of the fiber’s center of mass can be lower than
before tumbling occurred.

Figure 3: Fibers tumbling in a pure shear flow
${\boldsymbol{u}}=\sigma\,y\,\hat{\boldsymbol{x}}$. (a) Rigid fiber bouncing
on the wall. (b) Flexible fiber ending up in a hook position. The red dotted
lines are the center of mass trajectories. (c) Phase diagram showing the
fiber’s final distance to the wall $y_{\rm f}$ as a function of the
flexibility $\mathcal{F}$ and of its initial distance $y_{0}$. The dashed line
separates bouncing and hook configurations. The solid line separates
trajectories ending up further to the wall from those approaching it.

The outcome of a near-wall tumbling event depends upon two parameters: the
local non-dimensional flexibility $\mathcal{F}$ and the initial distance
$y_{0}$ of the fiber from the wall. The occurrence of a pole-vault effect is
intricately determined by these parameters, as demonstrated in the phase
diagram of Fig. 3(c). After a tumbling motion near the wall, fibers can end up
either closer to or further away from the boundary. This depends whether the
parameters are above or below the red solid line that behaves approximately as
$y_{0}\propto-\ell\,\log\mathcal{F}$. The gray dashed line is the separatrix
between the two types of near-wall motion: bounce or hook. It confirms that
the hook shape only occurs for fibers that are flexible enough. Yet, both
types of motion can bring fibers closer or further to the wall. In actual
turbulent boundary layers, the situation is even more complex due to local
fluctuations of the shear and to the presence of near-wall structures.
Furthermore, we consider fibers that are longer than the viscous scale
$\delta_{\nu}$, so their tumbling motion involves scales falling in the
turbulent logarithmic layer. Nevertheless, pole-vaulting during near-wall
tumbling seems a frequent occurrence in turbulent channel flow. Figure 4
serves as an illustrative example, wherein the fiber is bent in three
dimensions but still undergoes a strong enough impulse to migrate away from
the wall. As we will now see, these interactions with the boundaries are the
underlying cause of non-uniform fiber distributions in the channel, ultimately
leading to the enhancement of their transport.

Figure 4: Near-wall tumbling of a fiber of length $\ell=0.25\,h\approx
45\,\delta_{\nu}$ in the turbulent channel flow, leading to an instance of
pole-vaulting. Different colors correspond to various instant of time,
expressed in viscous units.

Building upon prior observations, we now present our findings on the spatial
distribution and velocities of fibers. Figure 5(a) depicts the average
concentration as a function of the wall-normal distance $y^{+}$. A depletion
in the boundary layer is readily apparent, which becomes more pronounced with
an increase in fiber length. As a result, fibers tend to accumulate far from
the boundaries, towards the bulk of the channel flow. The inset (b) of Fig. 5
shows the same data, but in semi-logarithmic coordinates and as a function of
the ratio between the wall-normal distance and the fiber length
$y^{+}/\ell^{+}$. Notably, the data seem to collapse onto a master curve,
suggesting that fiber length is indeed the appropriate scale to characterize
boundary depletion and corroborating the mechanisms described previously.

Figure 5:  Inhomogeneous fiber distributions and velocities. (a) Average
density as a function of the distance to the wall for the various fiber
lengths. (b) Same, plotted as a function of $y^{+}/\ell^{+}$. The dashed line
shows a behavior $\propto\log\ell^{+}$. (c) Average fiber velocity in the
streamwise direction, as a function of the wall-normal coordinate. (d) Fiber
fluxes normalised to that of the fluid. The dashed line represents a log
behavior.

A higher concentration of fibers in the bulk region, where the flow speed is
higher, may suggest that fibers are transported more efficiently by the flow.
However, as the fibers are not expected to follow exactly the fluid due to
their important length, it is essential to examine first their velocity in the
streamwise direction $x$. Figure 5(c) shows this mean velocity, with
statistics conditioned on $X_{2}(s,t)=y$ and therefore not carrying any
information on the concentration. One observes an increase in velocity close
to the boundary and a minor decrease in the bulk, both of which become more
pronounced with fiber length. The increase close to the walls could be
attributed to the concentration of fibers in high-velocity streaks [5, 2, 15].
However, a long fiber in a channel flow can sample different regions of the
flow simultaneously. Thus, segments of the fiber located at different
distances from the wall experience different average fluid velocities. As the
fiber is inextensible, the segments located in the channel center with a
higher mean fluid velocity tend to pull the whole fiber, while those near the
walls with a lower velocity slow it down. This mechanism results in a higher
fiber streamwise velocity than the fluid near the walls and lower in the
center, as observed in Fig. 5(c).

By combining the concentration and velocity measurements to obtain the average
flux in a $(y,z)$-plane, we find that these effects ultimately result in an
increase in fiber transport by several percent, as shown in Fig. 5(d).
Interestingly, the dominant effect on flux modification is the concentration
profile, despite the fact that the fibers are slightly slower than the fluid
in the center of the channel. In particular, we observe that the enhancement
only occurs for fibers whose lengths are within the log-layer of the wall
(typically $\ell^{+}>25$). When $\ell^{+}\in[25,270]$, the increase seems
proportional to $\log\ell^{+}$.

To conclude, we have observed an enhancement in the transport of inertialess
flexible fibers in the case when they can be modeled as slender-bodies and
interact with boundaries through a simple soft-boundary approach. However,
many questions still remain. One such question is the need for more realistic
models for interactions with boundaries to take into account the dissipation
that can occur during collisions due to wall-normal lubrication forces,
plastic deformations during contact and/or tangential friction. These
interactions could affect the transfer of momentum toward the wall-normal
direction and influence the migration of fibers into the bulk. This might
result in an increase of their concentration near the boundary, as observed in
the turbophoresis of inertial particles [29]. This actually brings out a
second question as to whether similar behavior can be observed for fibers
possessing inertia. Another issue that requires further exploration is the
role of two-way coupling [23] in pole-vaulting effects. As fibers deform, they
can store some of the fluid kinetic energy, which they can restore later to
the flow, potentially reducing the significance of pole-vaulting in the
migration of fibers away from the wall. Additionally, in practical
applications such as in the papermaking industry [1], fibers can display
internal properties that fluctuate along their arc-length, such as thickness
or flexibility. This could lead to mixed states upon collisions with
boundaries, with partial hook and pole-vaulting effects, or even to their
break-up into smaller fragments. All these aspects are kept for future work.

We benefited from discussions with Sofìa Allende, who is warmly acknowledged.
We are grateful to the OPAL infrastructure from Université Côte d’Azur for
providing computational resources. This work received support from the UCA-
JEDI Future Investments funded by the French government (grant no.
ANR-15-IDEX-01) and from the Agence Nationale de la Recherche (grant no.
ANR-21-CE30-0040-01).

## References

  * Lundell  _et al._ [2011] F. Lundell, L. D. Söderberg, and P. H. Alfredsson, Annu. Rev. Fluid Mech. 43, 195 (2011).
  * Voth and Soldati [2017] G. A. Voth and A. Soldati, Annu. Rev. Fluid Mech. 49, 249 (2017).
  * Ouchene  _et al._ [2018] R. Ouchene, J. I. Polanco, I. Vinkovic, and S. Simoëns, J. Turbul. 19, 827 (2018).
  * Michel and Arcen [2021] A. Michel and B. Arcen, Phys. Rev. Fluids 6, 114305 (2021).
  * Do-Quang  _et al._ [2014] M. Do-Quang, G. Amberg, G. Brethouwer, and A. V. Johansson, Phys. Rev. E 89, 013006 (2014).
  * Dearing  _et al._ [2013] S. S. Dearing, M. Campolo, A. Capone, and A. Soldati, Exp. Fluids 54, 1419 (2013).
  * Håkansson  _et al._ [2013] K. M. O. Håkansson, M. Kvick, F. Lundell, L. P. Wittberg, and L. D. Söderberg, Exp. Fluids 54, 1555 (2013).
  * Capone  _et al._ [2017] A. Capone, M. Miozzi, and G. P. Romano, Int. J. Multiphas. Flow 94, 189 (2017).
  * Shaik  _et al._ [2020] S. Shaik, S. Kuperman, V. Rinsky, and R. van Hout, Phys. Rev. Fluids 5, 114309 (2020).
  * Alipour  _et al._ [2021] M. Alipour, M. De Paoli, S. Ghaemi, and A. Soldati, J. Fluid Mech. 916, A3 (2021).
  * Alipour  _et al._ [2022] M. Alipour, M. De Paoli, and A. Soldati, J. Fluid Mech. 934, A18 (2022).
  * Baker and Coletti [2022] L. J. Baker and F. Coletti, J. Fluid Mech. 943, A27 (2022).
  * Shaik and van Hout [2023] S. Shaik and R. van Hout, Int. J. Multiphas. Flow 158, 104262 (2023).
  * Moses  _et al._ [2001] K. B. Moses, S. G. Advani, and A. Reinhardt, Rheol. Acta 40, 296 (2001).
  * Dotto and Marchioli [2019] D. Dotto and C. Marchioli, Acta Mech. 230, 597 (2019).
  * Dotto  _et al._ [2020] D. Dotto, A. Soldati, and C. Marchioli, Meccanica 55, 343 (2020).
  * Brouzet  _et al._ [2014] C. Brouzet, G. Verhille, and P. Le Gal, Phys. Rev. Lett. 112, 074501 (2014).
  * Gay  _et al._ [2018] A. Gay, B. Favier, and G. Verhille, Europhys. Lett. 123, 24001 (2018).
  * Picardo  _et al._ [2018] J. R. Picardo, D. Vincenzi, N. Pal, and S. S. Ray, Phys. Rev. Lett. 121, 244501 (2018).
  * Picardo  _et al._ [2020] J. R. Picardo, R. Singh, S. S. Ray, and D. Vincenzi, Phil. Trans. Roy. Soc. A 378, 20190405 (2020).
  * Rosti  _et al._ [2018] M. E. Rosti, A. A. Banaei, L. Brandt, and A. Mazzino, Phys. Rev. Lett. 121, 044501 (2018).
  * Olivieri  _et al._ [2022] S. Olivieri, A. Mazzino, and M. E. Rosti, J. Fluid Mech. 946, A34 (2022).
  * Elgobashi [2006] S. Elgobashi, in  _IUTAM Symposium on Computational Approaches to Multiphase Flow_ , edited by S. Balachandar and A. Prosperetti (Springer, 2006).
  * Allende  _et al._ [2018] S. Allende, C. Henry, and J. Bec, Phys. Rev. Lett. 121, 154501 (2018).
  * Emergent Complexity in Physical Systems ECPS team, Swiss Federal Institute of Technology Lausanne (2018) [EPFL] Emergent Complexity in Physical Systems ECPS team, Swiss Federal Institute of Technology Lausanne (EPFL), Channelflow 2.0 (2018).
  * Tornberg and Shelley [2004] A.-K. Tornberg and M. J. Shelley, J. Comput. Phys. 196, 8 (2004).
  * Parsa  _et al._ [2012] S. Parsa, E. Calzavarini, F. Toschi, and G. A. Voth, Phys. Rev. Lett. 109, 134501 (2012).
  * Pujara  _et al._ [2021] N. Pujara, J.-A. Arguedas-Leiva, C. C. Lalescu, B. Bramas, and M. Wilczek, J. Fluid Mech. 922, R6 (2021).
  * Brandt and Coletti [2022] L. Brandt and F. Coletti, Annu. Rev. Fluid Mech. 54, 159 (2022).


# Bayesian Optimization of Area-based Road Pricing

Renming Liu  DTU Management
Technical University of Denmark
Copenhagen, Denmark
<EMAIL_ADDRESS>Yu Jiang  DTU Management
Technical University of Denmark
Copenhagen, Denmark
<EMAIL_ADDRESS>Carlos Lima Azevedo  DTU Management
Technical University of Denmark
Copenhagen, Denmark
<EMAIL_ADDRESS>
###### Abstract

This study presents a Bayesian Optimization framework for area- and distance-
based time-of-day pricing (TODP) for urban networks. The road pricing
optimization problem can reach high level of complexity depending on the
pricing scheme considered, its associated detailed network properties and the
affected heterogeneous demand features. We consider heterogeneous travellers
with individual-specific trip attributes and departure-time choice parameters
together with a Macroscopic Fundamental Diagram (MFD) model for the urban
network. Its mathematical formulation is presented and an agent-based
simulation framework is constructed as evaluation function for the TODP
optimization problem. The latter becomes highly nonlinear and relying on an
expensive-to-evaluate objective function. We then present and test a Bayesian
Optimization approach to compute different time-of-day pricing schemes by
maximizing social welfare. Our proposed method learns the relationship between
the prices and welfare within a few iterations and is able to find good
solutions even in scenarios with high dimensionality in the decision variables
space, setting a path for complexity reduction in more realistic road pricing
optimization problems. Furthermore and as expected, the simulation results
show that TODP improves the social welfare against the no-pricing case.

###### Index Terms:

Bayesian Optimization, Road pricing, Macroscopic Fundamental Diagram, Demand
Management, Day-to-day Dynamics, Machine Learning

††publicationid: pubid: 978-1-7281-8995-6/21/$31.00 © 2021 IEEE

##  I Introduction

As one of the most promising traffic demand management policies for mitigating
traffic congestion [1], congestion pricing has been investigated extensively
over a century. Various pricing models and solution methods have been proposed
in theory and practice [2, 3, 4, 5]. These are often categorized based on the
type of infrastructure they are applied to [6]:

##### Facility-based schemes

Under such schemes tolls are imposed at a single or multiple points on a
specific facility such as a road link, bridge, tunnel, or levied on a part of
the facility, e.g.: High Occupancy Vehicle (HOV) lanes.

##### Cordon-based schemes

Travelers crossing a cordon in the inbound and/or outbound directions are
charged, while those travelling inside the cordoned area are not.

##### Area-based schemes

Travelers are required to pay a toll when entering, exiting or traveling
within a defined area.

##### Distance-based schemes

The charges are associated with the trip length linearly or nonlinearly [7],
mainly to enhance fairness and efficiency compared to the previous purely
access-based controls [8].

Compared to facility-based schemes, area-based schemes are able to intercept
more trips and less sensitive to traffic diversion [9]. It could simplify the
complexity of the problem in both theory and practice as it avoids link-based
pricing and, in some cases, does not require knowing the within-network
dynamics, link attributes, and detailed information of path selection for
different origins and destinations. Yet these simplifications also encompass
limitations in efficiency and applicability compared to more disaggregate
controls [6].

Despite the differences in the above schemes, finding an optimal toll value
for practice is often associated with a computationally expensive objective
function, high dimensional decision variables (due to a large number of
possible price values of a step-toll [10], dynamic link-specific tolls for
large networks [11] or realistic demand and supply modelling features [12] in
the optimization evaluation function). The optimization problems often do not
have a closed form solution that can be solved analytically. Alternatively,
simulation-based optimization approaches are useful to handle such toll
optimization problem with an expensive-to-evaluate objective function and high
non-linearity, as these approaches only require the paired data of decision
variables and objective function values to search for optimal decision
variables [13, 12, 14]. In the existing applications of simulation-based
optimization on pricing toll scheme design, there are two different
approaches: feedback control and surrogate-based optimization. [15] proposed a
discrete integral controller to adjust the cordon-based time-dependent charge
rates so that the accumulation of the network does not exceed the flow-
maximizing value. [16] further developed a proportional-integral controller to
iteratively control area-based toll rates, and the results demonstrated that
the later controller outperformed the former one due to higher flexibility.
Recently, [12] used feature-variant clustering methods for toll-area
definition under distance- and area-based schemes.

On the other hand, surrogate-based optimization aims to approximate the map
from the decision variables to the objective function values within a few
iterations. For example, [17] constructed the surrogate model using a radial
basis function to optimize the link-specific tolls for networks, and showed
that the proposed method converged faster than the genetic algorithm. [14]
compared different surrogate models for the link toll optimization problem and
found that kriging (or Gaussian process regression) performed best. [18]
applied the surrogate model with expected improvement sampling to optimized
the step toll parameters to control the accumulation of the network.

In this study, we exploit the potential of using Bayesian Optimization (BO)
for determining the optimum price in distance- and area-based pricing schemes.
We focus on the morning commute problem and on the design of time-of-day
pricing (TODP) profiles for networks with heterogeneous travelers. TODP is
usually preferable rather than flat tolls, as it can internalize the dynamic
congestion costs that travelers impose on each other [6], reduce the social
costs of trips in non-congested periods and help improving overall network
performance [19]. We also extend other area-based pricing optimization
frameworks to the objective of social welfare maximization instead of purely
network performance measures. We use a Gaussian (mixture) function to
parameterize the TODP scheme for the purpose of reducing the dimension of the
decision variables.

Our hypothesis is that BO can handle the stochasticity from relatively complex
simulation-based evaluation functions and find good solutions within few
evaluations. BO’s implicit surrogate-based method plays a role in mapping the
decision variables (i.e., the pricing toll scheme) to the objective function
(i.e., the social welfare) and its inherited uncertainty from the demand and
network simulation models. For the latter, we rely on the recent developments
in (network) Macroscopic Fundamental Diagram (MFD) modeling founded by [20,
21] for capturing the dynamics of travel speed, travel production and network
accumulation under different pricing setups. Furthermore, we employ the trip-
based MFD formulation proposed in [22, 23, 24] for allowing for the inclusion
of heterogeneous trip lengths in both the control feature and the network
performance function.

##  II Methodology

###  II-A Trip-based MFD Model

As defined in [23, 24], the trip-based MFD considers the trip distance
traveled by traveler $i$ as a function of trip departure time $t_{i}^{dep}$
and arrival time, $t_{i}^{dep}+T_{i}(t_{i}^{dep})$, within the network. It
models it as the integration of the network travel speed over the travel time
period, i.e.,

$L_{i}=\int_{t_{i}^{dep}}^{t_{i}^{dep}+T_{i}(t_{i}^{dep})}V(n(t))dt$ (1)

where $L_{i}$ is the trip length of traveler $i$, $V(\cdot)$ represents the
travel speed and $n(t)$ is the accumulation of the network at time $t$. In
this work, the area-based network is considered as a single-reservoir in which
the traffic congestion is assumed to be distributed in space with small
heterogeneity[20, 21]. The traditional MFD assumption that the travel speed
$V(n(t))$ is the same for all travelers in the network at time $t$ holds. An
event-based simulation is then employed to represent the trip-based MFD
process, as proposed by [24, 25]:

_Step 1._ Initialization: Specify $t_{i}^{dep}$, $L_{i}$, $V(n)$ and the total
number of travelers $N$; set $n=0$, calculate the initially estimated arrival
time for all travelers $\forall i,1...N$ by $L_{i}/V(0)$ using (1). Note that
$n(t)$ is simplified as $n$.

_Step 2._ Construct the event list in the order of time, which consists of $N$
departures and $N$ arrivals.

_Step 3._ Calculate the remaining trip distance for travelers who have entered
the network and not finished their trip by
$L_{i}=L_{i}-V(n)\cdot(t_{j}-t_{j-1})$

_Step 4._ For every event, if the next closest event is another traveler
$i^{\prime}$ departure, set $n=n+1$ and remove this event from the event list;
otherwise, set $n=n-1$, output the experienced travel time of $i$,
$T_{i}(t_{i}^{dep})$, and remove this event from the list.

_Step 5._ Update the current average traveling speed $V(n)$; update the
estimated arrival time for the rest of travelers in the network by using eq.
(1); and sort the event list in the order of time again.

_Step 6._ Go to _Step 3_ until all travelers have ended their trips.

Note that the heterogeneous trip length can be accommodated by the trip-based
MFD model as it handles the travel distances of all travelers separately.

###  II-B Day-to-Day Dynamics Under TODP Strategy

In a morning commute problem, travelers make departure time choices within a
predetermined time window based on perceived travel costs [26], often modeled
via the discrete choice model [27]. Let $t$ represent $t_{i}^{dep}$ for
simplicity, and $C_{i,d}(t)$ be the observable monetary utility of traveler
$i$ departing at time $t$ on day $d$, the probability of choosing departure
time $t$ on day $d$, $Pr_{i,d}(t)$, is written as,

$Pr_{i,d}(t)=\frac{\exp\big{(}\mu\cdot C_{i,d}(t)\big{)}}{\sum_{s\in
TW_{i}}\exp\big{(}\mu\cdot C_{i,d}(s)\big{)}}$ (2)

where $\mu$ is the scale parameter and $TW_{i}$ is the time window of traveler
$i$. The utility of the departure time choice for traveler $i$ then consists
of the observable utility and a random term which captures the unobserved
utility, written as follows,

$U_{i,d}(t)=C_{i,d}(t)+\epsilon_{i}$ (3)

where $\epsilon_{i}$ is assumed to be i.i.d., following the extreme value
distribution, Gumbel distribution.

Under a given pricing scheme, the monetary metric utility includes three
components, the travel time cost, schedule delay cost (the difference between
realized arrival time and desired arrival time) and the pricing toll payment,

$\displaystyle c_{i,d}(t)=$
$\displaystyle-\theta_{i}\cdot\Big{[}T_{i,d}(t)+\delta_{i}\cdot
SDE_{i}\cdot\big{(}T_{i}^{*}-t-T_{i,d}(t)\big{)}$ (4)
$\displaystyle+(1-\delta_{i})\cdot
SDL_{i}\cdot\big{(}t+T_{i,d}(t)-T_{i}^{*}\big{)}\Big{]}$ $\displaystyle-
Toll(t)\cdot L_{i}\cdot w$ $\displaystyle=$ $\displaystyle-\theta_{i}\cdot
tc_{i,d}(t)-Toll(t)\cdot L_{i}\cdot w$

The first term represents the sum of the travel time and schedule delay costs,
where $\theta_{i}$ is the value of time for traveler $i$, $T_{i,d}(t)$ is the
travel time for traveler $i$ on day $d$ departing at time $t$, $SDE_{i}$ and
$SDL_{i}$ are the schedule deviation penalty parameters for early and late
arrival for traveler $i$, and $\delta_{i}$ is a binary variable that equals 1
if $i$ arrives early and 0 otherwise. It is worth to note that the simulation
in Section II-A can only measure the travel time for the departure time chosen
by traveler $i$. [28] proposed to estimate the travel time for the not-chosen
departure time by assuming there are fictional travelers departing at these
departure times without being counted when update the number of travelers in
the network. The second term represents the trip length-specific pricing toll
payment, where $Toll(\cdot)$ is a time-dependent toll function and $w$ is a
constant used to scale down the magnitude of the trip length. In this study,
the toll function is parameterized by a Gaussian function,
$Toll(t|A,\xi,\sigma)=A\times e^{\frac{-(t-\xi)^{2}}{2\sigma^{2}}}$, which is
controlled by three parameters, mean $\xi$, variance $\sigma$ and amplitude
$A$. Without loss of generality, this toll function can be extended to a
linear Gaussian mixture function to provide asymmetry and flexibility. This
issue will be further investigated with numerical experiments in Section IV.

At the end of each day $d$, travelers update their perception of the
observable monetary utility for day $d+1$ by the linear weighted sum of the
initially perceived generalized costs on day $d$ and the experienced (chosen
alternative) or estimated (unchosen alternatives) utility on day $d$, as
follows:

$C_{i,d+1}(t)=\omega\cdot C_{i,d}(t)+(1-\omega)\cdot c_{i,d}(t)$ (5)

where the weight coefficient $0<\omega<1$ is the learning parameter, which is
here assumed to be identical among travelers in this study. This day-to-day
evolution process is expected to converge to an equilibrium state that the
perceived and the experienced (or estimated) utility remain the same
afterwards and there is no significant changes in departure time choices. See
Section IV-A for details on the equilibrium process.

It is relevant to note that the aforementioned day-to-day dynamics is
essentially a stochastic simulation, of which randomness comes from different
sources such as random choice parameters, random trip attributes or even the
probabilistic day-to-day departure time decision process.

##  III Simulation-based Optimization Framework

###  III-A Objective function

To evaluate the performance of the proposed pricing scheme, we compare the
social welfare per capita $W$ at the equilibrium states of scenarios with and
without pricing. For the no toll equilibrium (or NTE), combined with (3) and
(4), the social welfare per capita $W_{NTE}$ is computed as,

$\displaystyle W_{NTE}=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}U_{i,*}(t_{i,*}^{dep})$ (6)
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\Big{[}-\theta_{i}\cdot
tc_{i,*}(t_{i,*}^{dep})+\epsilon_{i}(t_{i,*}^{dep})\Big{]}$

where the subscript ’$*$’ denotes the equilibrium state.

For the scenario with the TODP scheme, the social welfare $W_{TODP}$ is the
consumer surplus (CS) plus the regulator revenue (RR) from charging travelers,
which can be written as follows,

$\displaystyle W_{TODP}=$ $\displaystyle CS+RR$ (7) $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\Big{(}c_{i,*}(t_{i,*}^{dep})+\epsilon_{i}(t_{i,*}^{dep})+Toll(t_{i,*}^{dep})\cdot
L_{i}\cdot w\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\Big{[}-\theta_{i}\cdot
tc_{i,*}(t_{i,*}^{dep})+\epsilon_{i}(t_{i,*}^{dep})\Big{]}$

Note that the welfare computations in both scenarios are equivalent and
consist of travel time cost, schedule delay cost and the random component.

###  III-B Bayesian Optimization

The BO has two core iterative steps [29]: 1) Update the Gaussian Process (GP)
that provides a posterior probability distribution of the objective function;
and 2) Determine where to sample the next decision variables by optimizing a
given acquisition function.

####  III-B1 Gaussian Process

With a mean function $\mu(\bm{x})$ and covariance function (or kernel)
$k(\bm{x},\bm{x^{\prime}})$, where $\bm{x}$ represents the decision variables,
which are the parameters of the pricing toll scheme, the GP regression can
approximate the objective function using the historical observations, by
taking the distribution of objective function values and decision variables to
be a multivariate Gaussian distribution.

Suppose we have $m$ observed objective values as
$\mathcal{D}_{m}=\\{\bm{x}_{1:m},W_{1:m}\\}$, where
$\bm{x}_{1:m}=[\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{m}]^{T}$ are the decision
variables and $W_{1:m}=[W_{1},W_{2},\dots,W_{m}]^{T}$ are the corresponding
objective values, and wish to infer the objective function value at some new
decision variables $\bm{x}_{m+1}$. We can first derive the joint distribution
of $W_{1:m}$ and $W_{m+1}$ as follows, wherein the mean function is
$\mu(\bm{\cdot})$ is assumed to be zero function for simplicity [30]:

$\left[\begin{array}[]{c}W_{1:m}\\\
W_{m+1}\end{array}\right]~{}\sim~{}\mathcal{N}\Bigg{(}\bm{0},~{}\left[\begin{array}[]{cc}\bm{K}&\bm{k}\\\
\bm{k}^{T}&k(x_{m+1},x_{m+1})\end{array}\right]\Bigg{)}$ (8)

where
$\bm{k}=[k(\bm{x}_{m+1},\bm{x}_{1}),k(\bm{x}_{m+1},\bm{x}_{2}),\dots,k(\bm{x}_{m+1},\bm{x}_{m})]^{T}$,
and $\bm{K}$ is the covariance matrix with entries
$\bm{K}_{i,j}=k(\bm{x}_{i},\bm{x}_{j})$ for $i,~{}j\in\\{1,2,\dots,m\\}$.

We can then compute the conditional distribution of $W_{m+1}$ using Bayes’
rule:

$W_{m+1}|W_{1:m}~{}\sim~{}\mathcal{N}\Big{(}\mu(\bm{x}_{m+1}),~{}\sigma^{2}(\bm{x}_{m+1})\Big{)}$
(9)

where $\mu(\bm{x}_{m+1})=\bm{k}^{T}\bm{K}^{-1}W_{1:m}$ and
$\sigma^{2}(\bm{x}_{m+1})=k(x_{m+1},x_{m+1})-\bm{k}^{T}\bm{K}^{-1}\bm{k}$.

The kernel encodes the correlation between two sets of decision variables via
a parametric structure, which is required to be a positive semi-definite
function [29]. In this work, we adopt the commonly used Matern kernel [31].

####  III-B2 Acquisition Function

Based on the inferred mean and variance of the objective function values,
acquisition function determines the next decision variables with the maximum
utility considering the trade-off between exploration and exploitation. This
study uses a popular acquisition function, the upper confidence bound (UCB)
[32], as follows:

$\alpha_{UCB}(\bm{x};\beta)=\mu(\bm{x})+\beta\sigma(\bm{x})$ (10)

where $\beta$ determines the balance between exploration and exploitation. A
larger $\beta$ leads to more exploration.

####  III-B3 Space-filling Experimental Design

To enhance the solution quality and optimization efficiency, it is useful to
start with initial sets of decision variables that has a large coverage of the
feasible region. The Latin Hypercube Sampling (LHS) is used in this study for
the initial space-filling design, which can generate independent sampled sets
of decision variables without overlap, providing a better representative of
the real variability than random generation and Monte Carlo methods.

####  III-B4 High dimensional BO Using Dropout

Let $D=\text{dim}(\bm{x})$ represent the dimension of the decision variables.
When $D$ increases, the size of the search space grows exponentially, the
maximization of the acquisition function becomes more complex and the
estimation of the GP is less accurate [33]. All these three issues impact the
performance of BO in high dimension. [34] tackles these challenges by using
the Dropout strategy from deep neural networks, that is, BO steps are
implemented on the subspace of different $d$ out of $D$. Its performance is
tested in Section IV via numerical experiments.

##  IV Numerical Experiments

This section presents the results of (1) the convergence of the day-to-day
departure time choices with and without the TODP scheme; (2) the performance
of the BO; and (3) the comparison between the optimized TODP against the NTE
case. The numerical settings are presented in Table I.

TABLE I: Numerical settings Parameters | Specification
---|---
Demand |  $N=3700$ [traveler]
Trip length |  $L_{i}=4600+\mathcal{N}(0,(0.2*4600)^{2})$ [m], $L_{i}>0$
Scale factor of trip length | $w=2\times 10^{-4}$
Value of time |  $\theta_{i}=1.1$ [DKK/min] [35]
Schedule deviation penalty | $\begin{bmatrix}SDE_{i}\\\ SDL_{i}\end{bmatrix}=\begin{bmatrix}0.5\\\ 4\end{bmatrix}+\mathcal{N}\Bigg{(}\begin{bmatrix}0.05^{2}&0.1^{2}\\\ 0.1^{2}&0.4^{2}\end{bmatrix}\Bigg{)}$
|  $SDE_{i}\in[0.3,0.7]$, $SDL_{i}\in[2.5,5.5]$
Network capacity |  $n_{jam}=4500$ [vehicle]
Free flow speed |  $v_{f}=9.78$ [m/s]
Speed function |  $V(n)=v_{f}(1-\frac{n}{n_{jam}})^{2}$ [m/s]
Learning parameter | $\omega=0.7$
Toll profile function | $Toll(t|A,\xi,\sigma)=A\times e^{\frac{-(t-\xi)^{2}}{2\sigma^{2}}}$

In this study, we assume the capacity of the single-reservoir network as 4500
travelers, using speed function in [36] and other parameters adopted in [25]
and [37]. We can derive the flow-maximizing accumulation value (or critical
value) $n_{cr}=1500$ [travelers]. To test the performance of the pricing
scheme under a high congestion scenario, the demand is set at 3700 travelers
so that the peak accumulation value at the no toll equilibrium state exceeds
the critical value (see Section IV-A for details). Besides, we form the
heterogeneity of travelers by drawing their trip lengths and schedule
deviation penalties from three truncated Gaussian distributions shown in Table
I.

###  IV-A Day-to-day evolution process

In this section, we illustrate that with and without the TODP scheme, the
proposed day-to-day evolution can reach the equilibrium states where the
departure time choices of all travelers keep unchanged. Besides, the perceived
generalized cost vector should equal to the experienced cost vector at the
equilibrium, i.e., $\bm{C}_{*}=\bm{c}_{*}$. Hence we compute the inconsistency
between $\bm{C}_{d}$ and $\bm{c}_{d}$ as
$\parallel\bm{C}_{d}-\bm{c}_{d}\parallel_{1}/N$ to indicate the degree of
convergence.

Fig. 1(a) presents the convergence process of the inconsistency in the no toll
scenario, which comes close to 0 and stays stable after 25 days. Fig. 1(b)-(c)
show the evolution process of the average consumer surplus and social welfare
per capita, respectively. Curves of these two plots are the same since the
social welfare equals to the consumer surplus when there is no pricing. Fig.
1(d) illustrate the within-day evolution of accumulation on different days, it
can be found that the accumulation profile becomes stable after 25 days, and
the peak accumulation exceeds the critical value $n_{cr}$. The observations
from Fig. 1 imply that the equilibrium state of the day-to-day evolution is
reached. Note that the inconsistency, consumer surplus and social welfare are
plotted from day 1 as travelers do not have the perception of cost on day 0.

Figure 1: The evolution process of the no toll scenario.

To testify the convergence, we impose a TODP with parameters $A=11$, $\xi=80$
and $\sigma=18$ of day-to-day dynamics. In this case, the toll scheme is not
optimized but arbitrarily given. As shown in Fig. 2(a)-(c), the day-to-day
evolution under a given TODP scheme also converges to the equilibrium after 25
days. In addition to accumulation profiles, Fig. 2(d) also plots the TODP toll
profile in a dashed grey line. It is found that travelers change their
departure time to avoid the period with higher charge rate.

Figure 2: The evolution process under a TODP scheme.

###  IV-B Algorithm Performance

In this section, we present the application of the simulation-based
optimization framework to optimize the TODP parameters and evaluate the
performance of BO. We first investigate the influence of the shape of TODP
toll profile on social welfare by parametrizing the TODP with a mixture
function, which vary from one (1G) to six (6G) Gaussian functions.

In this test, the domains of the TODP function parameters (decision variables)
are set as $A\in[4,30]$ (unit: DKK/m), $\xi\in[30,90]$ and $\sigma\in[10,50]$.
To be fair, we generate 30 initial sampled sets of decision variables to fill
the space for each case, and the number of function evaluations is fixed as
60. Fig. 3 illustrates the evolution of the best social welfare with the
number of function evaluations for BO using six different TODP functions. It
can be found that BO has a relatively stable performance when the dimension of
the decision variables is lower or equal to 12, and the obtained toll profiles
have similar shapes. However, for the cases of 5G and 6G (15 and 18
parameters, respectively), standard BO fails to reach the maximal value
obtained in lower dimension cases. The best solution found by BO is in the
case of 1G, where $A=26.2$, $\xi=67.1$ and $\sigma=28.8$ with
$W_{TODP}=-26.76$ [DKK].

Figure 3: BO’s performance with different dimensions of variables.

Moreover, we adopt the BO dropout method to enhance the performance of BO for
the case with high dimensions, e.g., for the case of 6G with 18 parameters.
First, four cases with $d=1,3,5,7$ are tested to find out the effect of $d$.
Since the dropout method chooses the decision variables in a random way, we
propose two promising variable choosing strategies for further improvement
with $d=5$: In strategy 1 (S1), we categorize the 18 variables into 3 groups
by physical meanings, e.g., all the $A$s are in one group, and at least one
variable of each group will be chosen in each iteration of BO; In strategy 2
(S2), we categorize the 18 variables into 6 groups by functions, e.g.,
$A_{1}$, $\xi_{1}$ and $\sigma_{1}$ are in a group, we randomly choose 5 out
of the 6 groups and draw one variable in each chosen group. We run each case
four times and summarize the average objective values and standard deviations
in Table II. The results demonstrate that all dropout cases reach higher
objective values with smaller standard deviations compared to the standard BO.
Besides, the proposed two strategies enhance the performance of the pure
random dropout method with $d=5$, and strategy 2 even outperforms the random
dropout method in all cases. This implies that performance improvement might
be realized by developing case-specific proper rules on choosing variables.
Nevertheless, more experiments are needed to investigate the influence of $d$
and dropout strategies in terms of uncertainty.

TABLE II: Average Performances of Different Dropout Methods Method | BO ($6G$) | $d=1$ | $d=3$ | $d=5$ | $d=7$ | S1 | S2
---|---|---|---|---|---|---|---
Mean | 27.00 | 26.81 | 26.82 | 26.85 | 26.85 | 26.82 | 26.79
Std.dev | 0.17 | 0.04 | 0.07 | 0.05 | 0.08 | 0.05 | 0.01

###  IV-C Optimization results

In this subsection we present the optimization results of the TODP scheme,
wherein the toll profile is represented by a single Gaussian function. Fig.
4(a) illustrates the evolution of social welfare per capita and Fig. 4(b)
plots the variations of accumulation and the optimized toll profile. Compared
to the no toll case, the peak accumulation is reduced from 1792 to 1281
[traveler], leading to 16.0% improvement in the social welfare per capita from
-31.86 [DKK] to -26.76 [DKK]. Moreover, the average travel time cost is
reduced by 3.4% and the average schedule deviation cost is increased by 37.2%.

Figure 4: The evolution process under optimized TODP.

##  V Conclusion

This paper focuses on the dynamics of TODP, under heterogeneous demand and a
trip-based MFD network simulation model, and presents a simulation-based BO
framework for optimizing the pricing scheme in terms of social welfare.

The performed numerical results show that BO performs worse with higher
dimensions in the decision variable space, as the needed number of function
evaluations increases exponentially to find the optimum. When testing the
dropout method proposed by [34], we were able to reach higher dimension
search, which is quite relevant for the problem at stake when implementing
detailed pricing schemes. We then propose a model-informed approach to
formulate problem-specific dropout strategies, and showcase how selecting
related toll profile features during the dropout process outperforms the
existing random dropout method in terms of optimization and stability. Such
insights encourages the testing of our proposed method to other practical but
possibly high-dimension TODP schemes, such as step tolls or multiple area
networks. Finally, our optimized TODP is compared against the no-control
strategy, showing clear benefits in terms of social welfare and network
performance.

## Acknowledgment

This research was carried out under the NEMESYS project funded by the DTU-NTU
Alliance.

## References

  * [1] T. Langmyhr, “Understanding innovation: the case of road pricing,” Transport Reviews, vol. 19, no. 3, pp. 255–271, 1999.
  * [2] H. Yang and H.-J. Huang, Mathematical and economic theory of road pricing.  Oxford: Elsevier, 2005.
  * [3] R. Lindsey, “Do economists reach a conclusion?,” Econ Journal Watch, vol. 3, no. 2, pp. 292–379, 2006.
  * [4] D. Vonk Noordegraaf, J. A. Annema, and B. van Wee, “Policy implementation lessons from six road pricing cases,” Transportation Research Part A: Policy and Practice, vol. 59, pp. 172 – 191, 2014.
  * [5] Z. Gu, Z. Liu, Q. Cheng, and M. Saberi, “Congestion pricing practices and public acceptance: A review of evidence,” Case Studies on Transport Policy, vol. 6, no. 1, pp. 94–101, 2018.
  * [6] A. de Palma and R. Lindsey, “Traffic congestion pricing methodologies and technologies,” Transportation Research Part C: Emerging Technologies, vol. 19, no. 6, pp. 1377–1399, 2011.
  * [7] Q. Meng, Z. Liu, and S. Wang, “Optimal distance tolls under congestion pricing and continuously distributed value of time,” Transportation Research Part E: Logistics and Transportation Review, vol. 48, no. 5, pp. 937–957, 2012\.
  * [8] C. F. Daganzo and L. J. Lehe, “Distance-dependent congestion pricing for downtown zones,” Transportation Research Part B: Methodological, vol. 75, pp. 89–99, 2015.
  * [9] P. Olszewski and L. Xie, “Modelling the effects of road pricing on traffic in singapore,” Transportation Research Part A: Policy and Practice, vol. 39, no. 7-9, pp. 755–772, 2005.
  * [10] C. Robin Lindsey, V. A. van den Berg, and E. T. Verhoef, “Step tolling with bottleneck queuing congestion,” Journal of Urban Economics, vol. 72, no. 1, pp. 46 – 59, 2012.
  * [11] D. Joksimovic, M. C. J. Bliemer, and P. H. L. Bovy, “Optimal toll design problem in dynamic traffic networks with joint route and departure time choice,” Transportation Research Record, vol. 1923, no. 1, pp. 61–72, 2005\.
  * [12] A. F. Lentzakis, R. Seshadri, A. Akkinepally, V.-A. Vu, and M. Ben-Akiva, “Hierarchical density-based clustering methods for tolling zone definition and their impact on distance-based toll optimization,” Transportation Research Part C: Emerging Technologies, vol. 118, p. 102685, 2020.
  * [13] S. Amaran, N. V. Sahinidis, B. Sharda, and S. J. Bury, “Simulation optimization: a review of algorithms and applications,” Annals of Operations Research, vol. 240, no. 1, pp. 351–380, 2016.
  * [14] X. Chen, L. Zhang, X. He, C. Xiong, and Z. Li, “Surrogate-based optimization of expensive-to-evaluate objective for optimal highway toll charges in transportation network,” Computer-Aided Civil and Infrastructure Engineering, vol. 29, no. 5, pp. 359–381, 2014.
  * [15] N. Zheng, R. A. Waraich, K. W. Axhausen, and N. Geroliminis, “A dynamic cordon pricing scheme combining the macroscopic fundamental diagram and an agent-based traffic model,” Transportation Research Part A: Policy and Practice, vol. 46, no. 8, pp. 1291–1303, 2012.
  * [16] N. Zheng, G. Rérat, and N. Geroliminis, “Time-dependent area-based pricing for multimodal systems with heterogeneous users in an agent-based environment,” Transportation Research Part C: Emerging Technologies, vol. 62, pp. 133–148, 2016.
  * [17] J. Y. Chow and A. C. Regan, “A surrogate-based multiobjective metaheuristic and network degradation simulation model for robust toll pricing,” Optimization and Engineering, vol. 15, no. 1, pp. 137–165, 2014.
  * [18] Z. Gu, S. T. Waller, and M. Saberi, “Surrogate-based toll optimization in a large-scale heterogeneously congested network,” Computer-Aided Civil and Infrastructure Engineering, vol. 34, no. 8, pp. 638–653, 2019.
  * [19] A. De Palma, R. Lindsey, E. Quinet, et al., “Time-varying road pricing and choice of toll locations,” Road pricing: Theory and evidence, research in transportation economics, vol. 9, pp. 107–131, 2004.
  * [20] C. F. Daganzo, “Urban gridlock: Macroscopic modeling and mitigation approaches,” Transportation Research Part B: Methodological, vol. 41, no. 1, pp. 49–62, 2007.
  * [21] N. Geroliminis, C. F. Daganzo, et al., “Macroscopic modeling of traffic in cities,” in Transportation Research Board 86th Annual Meeting, no. 07-0413, No. 07-0413, 2007.
  * [22] M. Fosgerau, “Congestion in the bathtub,” Economics of Transportation, vol. 4, no. 4, pp. 241–255, 2015.
  * [23] R. Lamotte and N. Geroliminis, “The morning commute in urban areas: Insights from theory and simulation,” tech. rep., Transportation Research Board 95th Annual Meeting, 2016.
  * [24] G. Mariotte, L. Leclercq, and J. A. Laval, “Macroscopic urban dynamics: Analytical and numerical comparisons of existing models,” Transportation Research Part B: Methodological, vol. 101, pp. 245–267, 2017\.
  * [25] M. Yildirimoglu and M. Ramezani, “Demand management with limited cooperation among travellers: A doubly dynamic approach,” Transportation Research Part B: Methodological, vol. 132, pp. 267–284, 2020.
  * [26] W. S. Vickrey, “Congestion theory and transport investment,” The American Economic Review, vol. 59, no. 2, pp. 251–260, 1969.
  * [27] M. Ben-Akiva, M. Cyna, and A. De Palma, “Dynamic model of peak period congestion,” Transportation Research Part B: Methodological, vol. 18, no. 4-5, pp. 339–355, 1984.
  * [28] R. Lamotte and N. Geroliminis, “Dynamic traffic modeling: Approximating the equi-librium for peak periods in urban areas,” 2015.
  * [29] P. I. Frazier, “A tutorial on bayesian optimization,” arXiv preprint arXiv:1807.02811, 2018.
  * [30] C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, vol. 2.  MIT press Cambridge, MA, 2006.
  * [31] B. Matérn, Spatial variation, vol. 36.  Springer Science & Business Media, 2013.
  * [32] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger, “Gaussian process optimization in the bandit setting: No regret and experimental design,” arXiv preprint arXiv:0912.3995, 2009.
  * [33] K. Kandasamy, J. Schneider, and B. Póczos, “High dimensional bayesian optimisation and bandits via additive models,” in International conference on machine learning, pp. 295–304, 2015.
  * [34] C. Li, S. Gupta, S. Rana, V. Nguyen, S. Venkatesh, and A. Shilton, “High dimensional bayesian optimization using dropout,” in Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17, pp. 2096–2102, 2017.
  * [35] M. Fosgerau, K. Hjorth, and S. V. Lyk-Jensen, “The danish value of time study,” 2007.
  * [36] R. Lamotte and N. Geroliminis, “The morning commute in urban areas with heterogeneous trip lengths,” Transportation Research Part B: Methodological, vol. 117, pp. 794–810, 2018.
  * [37] R. Liu, S. Chen, Y. Jiang, R. Seshadri, M. E. Ben-Akiva, and C. L. Azevedo, “Managing network congestion with tradable credit scheme: a trip-based mfd approach,” 2020.