# Coalitional Bargaining via Reinforcement Learning: An Application to
Collaborative Vehicle Routing

Stephen Mak
University of Cambridge
&Liming Xu
University of Cambridge
&Tim Pearce
Tsinghua University
&Michael Ostroumov
Value Chain Lab
&Alexandra Brintrup
University of Cambridge

###### Abstract

Collaborative Vehicle Routing is where delivery companies cooperate by sharing
their delivery information and performing delivery requests on behalf of each
other. This achieves economies of scale and thus reduces cost, greenhouse gas
emissions, and road congestion. But which company should partner with whom,
and how much should each company be compensated? Traditional game theoretic
solution concepts, such as the Shapley value or nucleolus, are difficult to
calculate for the real-world problem of Collaborative Vehicle Routing due to
the characteristic function scaling exponentially with the number of agents.
This would require solving the Vehicle Routing Problem (an NP-Hard problem) an
exponential number of times. We therefore propose to model this problem as a
coalitional bargaining game where – crucially – agents are _not_ given access
to the characteristic function. Instead, we _implicitly_ reason about the
characteristic function, and thus eliminate the need to evaluate the VRP an
exponential number of times – we only need to evaluate it once. Our
contribution is that our decentralised approach is both scalable and considers
the self-interested nature of companies. The agents learn using a modified
Independent Proximal Policy Optimisation. Our RL agents outperform a strong
heuristic bot. The agents correctly identify the optimal coalitions 79% of the
time with an average optimality gap of 4.2% and reduction in run-time of 62%.

##  1 Introduction

Heavy Goods Vehicles (HGVs) in the UK contributed 4.3% of the UK’s _total_
greenhouse gas emissions in 2019 [39]. Furthermore, HGVs are utilised
inefficiently at 61% of their total weight capacity. Moreover, 30% of the
distance travelled is empty [40, RFS0125].

Collaborative Vehicle Routing (CVR) has been proposed to improve HGV
utilisation. Here, delivery companies, or _carriers_ , need to share their
delivery information in order to achieve economies of scale. If companies
agree to work together, they are said to be in a _coalition_. The case where
all companies join the same coalition is called the _grand coalition_. As a
result of improved utilisation, total travel costs can be reduced resulting in
a _collaboration gain_. The remaining question is then how to allocate this
collaboration gain in a reasonable manner so that carriers are incentivised to
form coalitions. An example of Collaborative Vehicle Routing is given in
Figure 1.

(a) Pre-collaboration
Total Cost = 3.35

(b) Post-collaboration (with grand coalition, $\\{1,2,3\\}$)
Total Cost = 2.47

Figure 1: Three agents (denoted by colours) before and after collaboration.
Squares denote depots. Crosses denote customer locations. Node indices
(arbitrary) are denoted in black, with costs given in their respective
colours. The _collaboration gain_ is defined as the difference in social
welfare before and after collaboration. In Figure 1b, Agents 1, 2 and 3 decide
to cooperate forming the grand coalition and reduces the system’s total cost
by 0.88 or 26%. This results in a _collaboration gain per capita_ (assuming
agents split the gain equally) of 0.29. For detailed calculations, see
Appendix A.

Prior literature suggests that collaborative routing can reduce costs by
around 20-30% with further reductions to greenhouse gas emissions and road
congestion [8, 30, 16, 17, 14]. However, real-world adoption remains limited,
with only a few companies participating, according to a number of case-studies
conducted [8, 3, 10, 19]. One reason behind this is the complexity of fair
_gain sharing_ amongst a larger number of companies.

Whilst gain sharing has been studied in collaborative routing using
cooperative game theory [19], the solution concepts typically assumes that the
characteristic function is given. For a set of $n$ agents,
$N=\\{1,\dots,n\\}$, the characteristic function
$v:\mathbf{2}^{N}\to\mathbb{R}$ assigning a _value_ , or in our case
_collaboration gain_ , for every possible coalition that could be formed. Note
that there exists $\mathcal{O}(2^{n})$ possible coalitions. This is
intractable for settings with more than a few agents, because evaluating the
collaboration gain of even a single coalition, involves solving a Vehicle
Routing Problem (VRP) [11, 38] which is NP-Hard. For detailed calculations of
the collaboration gain, see Appendix A.

Our first contribution is modelling the collaborative routing problem as a
_coalitional bargaining game_ [34, 32]. Here, agents attempt to reach
agreement on how to divide a “pie” (e.g. $100) between them through multiple
rounds of bargaining (see Section 3 for a formal definition). A benefit of
this approach is that we consider both the routing problem (Who should deliver
which requests?) and the gain sharing problem (Who receives how much of the
added value?) simultaneously, whereas a key limitation of many previous
methods consider these sub-problems in isolation from one another [16].
Moreover, our approach is agnostic to the underlying routing problem – the
complexity of the VRP formulation could be increased by further constraints
such as time-windows, without further modification to the method.

Our second contribution is that agents do not need access to the
characteristic function explicitly. Instead, they can reason about the
characteristic function through only receiving a high-dimensional graph input
of delivery information (for example, latitudes and longitudes), as well as
other agents’ actions. This eliminates the need to fully evaluate the
characteristic function, which involves solving an NP-Hard problem an
exponential number of times. Instead, we only need to solve this NP-Hard
problem once when deployed in real-world settings, thus allowing our approach
to scale.

To evaluate performance, we measure the agents’ ability to select the optimal
coalitions. The optimal coalition for player $i$ is defined as the coalition
that maximises the value player $i$ receives after sharing gains amongst all
agents in the coalition. Firstly, we measure the accuracy of agents to propose
the optimal coalitions. Next, we also measure the average optimality gap
between the value of their proposed coalition and the value of the optimal
coalitions.

The work of [32] analyses this coalitional bargaining game (but not in a
collaborative routing setting), and proves that agents will cooperate by
sharing gains equally (i.e. the egalitarian split) in our setting. Due to the
early-stage nature of this paper, we first focus on developing strong RL
agents that can identify and join the optimal coalitions. This is already
challenging as it requires agents to reason about the characteristic function
through receiving a high-dimensional input. We simplify the action space by
assuming all agents adopt the egalitarian split as expected by [32]. We leave
it to future work where agents can also learn this aspect themselves.

##  2 Related Work

Prior collaborative routing literature tackles the partner selection sub-
problem (i.e. who should each carrier work with?) by estimating the
collaboration gain between different partners (or companies) [33, 1]. However,
a limitation of this approach is that they do not consider how much each agent
should be compensated, nor if agents even agree to join the same coalitions
i.e. if the coalitions are stable. Therefore, posing this problem as a
_coalitional bargaining_ game aims to rectify this limitation.

The majority of the Collaborative Routing literature is concerned with the
exchange of _individual_ delivery requests amongst the carriers. This can be
divided into three types of planning approaches: centralised; decentralised
without auctions; and decentralised with auctions [16, 17].

Centralised planning approaches desire to simply maximise social welfare.
Typically, this goal is achieved by using a form of mixed integer linear
programming or (meta)heuristics [8, 27, 30]. This can be viewed as a common-
payoff setting where all agents receive the same reward. However, assuming a
common-payoff setting in practice is restrictive as companies are _self-
interested_ – they mostly care only about their own profits [9]. Therefore, a
more realistic setting is that of decentralised control with self-interested
agents.

A key challenge in decentralised settings is managing the explosion in the
number of _bundles_. Consider Figure 1 where Agent 2 may desire to sell
delivery node 10 perhaps to Agent 1. However, if Agent 2 offers both nodes 10
and 11 as a _bundle_ , then Agent 2 may be able to command a higher price.
Indeed, the number of possible bundles scales $\mathcal{O}(2^{m})$ where $m$
is the number of deliveries. To manage this explosion, a heuristic is
typically implemented where agents can only submit or request a few bundles
(sometimes only one) which would severely impact optimality [4, 26, 45].

A second challenge is to also elicit other agents’ preferences over all
bundles. One approach is to invoke structure on the problem in the form of
_combinatorial auctions_ [7] which aids optimality [25, 16, 18]. However, this
additional structure comes at additional computational complexity. Moreover,
in auction mechanism design, there are four desirable properties: efficiency;
individual rationality; incentive compatibility; and budget balance. [18]
shows that the Vickrey-Clarke-Groves Auction and both of their approaches are
unable to satisfy all four properties simultaneously and there exists a trade-
off instead.

Coalition formation has also been extensively studied in cooperative game
theory [5, 44, 36]. However, much of the existing literature assumes that the
characteristic function is given. Alternatively, they aim to find more
succinct representations of the characteristic function, typically at a cost
of increased computational complexity when computing solution concepts [44].
Examples include Induced Subgraph Games and Marginal Contribution Nets [13,
21]; however, even these representation schemes require evaluating the value
of multiple coalitions (and thus solving multiple VRPs). We argue that many
real-world scenarios consist of the characteristic function being a function
of the agents’ assets or capabilities – in our case, a function of the
deliveries an agent possesses. We therefore ask:_“Can agents form optimal
coalitions from the delivery information instead of having access to the
characteristic function?”_. Therefore, our paper can be viewed as using an
alternative, succinct representation scheme and approximating a rational
outcome by using a function approximator.

The most similar work to ours is that of [2, 6]. In Bachrach et al., [2], they
apply Multi-Agent Reinforcement Learning (MARL) to a spatial and non-spatial
Weighted Voting Game, where agents are given access to the characteristic
function (in this case, all agents’ weights and the quota, $q$). In Chang,
[6], they apply MARL to the 2-player multi-issue bargaining game. Here, both
agents’ weights over the multiple issues are held constant and thus agents do
not require the weights as explicit input to the agents. In our work, we
consider the three agents case (which can be extended to the $n$-player
setting), where the joint policy space is richer, potentially allowing for
collusion. However, for simplicity, we focus only on the single-issue setting.

Finally, significant progress has been made in applying single-agent RL to
Travelling Salesman and Vehicle Routing Problems [42, 31, 24, 22, 28, 43]. We
build upon prior work by adopting the neural network design of [24] trained
via Independent Proximal Policy Optimisation (I-PPO) [35, 12] with a few
modifications due to the high computational cost of stepping our environment.
Further details on agent design are given in Section 4.

##  3 Background

#### Coalitional Games

We consider the $n$-player coalitional game (or cooperative game) with a set
of agents $N=\\{1,\dots,n\\}$. A _coalition_ is defined as a subset of N, i.e.
$S\subseteq N$. The set of all coalitions is denoted $\Sigma$. The _grand
coalition_ is where the coalition consists of all agents in N, i.e. $S=N$. A
_singleton coalition_ is where the coalition consists of only one agent, i.e.
$|S|=1$.

A (transferable utility) coalitional game is a pair $G=\langle N,v\rangle$.
The _characteristic function_ $v:\mathbf{2}^{N}\to\mathbb{R}$ represents the
_value_ (or collaboration gain in our setting) that a given coalition $S$
receives. Like [32], we assume that the characteristic function is
_0-normalised_ , _essential_ and _super-additive_. The characteristic function
is _0-normalised_ if the value of all singleton coalitions is 0, i.e.
$v(\\{i\\})=0,\forall i\in N$. It is _essential_ if the value of the grand
coalition is strictly positive, $v(N)>0$. It is _super-additive_ if $v(S\cup
T)\geq v(S)+v(T)$ for all coalition pairs $(S,T)\in\Sigma$ where $S\cap
T=\varnothing$.

The payoff vector $\mathbf{x}^{S}=\langle x_{1},\dots,x_{k}\rangle,\
\sum_{i}{x_{i}}=1,\ x_{i}\in[0,1]$ denotes how the value achieved by a
coalition $S$ is distributed amongst its $k$ members. The payoff vector is
_feasible_ if $\sum_{i\in S}x_{i}^{S}\leq v(S)$. The set of all feasible
payoff vectors for a given coalition S is $X^{S}$, and $X_{+}^{S}$ when all
the elements of $X^{S}$ is non-negative.

#### Coalitional Bargaining

[32] presents the $n$-player, random proposers, alternating offers coalitional
bargaining game which we adopt. For simplicity, we focus on the 3-player
setting. At every timestep $t=1,2,\dots$ an agent from N is selected uniformly
at random to be the _proposer_. The proposer, player $i$, has two actions.
Firstly, to propose a coalition $S$ with $i\in S\subseteq N$ and $v(S)>0$. Due
to the characteristic function being 0-normalised this implies $|S|\geq 2$.
Secondly, the proposer proposes a payoff vector $\mathbf{x}^{S}\in X_{+}^{S}$.
We currently assume the payoff vector to be an egalitarian split – all agents
within coalition $S$ distribute the value equally. Next, the remaining agents
called the _responders_ decide to either accept or reject the proposal. If all
agents in the proposed coalition $S$ accepts, then the episode terminates with
agent $i$ receiving a reward of $v(S)\cdot x_{i}$ if $i$ is in the coalition
$S$, and 0 otherwise (as the game is 0-normalised). If any responder agents in
$S$ rejects the proposal, then the next round of bargaining begins with a new
proposer selected uniformly at random and the timestep incremented by 1. Since
a discount factor of $\gamma$ is applied which decreases the value received,
this encourages agents to reach agreement within the first timestep as shown
in [32]. The discount factor in this setting is analogous to the _patience_ of
an agent, or the urgency of the delivery decision. The episode continues until
either agreement is reached, or the finite time horizon is reached. If the
episode is terminated due to reaching the finite time horizon, all agents
receive a reward of 0.

We note that whilst Okada,’s results assumes the infinite horizon case, like
in Bachrach et al.,, we assume a finite time horizon as Okada, shows that
agreement should be reached without delay, i.e. at time $t=1$. This is
primarily due to the use of the discount factor $\gamma$ which decreases the
overall collaboration gain that could be achieved as the number of rounds
increases. Therefore, agents are incentivised to reach an agreement as soon as
possible.

#### Multi-Agent Reinforcement Learning

We model this problem as an $n$-player _Stochastic Game_ i.e. assuming perfect
information over states, but not other agents’ preferences. This can be
defined as a tuple $\langle N,S,A,\mathcal{T},\mathcal{R}\rangle$ where:

  * •

$N$ denotes the set of $n$ agents

  * •

$S$ denotes the set of states

  * •

$A=A_{i}\times\dots\times A_{n}$ denotes the set of joint actions, where
$A_{i}$ is player $i$’s set of actions.

  * •

$\mathcal{T}:S\times A\to S$ denotes the transition probabilities

  * •

$\mathcal{R}:S\times A\times S\to\mathbb{R}$ denotes the reward function

For every timestep $t$, an agent $i$ observes the global state $s$ and outputs
an action $a_{i,t}$ sampled from its policy $\pi_{\theta_{i}}(a_{i,t}\mid
s_{t})$ which is parameterised by $\theta_{i}$. This results in a reward
$r_{i,t}=\mathcal{R}(s_{t},\mathbf{a},s_{t+1})$. The return $G_{i}$ is
discounted by a factor $\gamma$, given by
$G_{i}=\sum_{t=1}^{T}\gamma^{t-1}r_{i,t}$. Agent $i$’s objective is to
maximise its expected discounted sum of rewards
$\mathbb{E}[\sum_{t=1}^{T}\gamma^{t-1}r_{i,t}]$.

##  4 Agent Design

#### Input Space

The agents receive a variety of input from the state:

A _location_ can refer to either a depot location or a customer location. A
location is defined as a 4-tuple $\langle
x,y,\texttt{owner},\texttt{is\\_depot}\rangle$. Indeed, $x$ and $y$ refers to
the $x,y$ coordinates of the location. The owner is an integer denoting to
which agent the location belongs to. is_depot takes the value 1 if the
location is a depot, and 0 if the location is a customer.

Note that we consider perfect information, that is, each agent sees
information over all depots’ and customers’ locations (including those that
are owned by other agents). However, agents do not have access to other
agents’ preferences.

The agents also receive the current timestep $t$, or bargaining round, as
input. In addition, the actions taken by other agents are received as well,
denoted by actions_taken. This has a constant shape of
$\langle\texttt{max\\_time\\_horizon},\ 2\times\texttt{n\\_agents}\rangle$,
where max_time_horizon is the maximum number of timesteps before the episode
is forcefully terminated, and n_agents the number of agents. The actions_taken
matrix is initialised all as $-1$. The first three columns are populated with
the proposed pay-off vectors at each round of bargaining. The last three
columns are populated with the responses of each agent at each round of
bargaining.

#### Actor Feature Extractor Design

The feature extractor is based on a modified version of [24] as their model
achieves high performance across a range of routing problems. Similar to [24],
we use an embedding dimension of 128.

The location information is fed through three encoder layers of a Transformer
[41] to form our node embeddings. However, [24] finds that replacing layer
normalisation with batch normalisation improves performance and thus we adopt
batch normalisation as well. To obtain our graph embedding, we take a mean
over the node embeddings. We additionally embed the timestep $t$ into a
128-dimensional vector to provide context, which we sum with the node
embeddings. We then take a single-head attention over the node embeddings to
form a resultant embedding. The resultant embedding is then used to decide
which actions to take.

#### Action Space

The agents have three action heads: _coalitions_ , _proposals_ and _response_.

coalitions takes the form $\\{0,1\\}^{|N|}$ where $|N|$ is the total number of
agents, in this case, 3. This action denotes whether the respective agent is
part of the coalition $S$. Note that this game assumes that player $i$ is in
the coalition $S$, i.e. $S_{i}=1$. The coalitions action head takes the
resultant embedding followed by four dense layers with 256 hidden neurons and
a ReLU activation function. The output is then passed through $n$ independent
Bernoulli distributions to determine the probability that a given agent exists
within the coalition $S$.

proposals is a vector $\mathbf{x}\in\mathbb{R}_{+}^{|N|}$ where
$\sum_{i}{x_{i}}=1,\ x_{i}\in[0,1]$. This vector denotes how much of the
collaboration gain is assigned to each respective agent (as a percentage).
Note that this action head is not strictly necessary in this current paper as
we assume an egalitarian split; however, we would require this in future work
if agents are to learn the gain sharing mechanism themselves. Furthermore,
note that this is a continuous action space, as opposed to the other actions
which are discrete. To calculate proposals, the resultant embedding is passed
through a single dense layer with 64 hidden neurons and a $\tanh$ activation
function. From this, a mean and standard deviation is output and fed into a
multi-variate Guassian. This Guassian is then sampled from to form the logits
of the pay-off vector $\mathbf{x}$. This pay-off vector is then masked by the
_coalitions_ vector, i.e. if a player $i$ is not in the coalition $S$, it will
also receive 0 in the pay-off vector. An egalitarian split is forced by
masking the logits to be 1 if the respective agent is in the proposed
coalition $S$, and an arbitarily large negative number otherwise. Finally, a
softmax is taken over the logits to ensure the constraint $\sum_{i}{x_{i}}=1$.

responses is a binary action that denotes whether an agent accepts or rejects
a given proposal. It takes the resultant embedding followed by four dense
layers with 256 hidden neurons and a ReLU activation function. The output is
then fed through a Bernoulli distribution.

#### Loss Function

The agents are trained through independent learning [37] in the form of
Independent PPO [35, 12]. However, modifications had to be made. Firstly, we
use a baseline instead of a critic (i.e. no bootstrapping in value function
estimates) as the design of a good critic for vehicle routing problems is
“non-trivial” [24]. To perform the greedy roll-out proposed by [24] is
expensive due to the computational cost of stepping our environment. Thus, we
use the baseline in [31] which estimates the value only of the initial state,
$S_{0}$ (see below paragraph on baseline design 4); this has an appealing
interpretation as the _difficulty_ of a given problem instance. For
simplicity, we do not use shared parameters for the actor and baseline to
avoid tuning an additional hyperparameter, albeit potentially at the cost of
sample efficiency. Note that we also normalise the advantage function due to
the small magnitude in rewards.

The loss function is the clipped objective function of [35] with entropy
regularisation:

$\mathcal{L}_{t}^{\text{CLIP}}(\theta)=\hat{\mathbb{E}_{t}}[\min(r_{t}(\theta)\hat{A}_{t},\,\text{clip}(r_{t}(\theta),1-\epsilon,1+\epsilon)\,\hat{A}_{t})+\beta
S[\pi_{\theta}](s_{t})]$

where $\beta$ is the entropy regularisation coefficient and $S$ the entropy
bonus. Strong entropy regularisation was employed. The entropy regularisation
coefficient $\beta$ was linearly annealed from 0.75 to 0 over 10,000 epochs,
followed by $\max(\beta,0.2)$. Such high entropy regularisation was essential
in this setting, as a strong locally optimal policy would be to always propose
the grand coalition and always accept every proposal. Without entropy
regularisation, the agents would converge on this local optimum too quickly
and agents would not learn to extract useful features over the states.

Due to computational limitations, a hyperparameter search was not performed.
We suspect that small improvements to sample efficiency could be made by
finding more optimal hyperparameter settings and schedules.

#### (Reward) Baseline Design

A useful baseline helps reduce the variance in policy-gradient methods. The
baseline we use estimates the value of the initial state
$V_{\pi_{\theta}}(S_{0})$. This has an appealing interpretation as being
approximately the difficulty of a given problem instance [31]. The baseline
adopts the same neural network architecture as the actor to arrive at the
resultant embedding. However, we then feed this embedding through six dense
layers of 256 hidden neurons using a ReLU activation function. The baseline is
trained using a clipped mean squared error loss as in [35].

##  5 Experiments

Figure 2: A plot of the distribution of depot and customer locations. Depots
are fixed denoted by squares. Customers may be uniformly at random located
within any of agent’s respective circles.

#### Problem Setting

We base our problem setting on a modified version of [15]. We consider an
environment with three companies, each represented by an agent. Each agent has
one depot and three customers that it must deliver to. The depot $(x,y)$
locations for each agent are held fixed at
$\\{(-0.2,0.173),(0.2,0.173),(0,-0.173)\\}$ respectively. The customer $(x,y)$
locations are generated uniformly at random within a certain radius. The
radius for each instance is selected uniformly at random from the set
$\\{0.3,0.4,0.6\\}$. Varying the radius has the effect of varying the degree
of overlap between agents – and thus opportunity for collaboration. This can
be seen in Figure 2. Future work will consider an increased number of agents
and real-world delivery data.

To calculate the pre-collaboration and post-collaboration gains, the shortest
paths are calculated exactly using Gurobi [20]. The pre-collaboration shortest
paths can be calculated by solving three (un)Capacitated Vehicle Routing
Problems (one for each agent). The post-collaboration shortest paths are
calculated by solving a single Multi-Depot Vehicle Routing Problem. Problem
formulations for the Capacitated VRP and Multi-Depot VRP can be found in
Appendix B and Appendix C respectively. Capacity is effectively removed by
setting the capacity of each vehicle to an arbitrarily large number and the
weight of each delivery to 1.

#### Experimental Design

We perform 5 independent runs with different random seeds to train our agents.
Agents are trained for 10,000 epochs and evaluated every 100 epochs. Every
training instance is randomly generated on the fly – thus agents never see the
same instance twice, both during training and evaluation. We train using a
batch size of 256 and evaluate with a batch size of 512. All agents use a
discount factor $\gamma$ of 0.99. The maximum number of bargaining rounds $T$
is set to 10. The learning rate was held constant at $3\text{\times}{10}^{-4}$
and we use Adam [23]. All code is run on a single desktop with an Intel® Core™
i7-9700K processor with an NVIDIA® GeForce® RTX 2080Ti GPU. We use TensorFlow
2.4.1 [29] to implement our agents.

#### Heuristic Bot

To compare our algorithm with a simple baseline algorithm, we hand-craft a bot
that always forms the grand coalition and always accepts every proposal. This
bot has two interpretations. Firstly, due to the slight class imbalance, this
is simply a bot that always selects the majority class for each agent in the
proposed coalition. An alternative interpretation is that this bot can be
viewed as an algorithm with a common-payoff. This has the effect of maximising
the collaboration gain (as this game is super-additive), but may not always
provide stable coalitions (it may be rational for some agents to defect from a
given coalition). This can be most easily illustrated in settings where a
_dummy player_ exists. A dummy player is a player who adds no value to a given
coalition. An obvious example would be a player $i$ who is situated far away
from the other two players $j$ and $k$ in the collaborative routing setting.

###  5.1 Evaluation Metrics

#### Accuracy

A simple evaluation metric is to measure how often the agents propose the
correct coalition. For player $i$, the correct coalition $S_{i}^{*}$ is
defined to be the coalition $S$ which would maximise player $i$’s reward. This
involves brute forcing the characteristic function to evaluate the value of
each possible coalition which is only possible since we consider 3 agents. The
reward $R$ is the collaboration gain from agreeing to coalition $S$, $v(S)$,
multiplied by the $i$th element of the pay-off vector, $x_{i}$.

#### Optimality Gap

We denote the absolute and relative optimality gap of player $i$ by $\phi_{i}$
and $\eta_{i}$ respectively. The absolute optimality gap $\phi_{i}$ for player
$i$ is defined as $\phi_{i}=v(S_{i}^{*})-v(S)$, where $S_{i}^{*}$ is the
correct coalition, $S_{i}$ is player $i$’s proposed coalition, and $v(\cdot)$
is the characteristic function (i.e. the collaboration gain of a given
coalition). The relative optimality gap $\eta_{i}$, is calculated as
$\eta_{i}=\frac{v(S_{i}^{*})-v(S_{i})}{v(S_{i}^{*})}$.

Since the data is randomly generated, there could be scenarios where there is
no value in collaborating, i.e. even the value of the grand coalition is 0,
$v(N)=0$. Note that we exclude these scenarios when calculating the above
evaluation metrics; however, this only occurs 1.9% of the time when brute-
forcing 51,200 instances.

###  5.2 Results

We perform 5 independent runs comparing our RL bot to the heuristic bot, as
well as a random agent which simply proposes coalitions uniformly at random.
We evaluate the performance of our agents both in terms of accuracy,
optimality gap and run-time. This can be seen in Figure 3.

(a) Average Accuracy (%)

(b) Average Optimality Gap (%)

Figure 3: Learning curve taking an (a) average accuracy (b) average optimality
gap across all 3 agents for readability. For individual learning curves of
each agent, see Appendix D and Appendix E. Solid lines denote mean accuracy
across all 15 runs (five for each of the three agents). Shaded regions denote
$\pm$ two standard deviations. After training for 10,000 epochs, our RL agents
reach an average accuracy of 79%, outperforming the heuristic bot.

From Figures 3a and 3b, we conclude that our agents have learnt to outperform
the heuristic bot reaching an average accuracy of 79% and average optimality
gap of 0.01 (or 4.2%). For individual learning curves of absolute optimality
gap, see Appendix F. Furthermore, our RL agents are able to reach agreement in
512 parallel instances within an average of 9.2s (or 0.018s per instance). We
note that to brute force the characteristic function for 512 instances takes
24.3s (or 0.047s per instance). Thus, our RL agents achieve a 62% reduction in
computational time. Whilst 0.047s per instance may seem reasonable even when
applying brute force, we stress that this is due to the simplistic VRP setting
we consider – this will not scale with the number of agents nor problem
complexity via additional constraints such as time-windows. Moreover, our
agents reach agreement in a distributed and self-interested manner, which
overcomes the limitations of central orchestration methods mentioned in
Section 2. We hypothesise that a key remaining challenge is due to the low
magnitude in absolute optimality gap of 0.01. To overcome this, future work
could investigate this challenge from either an algorithmic improvement point
of view, or simply use a more realistic revenue and cost structure for the
underlying routing problem (where the magnitude is much greater). In addition,
there exists symmetries in our delivery information. Future work could
investigate learning more robust representations over the input, such that the
actions taken are invariant to linear transformations, such as rotation and
reflection.

It is also difficult to parallelise this workflow efficiently when calculating
the collaboration gains as it is all performed on CPU. It should be possible
to sacrifice the mathematical guarantee of optimality with approximate
solution methods, such as (meta)heuristics, or indeed – reinforcement
learning. Using RL in the environment would allow us to easily shift the
computation to the GPU and should enable a larger number of instances to be
solved in parallel. However, it would be important to ensure the robustness of
these RL models first and thus we leave it for future work.

##  6 Conclusion

In this paper, we propose to tackle the challenging real-world problem of
Collaborative Vehicle Routing through a coalitional bargaining lens with
reinforcement learning. The main challenge of our setting is the inability of
extant methods to fully evaluate the characteristic function due to high
computational complexity. The RL agents designed in this work are able to
correctly reason over a high-dimensional graph input to _implicitly_ reason
about the characteristic function instead. This eliminates the need to
evaluate the Vehicle Routing Problem (VRP) an exponential number of times and
increase its practicability as we only need to perform this once. Another
important point is that collaboration is not centrally orchestrated but
facilitated using distributed decision making. This marks an important step
towards real-world adoption which would aid transportation planners to
consider more possible collaboration scenarios.

The current focus of this work is to obtain strong RL agents that can identify
and join the optimal coalitions. Whilst we have achieved this, a key current
assumption in our work is the egalitarian split (where agents split the
collaboration gain equally). Instead, in future work agents should learn this
themselves. Furthermore, future work should investigate more realistic routing
scenarios. In particular, capacity is currently represented as a scalar value,
whereas real-world settings may require considering the 3-dimensional volume
of packages as well. This would lead to a “Three-dimensional Loading
Capacitated Vehicle Routing Problem” (3L-CVRP).

##  7 Broader Impact

We believe that tackling collaborative routing through coalitional bargaining
and MARL presents a promising research direction which allows reduction in
cost, greenhouse gas emissions and road congestion. More broadly, this could
be applied to _collaborative logistics_ where, for example, companies
cooperate by sharing warehouse space and where again gain sharing is
essential.

However, there still remains a research gap before this approach can be
deployed in the real-world. Firstly, can this system be gamed? Can companies
(un)intentionally mis-report the deliveries they possess in the form of
_phantom_ or _decoy_ delivery tasks? More generally, how can we verify
companies’ information in a sensitive manner such that all parties trust the
system? Whilst these are interesting and challenging research questions, we
believe that multi-agent reinforcement learning could be a useful tool (in
addition to other approaches, such as game theory) to help address these
questions.

## Acknowledgments and Disclosure of Funding

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) grant on “Intelligent Systems for Supply Chain Automation”
under Grant Number 2275316, as well as by the UK EPSRC Connected Everything
Network Plus under Grant EP/S036113/1.

We thank Eugene Vinitsky for his mentorship which improved the quality of this
paper and the Cooperative AI Organising Committee for facilitating the
mentorship scheme. We also thank the Manufacturing Analytics Group for their
insightful discussions regarding early drafts of this paper.

## References

  * Adenso-Díaz et al., [2014] Adenso-Díaz, B., Lozano, S., and Moreno, P. (2014).  Analysis of the synergies of merging multi-company transportation needs.  Transportmetrica A: Transport Science, 10(6):533–547.  Publisher: Informa UK Limited.
  * Bachrach et al., [2020] Bachrach, Y., Everett, R., Hughes, E., Lazaridou, A., Leibo, J. Z., Lanctot, M., Johanson, M., Czarnecki, W. M., and Graepel, T. (2020).  Negotiating Team Formation Using Deep Reinforcement Learning.  arXiv:2010.10380 [cs].  arXiv: 2010.10380.
  * Ballot and Fontane, [2010] Ballot, E. and Fontane, F. (2010).  Reducing transportation CO2 emissions through pooling of supply networks: perspectives from a case study in French retail chains.  Production Planning & Control, 21(6):640–650.  Publisher: Informa UK Limited.
  * Bo Dai and Chen, [2009] Bo Dai and Chen, H. (2009).  Mathematical model and solution approach for collaborative logistics in less than truckload (LTL) transportation.  In 2009 International Conference on Computers & Industrial Engineering. IEEE.
  * Chalkiadakis et al., [2011] Chalkiadakis, G., Elkind, E., and Wooldridge, M. (2011).  Computational Aspects of Cooperative Game Theory (Synthesis Lectures on Artificial Inetlligence and Machine Learning).  Morgan & Claypool Publishers, 1st edition.
  * Chang, [2020] Chang, H.-C. H. (2020).  Multi-Issue Bargaining With Deep Reinforcement Learning.  arXiv:2002.07788 [cs].  arXiv: 2002.07788.
  * Cramton et al., [2006] Cramton, P., Shoham, Y., and Steinberg, R. (2006).  Combinatorial Auctions.  The MIT Press.
  * [8] Cruijssen, F., Bräysy, O., Dullaert, W., Fleuren, H., and Salomon, M. (2007a).  Joint route planning under varying market conditions.  International Journal of Physical Distribution & Logistics Management, 37(4):287–304.  Publisher: Emerald.
  * [9] Cruijssen, F., Cools, M., and Dullaert, W. (2007b).  Horizontal cooperation in logistics: Opportunities and impediments.  Transportation Research Part E: Logistics and Transportation Review, 43(2):129–142.  Publisher: Elsevier BV.
  * [10] Cruijssen, F., Dullaert, W., and Fleuren, H. (2007c).  Horizontal cooperation in transport and logistics: A literature review.  Transportation Journal, 46.
  * Dantzig and Ramser, [1959] Dantzig, G. B. and Ramser, J. H. (1959).  The Truck Dispatching Problem.  Management Science, 6(1):80–91.
  * de Witt et al., [2020] de Witt, C. S., Gupta, T., Makoviichuk, D., Makoviychuk, V., Torr, P. H. S., Sun, M., and Whiteson, S. (2020).  Is Independent Learning All You Need in the StarCraft Multi-Agent Challenge?  arXiv:2011.09533 [cs].  arXiv: 2011.09533.
  * Deng and Papadimitriou, [1994] Deng, X. and Papadimitriou, C. H. (1994).  On the Complexity of Cooperative Solution Concepts.  Mathematics of Operations Research, 19(2):257–266.  Publisher: INFORMS.
  * Ferrell et al., [2020] Ferrell, W., Ellis, K., Kaminsky, P., and Rainwater, C. (2020).  Horizontal collaboration: opportunities for improved logistics planning.  International Journal of Production Research, 58(14):4267–4284.  Publisher: Informa UK Limited.
  * [15] Gansterer, M. and Hartl, R. F. (2018a).  Centralized bundle generation in auction-based collaborative transportation.  OR Spectrum, 40(3):613–635.  Publisher: Springer Science and Business Media LLC.
  * [16] Gansterer, M. and Hartl, R. F. (2018b).  Collaborative vehicle routing: A survey.  European Journal of Operational Research, 268(1):1–12.  Publisher: Elsevier BV.
  * Gansterer and Hartl, [2020] Gansterer, M. and Hartl, R. F. (2020).  Shared resources in collaborative vehicle routing.  TOP, 28(1):1–20.  Publisher: Springer Science and Business Media LLC.
  * Gansterer et al., [2019] Gansterer, M., Hartl, R. F., and Vetschera, R. (2019).  The cost of incentive compatibility in auction-based mechanisms for carrier collaboration.  Networks, 73(4):490–514.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/net.21828.
  * Guajardo and Rönnqvist, [2016] Guajardo, M. and Rönnqvist, M. (2016).  A review on cost allocation methods in collaborative transportation.  International Transactions in Operational Research, 23(3):371–392.  Publisher: Wiley.
  * Gurobi Optimization, LLC, [2021] Gurobi Optimization, LLC (2021).  Gurobi Optimizer Reference Manual.
  * Ieong and Shoham, [2005] Ieong, S. and Shoham, Y. (2005).  Marginal contribution nets: a compact representation scheme for coalitional games.  In Proceedings of the 6th ACM conference on Electronic commerce - EC ’05, pages 193–202, Vancouver, BC, Canada. ACM Press.
  * Joshi et al., [2019] Joshi, C. K., Laurent, T., and Bresson, X. (2019).  An Efficient Graph Convolutional Network Technique for the Travelling Salesman Problem.  arXiv:1906.01227 [cs, stat].  arXiv: 1906.01227.
  * Kingma and Ba, [2017] Kingma, D. P. and Ba, J. (2017).  Adam: A Method for Stochastic Optimization.  arXiv:1412.6980 [cs].  arXiv: 1412.6980.
  * Kool et al., [2019] Kool, W., van Hoof, H., and Welling, M. (2019).  Attention, Learn to Solve Routing Problems!  arXiv:1803.08475 [cs, stat].  arXiv: 1803.08475.
  * Krajewska et al., [2008] Krajewska, M. A., Kopfer, H., Laporte, G., Ropke, S., and Zaccour, G. (2008).  Horizontal cooperation among freight carriers: request allocation and profit sharing.  Journal of the Operational Research Society, 59(11):1483–1491.  Publisher: Informa UK Limited.
  * Li et al., [2015] Li, J., Rong, G., and Feng, Y. (2015).  Request selection and exchange approach for carrier collaboration based on auction of a single request.  Transportation Research Part E: Logistics and Transportation Review, 84:23–39.  Publisher: Elsevier BV.
  * Lin, [2008] Lin, C. (2008).  A cooperative strategy for a vehicle routing problem with pickup and delivery time windows.  Computers & Industrial Engineering, 55(4):766–782.  Publisher: Elsevier BV.
  * Lu et al., [2020] Lu, H., Zhang, X., and Yang, S. (2020).  A Learning-based Iterative Method for Solving Vehicle Routing Problems.  In International Conference on Learning Representations.
  * Martín Abadi et al., [2015] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Jia, Y., Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng (2015).  TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems.
  * Montoya-Torres et al., [2016] Montoya-Torres, J. R., Muñoz-Villamizar, A., and Vega-Mejía, C. A. (2016).  On the impact of collaborative strategies for goods delivery in city logistics.  Production Planning & Control, 27(6):443–455.  Publisher: Informa UK Limited.
  * Nazari et al., [2018] Nazari, M., Oroojlooy, A., Snyder, L. V., and Takáč, M. (2018).  Reinforcement Learning for Solving the Vehicle Routing Problem.  arXiv:1802.04240 [cs, stat].  arXiv: 1802.04240.
  * Okada, [1996] Okada, A. (1996).  A Noncooperative Coalitional Bargaining Game with Random Proposers.  Games and Economic Behavior, 16(1):97–108.  Publisher: Elsevier BV.
  * Palhazi Cuervo et al., [2016] Palhazi Cuervo, D., Vanovermeire, C., and Sörensen, K. (2016).  Determining collaborative profits in coalitions formed by two partners with varying characteristics.  Transportation Research Part C: Emerging Technologies, 70:171–184.  Publisher: Elsevier BV.
  * Rubinstein, [1982] Rubinstein, A. (1982).  Perfect Equilibrium in a Bargaining Model.  Econometrica, 50(1):97.  Publisher: JSTOR.
  * Schulman et al., [2017] Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O. (2017).  Proximal Policy Optimization Algorithms.  arXiv:1707.06347 [cs].  arXiv: 1707.06347.
  * Shoham and Leyton-Brown, [2008] Shoham, Y. and Leyton-Brown, K. (2008).  Multiagent Systems: Algorithmic, Game-Theoretic, and Logical Foundations.  Cambridge University Press, USA.
  * Tan, [1993] Tan, M. (1993).  Multi-Agent Reinforcement Learning: Independent versus Cooperative Agents.  In Proceedings of the Tenth International Conference on International Conference on Machine Learning, ICML’93, pages 330–337, San Francisco, CA, USA. Morgan Kaufmann Publishers Inc.  event-place: Amherst, MA, USA.
  * Toth and Vigo, [2014] Toth, P. and Vigo, D., editors (2014).  Vehicle Routing: Problems, Methods, and Applications, Second Edition.  Number 18 in MOS-SIAM Series on Optimization. SIAM.
  * UK Department for Business, Energy & Industrial Strategy, [2021] UK Department for Business, Energy & Industrial Strategy (2021).  Final UK greenhouse gas emissions national statistics.
  * UK Department for Transport, [2020] UK Department for Transport (2020).  Road freight statistics: 2019.
  * Vaswani et al., [2017] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. (2017).  Attention Is All You Need.  arXiv:1706.03762 [cs].  arXiv: 1706.03762.
  * Vinyals et al., [2017] Vinyals, O., Fortunato, M., and Jaitly, N. (2017).  Pointer Networks.  arXiv:1506.03134 [cs, stat].  arXiv: 1506.03134.
  * Wang et al., [2021] Wang, H., Zong, Z., Xia, T., Luo, S., Zheng, M., Jin, D., and Li, Y. (2021).  Rewriting by Generating: Learn Heuristics for Large-scale Vehicle Routing Problems.
  * Wooldridge, [2009] Wooldridge, M. (2009).  An Introduction to MultiAgent Systems.  Wiley Publishing, 2nd edition.
  * Xu et al., [2017] Xu, S. X., Huang, G. Q., and Cheng, M. (2017).  Truthful, Budget-Balanced Bundle Double Auctions for Carrier Collaboration.  Transportation Science, 51(4):1365–1386.  Publisher: INFORMS.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes] See Abstract and Section
1

    3. (c)

Did you discuss any potential negative societal impacts of your work? [Yes]
See Section 7

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [N/A]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [No]
We will publicly release code after the extension of this work to a journal
paper.

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes] See Section 5

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes] See Figure 3, Figure 5 and Figure 6

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] All code is
run on a single desktop with an Intel® Core™ i7-9700K processor with an
NVIDIA® GeForce® RTX 2080Ti GPU. A single independent run takes 31.5 hours to
train, but we note that after training, it only takes 9 seconds for all agents
to reach agreement.

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes] We cite
the use of the neural network design of [24] and the benchmark design of [15].
We also thank TensorFlow [29] and Gurobi [20].

    2. (b)

Did you mention the license of the assets? [N/A]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [No]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A] We do not work with human subjects at this
early stage of work.

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Appendix A - Detailed Calculation of an Optimal Coalition,
Collaboration Gain and the Characteristic Function

(a) Pre-collaboration
Total Cost = 3.35

(b) Post-collaboration (with grand coalition, $\\{1,2,3\\}$)
Total Cost = 2.47

(c) Post-collaboration (with a coalition $\\{1,2\\}$ i.e. Agent 3 is excluded
from the coalition)
Total Cost = 2.59

Figure 4: Three agents, Agents 1, 2 and 3 are denoted by the colours green,
orange and purple respectively. Squares denote depots. Crosses denote customer
locations. Node indices (arbitrary) are denoted in black, with costs given in
their respective colors. The _collaboration gain_ is defined as the difference
in social welfare before and after collaboration. Figures 4b and 4c refer to
two possible post-collaboration scenarios with collaboration gains per capita
of 0.29 and 0.38 respectively. Thus, it would be rational for the coalition
$\\{1,2\\}$ to form instead of the grand coalition $\\{1,2,3\\}$.

We denote the set of $n$ agents as $N=\\{1,\dots,n\\}$. A coalition is a
subset of $N$, i.e. $S\subseteq N$. The grand coalition is where all agents
are in the coalition, i.e. $S=N$.

Pre-collaboration Profit and Social Welfare: The _pre-collaboration profit_ of
Agent 1 is calculated as follows: the _Revenue_ is 3 (1 for each delivery);
the _Cost_ is 1.42 (sum of the edge distances); thus the _Profit_ is 1.58
(Revenue subtract Cost). Similarly, the pre-collaboration profit of Agents 2
and 3 is 2 and 2.07. The _pre-collaboration social welfare_ is the sum of the
pre-collaboration profits, thus $1.58+2+2.07=5.65$.

Post-collaboration Profit and Social Welfare: Similarly, the post-
collaboration profit of Agent 1 is $1-(0.06+0.06)=0.88$. For Agents 2 and 3 it
is 2.19 and 3.46 respectively. Thus a _post-collaboration social welfare_ of
$0.88+2.19+3.46=6.53$.

Collaboration Gain: The _collaboration gain_ is defined as the difference in
social welfare before and after collaboration, in this case $6.53-5.65=0.88$.
The _value per capita_ is $\frac{0.88}{3}=0.29$. Note that if only Agents 1
and 2 form a coalition (and exclude Agent 3), then the collaboration gain is
divided by 2 instead – thus making it rational to object and form the
coalition $\\{1,2\\}$ (the value per capita of this coalition is 0.38). Due to
our assumption of the egalitarian split, i.e. agents share all gains equally,
then the reward the agents receive is equal to the value per capita.

Characteristic Function: The characteristic function,
$v:\mathbf{2}^{N}\to\mathbb{R}$ calculates for every possible coalition the
collaboration gain. Importantly, to fully evaluate the characteristic function
would require solving a variant of the Vehicle Routing Problem for every
possible coalition which scales $\mathcal{O}(2^{n})$.

Following the example in Figure 4:

$v(\\{1,2,3\\})=0.88$ | Value per Capita = $\frac{0.88}{3}=0.29$
---|---
$v(\\{1,2\\})=0.76$ | Value per Capita = $\frac{0.76}{2}=0.38$
$v(\\{1,3\\})=0.24$ | Value per Capita = $\frac{0.24}{2}=0.12$
$v(\\{2,3\\})=0.01$ | Value per Capita = $\frac{0.01}{2}=0.005$

Whilst the grand coalition $\\{1,2,3\\}$ maximises the collaboration gain and
thereby social welfare, it is rational for Agents 1 and 2 to object and form
the coalition $\\{1,2\\}$ instead, as this maximises their own value per
capita (and thus reward). Thus, we say that the coalition $\\{1,2\\}$ is
optimal.

Therefore, if Agents 1 or 2 is proposing a coalition, we measure the accuracy
of them correctly proposing the coalition $\\{1,2\\}$. For Agent 3, since it
is required to propose a coalition that contains itself, yet the optimal
coalition does not contain Agent 3, Agents 1 and 2 should object to any
proposal from Agent 3. Thus, we exclude these scenarios when calculating the
accuracy for Agent 3.

##  Appendix B Appendix B - Capacitated Vehicle Routing Problem

In our paper, the pre-collaboration social welfare can be calculated by first
solving three independent Capacitated Vehicle Routing Problems, where we
assume an arbitrarily high capacity for each vehicle.

The Capacitated Vehicle Routing Problem (CVRP) was first proposed by [11].
Since then, many different variants and formulations of the Vehicle Routing
Problem has been proposed [38]. Here we show the _three-index (vehicle-flow)
formulation_.

The CVRP considers the setting where goods are distributed to $n$ customers.
The goods are initially located at the _depot_ , denoted by nodes (or
vertices) $o$ and $d$. Node $o$ refers to the starting point of a route, and
node $d$ the end point of a route. The customers are denoted by the set of
nodes $N=\\{1,2,\dots,n\\}$. Each customer $i\in N$ has a _demand_ $q_{i}\geq
0$. In our setting, we consider $q_{i}=1$ for all customers. A _fleet_ of
$|K|$ vehicles $K=\\{1,2,\dots,|K|\\}$ are said to be _homogeneous_ if they
all have the same capacity $Q>0$. In our setting, we consider only one vehicle
and set its capacity $Q$ to an arbitrarily high number to remove the capacity
constraint. A vehicle must start at the depot, and can deliver to a set of
customers $S\subseteq N$ before returning to the depot. The _travel cost_
$c_{i,j}$ is associated for a vehicle travelling between nodes $i$ and $j$
which we assume to be the Euclidean distance.

This problem can be modelled as a complete directed graph $G=(V,A)$, where the
vertex set $V\coloneqq N\cup\\{o,d\\}$ and the arc set
$A\coloneqq(V\setminus\\{d\\})\times(V\setminus\\{o\\})$. We define the _in-
arcs_ of $S$ as $\delta^{-}(S)=\\{(i,j)\in A:i\notin S,j\in S\\}$. The _out-
arcs_ of $S$ is $\delta^{+}(S)=\\{(i,j)\in A:i\in S,j\notin S\\}$.

The binary decision variables $x_{ijk}$ denotes whether a vehicle $k\in K$
travels over the arc $(i,j)\in A$. The binary decision variables $y_{ik}$
denotes whether a vehicle $k\in K$ visits node $i\in V$. $u_{ik}$ denotes the
load in vehicle $k$ before visiting node $i$. We define the demand at the
depot nodes $o$ and $d$ to be 0, i.e. $q_{o}=q_{d}=0$. This yields:

∑_k ∈K c^T x_k ∑_k ∈K y_ik =1, ∀i ∈N x_k(δ^+(i)) - x_k(δ^-(i)) ={1, i = o,0, i
∈N, ∀i ∈V ∖{d}, k ∈K y_ik= x_k(δ^+(i)) ∀i ∈V ∖{d}, k ∈K y_dk= x_k(δ^-(d)) ∀k
∈K u_ik - u_jk + Qx_ijk≤Q - q_j ∀(i, j) ∈A, k ∈K q_i≤u_ik ≤Q ∀i ∈V, k ∈K
x=(x_k) ∈{0, 1}^K ×A y=(y_k) ∈{0, 1}^K ×V .

  * •

The objective function (B) minimises the Euclidean distance travelled by the
vehicle.

  * •

Constraint (B) ensures the vehicle only visits each customer once.

  * •

Constraint (B) ensures that the sum of vehicles entering node $d$ and exiting
node $d$ is $-1$. This ensures that a vehicle $k$ performs a route starting at
$o$ and ending at $d$.

  * •

Constraint (B and B) couples variables $x_{ijk}$ and $y_{ik}$.

  * •

Constraint (B) is the Miller-Tucker-Zemlin constraint which helps eliminate
subtours.

  * •

Constraint (B) is the capacity constraint.

##  Appendix C Appendix C - Multi-Depot Vehicle Routing Problem

In our paper, the post-collaboration social welfare can be calculated by
solving the Multi-Depot Vehicle Routing Problem (MDVRP) once. The number of
depots corresponds to the number of agents within the accepted coalition.
Again, we remove capacity constraints by setting the capacity of each vehicle
to an arbitrarily large number. However, we add the additional constraint that
each vehicle has to visit at least one customer.

The MDVRP is a simple extension of the CVRP formulation provided in Appendix
B. Instead of having the depot simply represented by nodes $o$ and $d$, the
depots are extended to belong to a specific vehicle $k$ through nodes $o_{k}$
and $d_{k}$. Doing so yields:

∑_k ∈K c^T x_k ∑_k ∈K y_ik =1, ∀i ∈V x_k(δ^+(i)) - x_k(δ^-(i)) ={1, i = ok,0,
i ∈N, ∀i ∈V ∖{d_k}, k ∈K y_ik= x_k(δ^+(i)) ∀i ∈V ∖{d_k}, k ∈K y_d_kk=
x_k(δ^-(d_k)) ∀k ∈K y_d_kk= 1 ∀k ∈K u_ik - u_jk + Qx_ijk≤Q - q_j ∀(i, j) ∈A, k
∈K q_i≤u_ik ≤Q ∀i ∈V, k ∈K x=(x_k) ∈{0, 1}^K ×A y=(y_k) ∈{0, 1}^K ×V .

  * •

The objective function (C) minimises the Euclidean distance travelled by all
vehicles.

  * •

Constraint (C) ensures that each vehicle only visits each customer once.

  * •

Constraint (C) ensures that the sum of vehicles entering node $d_{k}$ and
exiting node $d_{k}$ is $-1$. This ensures that a vehicle $k$ performs a route
starting at $o_{k}$ and ending at $d_{k}$.

  * •

Constraint (C and C) couples variables $x_{ijk}$ and $y_{ik}$.

  * •

Constraint C ensures that each vehicle performs at least one delivery.

  * •

Constraint (C) is the Miller-Tucker-Zemlin constraint which helps eliminate
subtours.

  * •

Constraint (C) is the capacity constraint.

##  Appendix D Appendix D - Individual Learning Curves - Accuracy

(a)  (b)

(c)

Figure 5: Accuracy learning curves for all 3 agents. Solid lines denote mean
accuracy across all 5 respective runs. Shaded regions denote $\pm$ two
standard deviations. After training for 10,000 epochs, Agents 1, 2 and 3 reach
an average accuracy of 79.6%, 80.5% and 77.7% respectively with all RL agents
outperforming the heuristic bot.

##  Appendix E Appendix E - Individual Learning Curves - Relative Optimality
Gap

(a)  (b)

(c)

Figure 6: Relative optimality gap learning curves for all 3 agents. Solid
lines denote mean accuracy across all 5 respective runs. Shaded regions denote
$\pm$ two standard deviations. After training for 10,000 epochs, Agents 1, 2
and 3 reach an average relative optimality gap of 3.6%, 3.2% and 5.9%, with
all RL agents outperforming the heuristic bot.

##  Appendix F Appendix F - Individual Learning Curves - Absolute Optimality
Gap

(a)  (b)

(c)

Figure 7: Absolute optimality gap learning curves for all 3 agents. Solid
lines denote mean accuracy across all 5 respective runs. Shaded regions denote
$\pm$ two standard deviations. After training for 10,000 epochs, Agents 1, 2
and 3 reach an average absolute optimality gap of 0.010, 0.009, 0.014 with all
RL agents outperforming the heuristic bot.


[1]Matthew J. Doyle

1]Department of Physics and Astronomy, University of Manchester, Schuster
Building, Manchester, M13 9PL, United Kingdom

2]Joint Quantum Centre (JQC) Durham–Newcastle, School of Mathematics and
Statistics, Newcastle University, Newcastle upon Tyne, NE1 7RU, United Kingdom

# Modelling turbulent flow of superfluid 4He past a rough solid wall in the
$T=$ 0 limit

<EMAIL_ADDRESS>Andrei I. Golov<EMAIL_ADDRESS>Paul M. Walmsley<EMAIL_ADDRESS>Andrew W. Baggaley
<EMAIL_ADDRESS>[  [

###### Abstract

We present a numerical study, using the vortex filament model, of vortex
tangles in a flow of pure superfluid 4He in the $T=0$ limit through a channel
of width $D=1$ mm for various applied velocities $V$. The flat channel walls
are assumed to be microscopically rough such that vortices terminating at the
walls are permanently pinned; vortices are liberated from their pinned ends
exclusively through self-reconnection with their images. Sustained tangles
were observed, for a period of 80 s, above the critical velocity $V_{c}\sim
0.20$ cm s${}^{-1}=20\kappa/D$. The coarse-grained velocity profile was akin
to a classical parabolic profile of the laminar Poiseuille flow, albeit with a
non-zero slip velocity $\sim$ 0.20 cm s-1 at the walls. The friction force was
found to be proportional to the applied velocity. The effective kinematic
viscosity was $\sim 0.1\kappa$, and effective Reynolds numbers within
$\mathrm{Re’}<15$. The fraction of the polarized vortex length varied between
zero in the middle of the channel and $\sim$ 60% within the shear flow regions
$\sim D/4$ from the walls. Therefore, we studied a state of polarized
ultraquantum (Vinen) turbulence fuelled at short lengthscales by vortex
reconnections, including those with vortex images due to the relative motion
between the vortex tangle and the pinning rough surface.

###### keywords:

superfluid helium, quantum turbulence, dissipative flow, vortex tangle,
pinning, friction

##  1 Introduction

Flow of superfluid 4He through a channel is only non-dissipative when its
velocity is below the channel-specific critical velocity $v_{c}$. Above it,
the dissipative regime sets in, in which the chaotic motion of quantum
vortices (called Quantum Turbulence) results in the transfer of the flow
momentum to the channel’s walls, i. e. a friction force.

Experimental evidence has shown that this friction force becomes greatly
reduced for $T<$ 0.7 K [1] and that vortex pinning becomes much weaker with
lowering temperature below $T<0.4$ K [2]. Approaching $T=0$, the density of
viscous normal component vanishes and the mutual friction which couples it to
the vortex tangle of turbulent superfluid is negligible, thus the direct
interaction between vortex lines and irregularities of the walls of the
channel (vortex pinning) must be considered.

Early evidence of quantized vorticity sticking to rough surfaces was observed
in 1958 by Hall and Shoenberg in torsional-oscillator experiments [3]. Further
studies involving rotating quantum turbulence [4], thermal counterflow [5] and
vortex capture probes such as wires [6] and MEMS devices [7] aid in providing
confirmation of the influence of surface roughness and vortex pinning on
superfluid flows. Numerical studies have examined vortex motion in the
presence of solid boundaries through use of the vortex filament model (VFM)
[8, 9, 10, 11, 12, 13, 14, 15, 16] and the Gross-Pitaevskii model [17].

In these computer simulations, we model a flow of superfluid 4He at $T=0$
between two parallel solid walls. We assume the limit of extremely rough
walls, where the areal density of sharp protuberances is greater than the
density of vortex lines in the vortex tangle. Then the processes of un-pinning
(due to a self-reconnection with an image vortex) of each vortex line from the
pinning protuberance and re-pinning at the nearest protuberance occur on the
lengthscale smaller than the characteristic scale of the vortex tangle (of
order mean inter-vortex distance $\ell$), i. e. independently of other
vortices. This is different from the the generation of dense vortex tangles by
larger-scale irregularities of the profile of solid wall simulated by Stagg et
al. [17].

In numerical simulations vortex lines are represented by chains of discrete
points with the set inter-point scale $\delta$, the smallest resolved
lengthscale is $\sim\delta$. Within the mechanism of pinning-unpinning
mentioned above, the elementary distance between the unpinned-pinned ends of a
vortex line becomes of order $\sim\delta$. Hence, in our framework, $\delta$
is a proxy of the scale of roughness of the solid wall.

##  2 Numerical Methods

The VFM [18] has been used to great effect to simulate and visualise the
dynamics of vortices in superfluid helium. In the limit of zero temperature,
using the VFM, the local self-induced vortex velocity $\dot{\mathbf{s}}$
evaluated at the point $\mathbf{r}$ along a vortex line may be described
entirely by the Biot-Savart integral

$\dot{\mathbf{s}}(\mathbf{r},t)=\frac{\kappa}{4\pi}\int_{\mathcal{L}}\frac{(\mathbf{s}-\mathbf{r})\times
d\mathbf{s}}{|\mathbf{s}-\mathbf{r}|^{3}},$ (1)

where $\kappa=h/m_{4}$ is the quantum of circulation in superfluid 4He. The
line integral over $\mathcal{L}$ represents inclusion of the complete vortex
configuration, which is discretized into points $\mathbf{s}_{i}$ for
$(i=1,...,N)$. The discretized integral, after removing the singularity by
separating the local and non-local contributions, becomes

$\dot{\mathbf{s}}_{i}=\frac{\kappa}{4\pi}\mathbf{s}_{i}^{\prime}\times\mathbf{s}_{i}^{\prime\prime}\text{ln}\bigg{(}\frac{2\sqrt{l_{+}l_{-}}}{e^{1/2}a}\bigg{)}+\frac{\kappa}{4\pi}\int_{r_{c}}\frac{(\mathbf{s}_{1}-\mathbf{s}_{i})\times
d\mathbf{s}_{1}}{|\mathbf{s}_{1}-\mathbf{s}_{i}|^{3}},$ (2)

where $l_{+}$ and $l_{-}$ represent the arc lengths to adjacent vortex points
$\mathbf{s}_{i+1}$, and $\mathbf{s}_{i-1}$, $\mathbf{s}_{i}^{\prime}$ and
$\mathbf{s}_{i}^{\prime\prime}$ are the local tangent and curvature
respectively, and $a\sim$ 1 Å is the size of the vortex core in He-II. Here
the first term is akin to the local induction approximation (LIA) and the
second term describes the non-local contributions, including the effect of
image-vortices due to any solid boundaries present [9]. To reduce the number
of computations required to simulate vortex tangles in a timely manner, the
integral has been reduced to only include contributions from points within a
cut-off radius $r_{c}$ of the target point $\mathbf{s}_{i}$, defining a sphere
of contributing non-local vortex points as shown in Figure 1 (right). The cut
radius was chosen to be half the container size $r_{c}=D/2$.

The simulation volume was a cubic cell $D\times D\times D$ with periodic
boundary conditions in the $x,y$ directions. The boundaries at $z=0$ and $z=D$
were modelled as a rough solid plane surface with strong pinning. To satisfy
the solid boundaries the method of images allows the duplication and
reflection of the vortex configuration across each solid boundary. This
combined with a periodic wrapping of the system at the $x,y$ boundaries gives
26 copies of the original volume: 2 reflected at solid boundaries, 8 simply
periodic and 16 periodic-reflected copies. An illustration of the boundary
conditions in the $x$ and $z$ directions is given in Fig. 1 (left). All 27
cubes were considered in evaluating Equation 2.

Figure 1: (Left) Visualisation of simulated boundary conditions. The
periodicity of the boundaries is shown only in the $x$ direction such that 9
out of 27 cubes are displayed. The original computational volume in the centre
is shown by the black box outline. Flat grey surfaces represent the solid
boundaries and are extended into the periodic boundaries in the $x$-direction.
The vortex tangle in each cube is coloured to indicate the type of
transformation applied: periodic-translation (red), solid-refelction (blue)
and a combination of the two (purple). (Right) Illustration of the sphere of
velocity contributions centred on a point, $\mathbf{s}_{i}$, near the corner
of the volume. Red vortex segments, within the sphere, contribute to
$\mathbf{\dot{s}}_{i}$ and vortex lines outside of the sphere are neglected.

The spatial resolution $\delta$ of the simulation is preserved by removing
points on the same vortex line that move within $\delta/2$ of each other.
Similarly, new points are introduced, according to the local curvature of the
line, when adjacent points are more than $\delta$ apart. Additionally, if two
non-adjacent points on the same vortex line or otherwise, become closer than
the critical distance $\delta/2$, a line-line reconnection occurs, altering
the vortex configuration and causing dissipation from the removal of vortex
line length. The reconnection procedure used was of Type-III [19].

To produce channel flow in the $T=0$ limit a positive finite superfluid
velocity $V$ is applied in the $x$ direction such that instantaneous velocity
of the discrete vortex point is

$\mathbf{u}_{i}=\dot{\mathbf{s}}_{i}+\begin{bmatrix}V\\\ 0\\\ 0\end{bmatrix}.$
(3)

The time-evolution followed a third order Adams-Bashforth scheme:

$\mathbf{s}_{i}^{n+1}=\mathbf{s}_{i}^{n}+\frac{\Delta
t}{12}\bigg{(}23\mathbf{u}_{i}^{n}-16\mathbf{u}_{i}^{n-1}+5\mathbf{u}_{i}^{n-2}\bigg{)}+\mathcal{O}(\Delta
t^{4})$ (4)

where $\Delta t$ is the size of a time step and $n=t/\Delta t$ is the current
time in integer steps [20].

The presence of solid wall boundaries alters the vortex motion since no flow
is allowed through the walls. The solid boundaries are assumed to be
microscopically rough such that the ends of vortex lines, terminated at the
walls, are permanently pinned and fixed in position with $\mathbf{u}(z=\pm
D/2)$ = 0. Liberation of the vortex line is only possible through its self-
reconnection with its image following the same reconnection criterion that
triggers a line-line event. Thus, a vortex point, which comes within a
distance $\delta$/4 of the wall, reconnects with its image, thus liberating
the line before it becomes pinned again $\sim 0.7\delta$ away, as demonstrated
in Figure 2. Here, the reconnection was triggered by the blue line segment,
which is then removed from the simulation alongside the adjacent red line
segment. The new end of the vortex is then instantaneously pinned by
artificially adjusting the $z$ coordinate to coincide with the surface, which
generates a kink in the filament just above the wall. The model mimics a
vortex line “walking” along a flat rough surface, with the vortex end jumping
between sharp protuberances spaced on the scale of the spatial resolution
$\delta$ [21].

Figure 2: Progressive snapshots of a single vortex line, which is initially
pinned vertically between two flat rough surfaces, evolving under an applied
velocity $V$ = 0.05 cm s-1. Vortex lines drawn in grey are shown at times: 50
ms, 100 ms, 150 ms, 200 ms, 255 ms and 260 ms. The region displayed is the
lower 5% of the channel with the solid green line marking the solid boundary.
The dashed green line shows the critical distance for vortex-wall
reconnection. The solid black lines show the filament in the instant before
and after the first time the vortex reconnected with the wall at $t=253.69$
ms. Coloured line segments (red and blue) show the length removed by the
event, and the blue segment is responsible for the reconnection. The size of
arrows indicates the spatial resolution, $\delta$ and the distance between new
and old pinning sites. The inset illustrates the moment of reconnection, with
the image-vortex drawn symmetrically behind the boundary.

Effectively, this mechanism assumes the existence of a critical angle between
the vortex and the wall for unpinning. Such an assumption was recently used in
simulations of vortex lines attached to a microelectromechanical (MEMS)
oscillator, with a chosen critical angle of $\theta_{c}=\pi/6$ [16].

The potential benefits of our approach are in that one could investigate the
role of Kelvin waves emitted on vortex lines by discrete steps; although to
assure an adequate resolution of those waves one might need to introduce a
separate lengthscale for the kink generated by an image-reconnection such that
the wavelength of injected distortions is larger than the resolution.

To examine the rate of momentum transfer to the walls, two equivalent methods
of determining the friction force were used. Both have been validated by
comparison with the analytical solution for a vortex semi-ring travelling
along a flat rough surface. The rate of change of the total impulse of all
vortex lines in the simulation volume, termed the integral method:

$\mathbf{F}=\rho_{s}\kappa\int(\dot{\mathbf{s}}\times\mathbf{s}^{\prime})d\xi,$
(5)

where $\rho_{s}$ is the superfluid density and $\xi$ is the arc length along a
vortex line. The second method assumes that the vortex line tension $f_{t}$
(energy per unit line length) is constant. The angle $\theta$ between the
surface and the pinned segment can be used to obtain the components of the
line tension for a single terminated line. Thus, summing over all the pinned
ends gives the total friction force in the streamwise direction

$F_{x}=-f_{t}\sum_{i,pinned}\text{cos}\theta_{i}=-\frac{\rho\kappa^{2}}{4\pi}\text{ln}\Big{(}\frac{b}{a}\Big{)}\sum_{i,pinned}\text{cos}\theta_{i},$
(6)

where $b$ was a cut-off length scale chosen as $\delta/2$ which is the
effective size of pinning sites in the vortex walking and $a=1$ Å. We term
this the tension method.

To examine the velocity profile, the instantaneous mean velocity as a function
of the wall-normal direction $z$ must be calculated. This was achieved by
slicing the simulation volume into $N$ slabs with thickness $D/(N-1)$, the
mean velocity of all line segments within a single slab was then time-averaged
during the steady state to obtain the coarse-grained velocity profile
$\langle\mathbf{u}(z)\rangle$. The component in the direction of flow,
$\langle u_{x}\rangle$ then gives the cross-channel velocity profile. Fixed
segments of the vortex line were omitted from such calculations to avoid
skewing the average velocity, with $u_{pinned}=0$, in the top and bottom
slabs. The variation in the coarse-grained vortex line density across the
height of the channel $\langle\mathcal{L}(z)\rangle$ was obtained similarly,
but without the need to exclude fixed segments.

##  3 Results

The parameters for simulations were: $D=0.1$ cm, $\delta=2\times 10^{-3}$ cm,
with temporal resolution $dt=8\times 10^{-5}$ s. Initially, at $t=0$, the
volume was populated with 80 randomly placed and randomly oriented vortex
rings with equivalent radii of 0.012 cm; with an instantly acting applied
superflow of $V=0.30$ cm s-1. Snapshots of the early development of the tangle
is shown by the inset panels of Fig. 3. At early times ($t<1$ s) frequent
reconnections between vortex rings in the centre of the channel cause an
immediate drop in the total vortex line length, $\Lambda$, of the system,
which is also shown in Fig. 3. Gradually vortex rings become attached to the
walls and begin “walking” across the surface in the direction of flow. Within
$t=4$ s both surfaces are decorated with several terminating vortices.

Figure 3: Vortex line length as a function of time shown for 5 different
simulations. Solid lines show data and dashed lines show fitted curves. The
inset plots are early-time snapshots of the developing vortex tangle for
$V=0.30$ cm s-1.

As the flow develops further vortex lines become stretched along it causing a
gradual increase in the total line length as a vortex tangle forms between the
walls. Both line-line and image-line reconnections occur frequently, each
process removing a small portion of vortex length, which becomes approximately
stable when the rate of dissipating length from reconnections is similar to
the rate of vortex growth at the walls. At $t=80$ s the vortex tangle is
comfortably in the steady state with the total vortex length in the simulation
volume $\Lambda=4.14$ cm, and the walls are decorated with 14 pinned vortex
lines. The tangle at $t=80$ s is the same as the configuration shown in Fig.
1.

To simulate flows at any other applied velocity in the range 0.08 cm
s${}^{-1}\leq V\leq$ 0.34 cm s-1 for further 80 s, the same initial
configuration of the vortex tangle was used, namely, the one for $V=0.30$ cm
s-1 at $t=80$ s.

The evolution of $\Lambda$ for a selection of $V$ is shown in Fig. 3 in the
range 80 s $\leq t\leq$ 160 s. Vortex tangles, sustained for at least 80 s as
in Fig. 4, were observed above the critical applied flow velocity
$V_{c}=0.19\pm 0.01$ cm s-1, i. e. $V_{c}D/\kappa\sim 20$. Below $V_{c}$, the
tangle is not sustainable; all vortex lines eventually detach from the walls
at some point and move with the flow without gaining any energy – hence the
tangle decays.

Figure 4: Snapshots of vortex tangles in the steady state at $t=$ 160 s for
all $V>V_{c}$. The flow is directed towards the right.

For each applied superfluid velocity $V$, the mean flow velocity was
calculated as $\langle u_{x}\rangle=D^{-1}\int_{0}^{D}u_{x}(z)dz$ and is
plotted against $V$ in Fig. 5a. They are proportional to each other with
relation $\langle u_{x}\rangle\approx 0.96V$.

Figure 5: (a) Time-averaged mean flow velocity as a function of the applied
superflow shown for both sustained (closed symbols) and unsustained (open
symbols) turbulence. (b) The coarse-grained friction force exerted by vortex
lines on the channel walls, normalised by a constant vortex line tension,
shown to be proportional to the mean flow. (c) Time-averaged vortex line
density as a function of the mean flow. Open circles represent unsustainable
line densities.

The stream-wise component of friction force was evaluated from equations (5)
and (6), time-averaged over the steady state and normalised by the reduced
vortex line tension
$f_{t}=\frac{\rho\kappa^{2}}{4\pi}\text{ln}\Big{(}\frac{b}{a}\Big{)}$, where
$b=\delta/2$ is the effective size of a pinning site. The observed force-
velocity relation, shown by Fig. 5b, was approximately linear within the range
of flow velocities selected. Friction force per number of pinned vortex ends
was constant across all velocities at $(0.626\pm 0.007)f_{t}$ demonstrating
that the total friction force is proportional to the number of pinned ends in
the simulation box and the effective critical angle in the direction of flow
is $\gtrsim\cos^{-1}(0.626)=51^{\circ}$. Both integral and tension methods of
determining the friction force consistently agreed to within 1%.

The total vortex line length $\Lambda$ fluctuates with the dominant period
$\sim 10$ s (corresponding to the time constant $\sim 2$ s). This could be
compared with the time $D/V_{c}\sim 0.5$ s for the mean flow to pass the cell
size, the intrinsic quantum time constant $D^{2}/\kappa=10$ s and our chosen
observation time of 80 s. Fig. 5c demonstrates the variation of time-averaged
vortex line density, $\langle\mathcal{L}\rangle=\langle\Lambda\rangle/D^{3}$,
which appears to be cubic with velocity. Interestingly, the values of
$\langle\mathcal{L}\rangle$ for unsustained tangles do not fall far from the
$V^{3}$ line.

The coarse-grained stream-wise velocity profile $\langle u_{x}(z)\rangle$ was
steady and nearly parabolic albeit with a non-zero slip velocity, as seen in
Fig 6. We fit the profile with the form of a classical Poiseuille profile,
adapted for finite slip boundary conditions at the walls,

$\langle
u_{x}(z)\rangle=u_{max}-4(u_{max}-u_{min})\left(\frac{z}{D}-\frac{1}{2}\right)^{2}\text{,}$
(7)

where $u_{min}$ represents the slip-velocity at the walls and $u_{max}$ is the
maximum local velocity in the centre of the channel. The agreement is
relatively good, but the fits consistently underestimate $u_{max}$. The slip-
velocity varied only weakly across the simulations with values in the region
of $V_{c}$ as seen in Fig. 7 (left). The universality of the flow profiles is
shown in Fig. 6 (right). Deviations from the universal profile were only
observed for the lowest $V$ which still produced a steady flow - possibly due
to the smaller number of line segments leading to less averaging off of
fluctuations.

Figure 6: Cross-channel coarse-grained velocity profiles. (Left) Absolute
values, reduced by the slip-velocity at the walls. (Right) Normalised by the
range of velocity present in the flow. Figure 7: (Left) Results of fitting
velocity profiles with equation (7) as function of $\langle u_{x}\rangle$.
Dashed line follows $u_{max}\propto\langle u_{x}\rangle$. (Right) The
effective kinematic viscosity as function of $\langle u_{x}\rangle$. The
vicinity of $V_{c}$ is shown by the shaded area in both graphs.

The effective kinematic viscosity $\nu^{\prime}=\frac{\langle
F\rangle}{8D(u_{max}-u_{min})\rho}$ was in the range 0.3$\kappa$ – 0.1$\kappa$
– shown in Fig. 7 (right) – meaning the effective Reynolds number
$\text{Re}^{\prime}=D(V-V_{c})/\nu^{\prime}$ was between 0 and $\sim 15$.
Clearly this was insufficient to support quasi-classical turbulence in the
coarse-grained velocity field. Thus, we have polarised ultra-quantum
turbulence driven by injections of short-wavelength Kelvin waves, enabled by
the frequent reconnection of vortex ends due to the relative motion between
the vortex tangle and the rough wall.

The coarse-grained vortex line density profile, shown in Fig. 8a, demonstrates
the tangle formed maximally dense regions $\sim D/4$ away from each surface,
with minima in the centre of the channel and at the walls. Such a spatial
variation in $\langle\mathcal{L}\rangle$ is similar with simulations at finite
temperature where flows were driven by a normal fluid flow described by a
Poiseuille profile [22]. The peaks are mirrored by the polarised line density
$\langle\mathcal{L}_{y}(z)\rangle$ (Fig. 8b) and their ratio
$\frac{\langle\mathcal{L}_{y}(z)\rangle}{\langle\mathcal{L}(z)\rangle}$ (Fig.
8c) shows the polarisation of the vortex tangle, whose total was between 0 and
60% depending on the distance from walls, in all simulations.

Figure 8: (a) Coarse-grained vortex line density profile. (b) Polarised vortex
line density profile. (c) Time-averaged polarisation of the vortex tangle.

##  4 Conclusions

These simulations demonstrate a critical velocity $V_{c}=0.19\pm 0.01$ cm
s${}^{-1}\simeq 20\frac{\kappa}{D}$ for pure superfluid flow in channels,
above which the turbulent state of the vortex tangle sustained for at least 80
s. There is a finite slip velocity at the solid boundary, whose value is
similar to $V_{c}$, resultant from the “walking” of terminated vortex lines
along the rough surface due to frequent reconnections with their images. The
coarse-grained velocity profiles $u_{x}(z)$ were similar to classical laminar
profiles (albeit with a non-zero velocity at the walls) – consistent with the
effective Reynolds number $\text{Re}^{\prime}\leq 15$ too small to sustain
quasi-classical turbulence that typically requires $\text{Re}>3000$. We thus
had ultra-quantum (Vinen’s) regime of quantum turbulence. The fraction of the
polarised vortex length reached $\sim 60\%$ within the shear regions near both
walls, and decreased to zero in the middle of the channel. The effective
kinematic viscosity $\nu^{\prime}$, computed for the shear regions, tends
towards $\sim 0.1\kappa$ with increasing the applied flow velocity – in
quantitative agreement with experimental observations for quantum turbulence
in superfluid 4He [1]. The flows exhibited an approximately proportional
dependence of the friction force $F_{x}(V)$ and near cubic dependence of the
vortex line length $\Lambda(V)$ on the mean flow velocity $V$. These
dependences held even below $V_{c}$ – for metastable vortex tangles of
lifetime shorter than 80 s.

Acknowledgments The authors would like to acknowledge the assistance given by
Research IT and the use of the Computational Shared Facility at The University
of Manchester. This work was supported by the Engineering and Physical
Sciences Research Council [grant number EP-T517823-1].

Data Availability The code used to generate these simulations and guidance in
its use can be accessed from the authors upon reasonable request.

## References

  *   * Zmeev et al. [2015] Zmeev, D.E., Walmsley, P.M., Golov, A.I., McClintock, P.V.E., Fisher, S.N., Vinen, W.F.: Dissipation of quasiclassical turbulence in superfluid ${}^{4}\mathrm{He}$. Phys. Rev. Lett. 115, 155303 (2015) https://doi.org/10.1103/PhysRevLett.115.155303
  * Donev et al. [2001] Donev, L.A., Hough, L., Zieve, R.J.: Depinning of a superfluid vortex line by $\mathrm{K}$elvin waves. Phys. Rev. B 64, 180512 (2001) https://doi.org/%****␣sn-article.bbl␣Line␣75␣****10.1103/physrevb.64.180512
  * Hall and Shoenberg [1958] Hall, H.E., Shoenberg, D.: An experimental and theoretical study of torsional oscillations in uniformly rotating liquid helium $\mathrm{II}$. Proc. Math. Phys. Eng. Sci. Series A. Mathematical and Physical Sciences 245(1243), 546–561 (1958) https://doi.org/10.1098/rspa.1958.0100
  * Adams et al. [1985] Adams, P.W., Cieplak, M., Glaberson, W.I.: Spin-up problem in superfluid ${}^{4}\mathrm{He}$. Phys. Rev. B 32, 171–177 (1985) https://doi.org/10.1103/PhysRevB.32.171
  * Hegde and Glaberson [1980] Hegde, S.G., Glaberson, W.I.: Pinning of superfluid vortices to surfaces. Phys. Rev. Lett. 45, 190–193 (1980) https://doi.org/10.1103/PhysRevLett.45.190
  * Neumann and Zieve [2014] Neumann, I.H., Zieve, R.J.: Vortex pinning by surface geometry in superfluid helium. Phys. Rev. B 89, 104521 (2014) https://doi.org/10.1103/PhysRevB.89.104521
  * Barquist et al. [2020] Barquist, C.S., Jiang, W.G., Gunther, K., Eng, N., Lee, Y., Chan, H.B.: Damping of a microelectromechanical oscillator in turbulent superfluid ${}^{4}\mathrm{He}$: A probe of quantized vorticity in the ultralow temperature regime. Phys. Rev. B 101, 174513 (2020) https://doi.org/%****␣sn-article.bbl␣Line␣150␣****10.1103/PhysRevB.101.174513
  * Schwarz [1981] Schwarz, K.W.: Vortex pinning in superfluid helium. Phys. Rev. Lett. 47, 251–254 (1981) https://doi.org/10.1103/PhysRevLett.47.251
  * Schwarz [1985] Schwarz, K.W.: Three-dimensional vortex dynamics in superfluid ${}^{4}\mathrm{He}$: Line-line and line-boundary interactions. Phys. Rev. B 31, 5782 (1985) https://doi.org/10.1103/PhysRevB.31.5782
  * Schwarz [1988] Schwarz, K.W.: Three-dimensional vortex dynamics in superfluid ${}^{4}\mathrm{He}$: Homogeneous superfluid turbulence. Phys. Rev. B 38, 2398–2417 (1988) https://doi.org/10.1103/PhysRevB.38.2398
  * Schwarz [1992] Schwarz, K.W.: Effect of surface roughness on the critical velocities of superfluid ${}^{4}\mathrm{He}$. Phys. Rev. Lett. 69, 3342–3345 (1992) https://doi.org/10.1103/PhysRevLett.69.3342
  * Schwarz [1993] Schwarz, K.W.: Unwinding of a single quantized vortex from a wire. Phys. Rev. B 47, 12030–12039 (1993) https://doi.org/10.1103/PhysRevB.47.12030
  * Tsubota and Maekawa [1993] Tsubota, M., Maekawa, S.: Pinning and depinning of two quantized vortices in superfluid ${}^{4}\mathrm{He}$. Phys. Rev. B 47, 12040–12050 (1993) https://doi.org/10.1103/PhysRevB.47.12040
  * Tsubota [1994] Tsubota, M.: Capacity of a pinning site for trapping quantized vortices in superfluid ${}^{4}\mathrm{He}$. Phys. Rev. B 50, 579–581 (1994) https://doi.org/10.1103/PhysRevB.50.579
  * Fujiyama and Tsubota [2009] Fujiyama, S., Tsubota, M.: Drag force on an oscillating object in quantum turbulence. Phys. Rev. B 79, 094513 (2009) https://doi.org/10.1103/PhysRevB.79.094513
  * Nakagawa et al. [2023] Nakagawa, T., Tsubota, M., Gunther, K., Lee, Y.: Dynamics of pinned quantized vortices in superfluid ${}^{4}\mathrm{He}$ in a microelectromechanical oscillator. Phys. Rev. B 108, 144110 (2023) https://doi.org/10.1103/PhysRevB.108.144110
  * Stagg et al. [2017] Stagg, G.W., Parker, N.G., Barenghi, C.F.: Superfluid boundary layer. Phys. Rev. Lett. 118, 135301 (2017) https://doi.org/10.1103/PhysRevLett.118.135301
  * Hänninen and Baggaley [2014] Hänninen, R., Baggaley, A.W.: Vortex filament method as a tool for computational visualization of quantum turbulence. Proc. Natl. Acad. Sci. 111, 4667–4674 (2014) https://doi.org/10.1073/pnas.1312535111 https://www.pnas.org/doi/pdf/10.1073/pnas.1312535111
  * Baggaley [2012] Baggaley, A.W.: The sensitivity of the vortex filament method to different reconnection models. J. Low Temp. Phys. 168(1-2), 18–30 (2012) https://doi.org/10.1007/s10909-012-0605-8
  * Baggaley and Barenghi [2012] Baggaley, A.W., Barenghi, C.F.: Tree method for quantum vortex dynamics. J. Low Temp. Phys. 166, 3–20 (2012) https://doi.org/10.1007/s10909-011-0405-6
  * Golov et al. [2021] Golov, A.I., Doyle, M.J., Walmsley, P.M., Baggaley, A.W.: Modeling of the interaction of a vortex line with a solid wall in the $\mathrm{T}=0$ limit. Poster presented at the International Conference on Quantum Fluids and Solids (2021)
  * Baggaley and Laurie [2015] Baggaley, A.W., Laurie, J.: Thermal counterflow in a periodic channel with solid boundaries. J. Low Temp. Phys. 178(1-2), 35–52 (2015) https://doi.org/10.1007/s10909-014-1226-1