# RankSHAP: a Gold Standard Feature Attribution Method for the Ranking Task

Tanya Chowdhury  Yair Zick  James Allan

###### Abstract

Several works propose various post-hoc, model-agnostic explanations for the
task of ranking, i.e. the task of ordering a set of documents, via feature
attribution methods. However, these attributions are seen to weakly correlate
and sometimes contradict each other. In classification/regression, several
works focus on _axiomatic characterization_ of feature attribution methods,
showing that a certain method uniquely satisfies a set of desirable
properties. However, no such efforts have been taken in the space of feature
attributions for the task of ranking. We take an axiomatic game-theoretic
approach, popular in the feature attribution community, to identify candidate
attribution methods for ranking tasks. We first define desirable axioms: Rank-
Efficiency, Rank-Missingness, Rank-Symmetry and Rank-Monotonicity, all
variants of the classical Shapley axioms. Next, we introduce Rank-SHAP, a
feature attribution algorithm for the general ranking task, which is an
extension to classical Shapley values. We identify a polynomial-time algorithm
for computing approximate Rank-SHAP values and evaluate the computational
efficiency and accuracy of our algorithm under various scenarios. We also
evaluate its alignment with human intuition with a user study. Lastly, we
theoretically examine popular rank attribution algorithms, EXS and Rank-LIME,
and evaluate their capacity to satisfy the classical Shapley axioms.

Machine Learning, ICML

##  1 Introduction

A user submits a query to a search engine, wanting to know, say, what is the
best car to purchase. The search engine processes the query and outputs a
ranked list of documents (web pages) in response to the query. How did the
search engine decide on the order in which the documents are presented? What
features cause a certain document to be ranked higher than another? It is
reasonable to assume we have neither the capacity to access the internal
workings of the predictive model, nor that we could even understand the
underlying decision-making process were we given such access. Thus, a standard
assumption in the literature is that the model is a _black-box_ , and we only
have limited query access to it.

The task of ranking entails a query $\vec{q}$ and a set of documents $D$ as
input to a ranking model $f_{R}$, which in turn generates an ordering of these
documents. We assume that the ranking model is available only as a black-box
which can be queried using an API. For a particular query, we wish to
understand the logic behind the generated ordering, to verify the
functionality of the model as well as to help end users. Ranking feature
attributions are critical for building trust in systems which are used for
high stakes decision making (Anand et al., 2022; Chowdhury et al., 2023). Such
feature attributions are also useful in ruling out prejudices in widely used
product search and recommendation systems (Singh & Anand, 2019; Singh et al.,
2020).

Additive feature attributions have become a well studied topic in recent
years, especially for regression and classification tasks. To maintain a
consistency between different feature attribution methods, various authors
(Lundberg & Lee, 2017; Datta et al., 2016) proposed a unique solution based on
the Shapley value (Shapley, 1953). This solution is considered the gold
standard of feature attributions, but is NP-hard to compute. As a result,
different approximations such as KernelSHAP (Lundberg & Lee, 2017) and
DeepSHAP (Lundberg & Lee, 2017; Shrikumar et al., 2017) are used for popular
real-world applications.

Quite a few attempts have been made to generate feature attributions for
ranking tasks, often by extending attribution methods from
classification/regression. EXS (Singh & Anand, 2019) and Rank-LIME (Chowdhury
et al., 2023) extend LIME in different manners to obtain ranking feature
attributions. However, these attributions contradict each other across the
different settings in the very same algorithm. Fernando et al. (2019) also
generate explanations for ranking models using DeepSHAP and report little
correlation between the attributions generated by different reference values
for the very same query. This shows us that extensions of popular feature
attribution methods for classification and regression are not good attribution
candidates for ranking. As a result, defining a gold standard in ranking
feature attributions becomes a necessity.

##### Our Contributions.

In Section 2, we briefly discuss the Shapley value, define the ranking
attribution problem, and describe the NDCG metric for measuring the
effectiveness of an ordering. In Section 3, we first define the axioms of
Rank-Efficiency, Rank-Missingness, Rank-Symmetry and Rank-Monotonicity, based
on the popular NDCG rank-evaluation metric, which we believe a standard
attribution method for the ranking task must satisfy. Next, we introduce
RankSHAP, which is Shapley values with the NDCG characteristic function, that
satisfies the above properties. We then find computationally efficient
approximations to RankSHAP for different scenarios. In Section 4, we evaluate
the fidelity based performance of our algorithm on the BM25, BERT and T5
ranking models with query-document sets from the MS Marco dataset in a variety
of settings. We observe that RankSHAP outperforms the best competing system by
$30.78\%$ on Fidelity and $23.68\%$ on wFidelity. We also conduct a user study
to understand how effective RankSHAP is towards helping end users and find
that RankSHAP feature attributions can guide participants to reorder documents
and estimate the query quite effectively. Finally, in Section 5, we discuss
alternatives to our chosen value function and weigh their pros and cons. We
follow up by analyzing whether the rank attribution algorithms of EXS and
RankLIME comply to the Shapley axiom framework. Ultimately, we suggest that
RankSHAP be used as the gold standard for feature attributions in ranking
related tasks. In the spirit of reproducibility, all our experiment scripts
and user study data are made available.

##  2 Background and Preliminaries

We first briefly describe the Shapley axioms and the Shapley value, as used
for credit assignment in coalitional game theory. Next we introduce the
ranking feature attribution problem and discuss two existing solutions to it.

###  2.1 Basic Axioms and the Shapley Value

The Shapley axioms for feature attributions have been borrowed from
coalitional game theory, where they are used to fairly distribute
costs/rewards among a set of players. We are given a classification/regression
model $f$. The feature attribution of feature $i$ (can be words, human
engineered features etc. ) in the decision made by $f$, towards the model
input $\vec{x}$ is represented by $\phi_{i}(f,\vec{x})$.

We adopt Young’s characterization of the Shapley value, also used by Datta et
al. (2016). Let us assume there are a total of $m$ features. For a feature
attribution $\phi_{i}(f,\vec{x})$ to be considered valid, it is necessary for
it to satisfy certain axioms. We briefly describe the axioms below:

Efficiency:

The sum of feature attributions over a set of features should be equal to the
difference between the model output containing all features and the output
containing no features.
$\sum_{i=1}^{m}\phi_{i}(f,\vec{x})=f(\vec{x})-f(\phi)$.

Symmetry:

If for every subset of features $S$ s.t features $i,j\not\in S$, the effects
of adding $i$ to $S$ is equal to the effect of adding $j$ to $S$, then
features $i$ and $j$ should be assigned equal attribution towards
$f(\vec{x})$.

Missingness:

If for every subset of features $S$ s.t features $i\not\in S$, $f(S\cup
i,\vec{x})$ is equal to $f(S,\vec{x})$, then $\phi_{i}(f,\vec{x})$ must be
given an attribution of $0$. $f(S,\vec{x})$ represents the case where features
not in set S, have been replaced by their reference value, and model input
$\vec{x}$ has been reconstructed.

Monotonicity:

Given two models $f$ and $f^{\prime}$, if for all feature coalitions
$S,i\not\in S$, $f(S\cup i,\vec{x})-f(S,\vec{x})$ is greater than or equal to
$f^{\prime}(S\cup i,\vec{x})-f^{\prime}(S,\vec{x})$, then
$\phi_{i}(f,\vec{x})$ is greater than or equal to
$\phi_{i}(f^{\prime},\vec{x})$.

The _Shapley value_ (Young, 1985; Shapley, 1953) uniquely satisfies the four
axioms above, and has been widely adopted in feature attribution for
regression/classification tasks. It is defined as:

$\displaystyle\phi_{i}(f,\vec{x})$
$\displaystyle=\sum_{\vec{z}\subseteq\vec{x}}\frac{|\vec{z}|!(m-|\vec{z}|-1)!}{m!}[f_{\vec{x}}(\vec{z})-f_{\vec{x}}(\vec{z}\setminus
i)]$

where $|\vec{z}|$ is the number of non-zero entries in $\vec{z}$, and
$\vec{z}\subseteq\vec{x}$ represents all $\vec{z}$ vectors where the non-zero
entries are a subset of the non-zero entries in $\vec{x}$.

###  2.2 The Ranking attribution Problem

Let $f_{R}$ represent a ranking model which, given a query $\vec{q}$ and a set
of documents $D=\\{\vec{d}_{1},\dots,\vec{d}_{n}\\}$, returns an ordering of
documents in $D$, i.e., $f_{R}(\vec{q},D)$. The ordering model $f_{R}$ is only
accessible to us as a black-box. The query and set of documents $(\vec{q},D)$
jointly form the instance ($\vec{x}$), and the goal is to generate post-hoc
feature attributions for $f_{R}(\vec{x})$. Documents here can be represented
in any manner e.g.: bag of words, human engineered features etc. The classical
Shapley framework cannot be directly applied here, since $f_{R}(\vec{x})$ is
not numerical, but a permutation of $1\dots n$.

###  2.3 Normalized Discounted Cumulative Gain (NDCG)

Normalized Discounted Cumulative Gain (NDCG) (Järvelin & Kekäläinen, 2002) is
the most popular metric to evaluate the effectiveness of an ordering, often
used in search engines and e-commerce. Later in this work, we use the NDCG
score as our characteristic function, to quantify the value of a coalition of
features. It assumes that some relevance score $\textit{rel}_{j}$ for each
document $d_{j}$ is known, and that in an effective ordering highly relevant
documents appear at the beginning of the list. For an ordered list of $n$
documents, the Discounted Cumulative Gain is computed as:

$\displaystyle\textit{DCG}_{n}$
$\displaystyle=\sum_{j=1}^{n}\frac{\textit{rel}_{j}}{\log_{2}{(j+1)}}$

The above expression is dependent on the number of documents in the set being
evaluated. This metric is thus difficult to use as-is, when comparing document
lists of different lengths. As a result, we normalize it by dividing it with
the maximum possible $\textit{DCG}_{n}$ score, also known as
Ideal-$\textit{DCG}_{n}$($\textit{IDCG}_{n}$). This is obtained by sorting
documents in decreasing order of relevance and computing $\textit{DCG}_{n}$
using the so obtained ideal ordering. To compute $\textit{NDCG}_{n}$, we
simply normalize $\textit{DCG}_{N}$:
$\textit{NDCG}_{N}=\frac{\textit{DCG}_{n}}{\textit{IDCG}_{n}}$. Note that
$\textit{NDCG}_{n}\in[0,1]$, and that the higher the value, the more effective
the ranking. In situations where document relevance scores $\textit{rel}_{j}$
are unknown, they can be estimated using implicit measures such as clicks,
views or the time spent on a page.

##  3 The RankSHAP Framework

We first present ranking-based variants of the Shapley axioms: Rank-
efficiency, Rank-missingness, Rank-symmetry and Rank-monotonicity. We then
show that our method, RankSHAP, a redefinition of Shapley values with the NDCG
characteristic function, satisfies the axioms listed above. We also provide
computationally feasible approximations of the RankSHAP expression (extending
KernelSHAP), such that it can be used in practice.

###  3.1 Ranking Axioms

The ranking feature attribution problem is the process of assigning $n$
documents with $m$ features each, a $m$ dimensional feature attribution vector
$\phi_{R}(f_{R},\vec{x})$, explaining the ordering output by a ranking
function $f_{R}$ over a query $\vec{q}$ and set of documents $D$. Let
$\vec{x}$ represent the input to the ranking model $f_{R}$, which is some
simplified representation of ($\vec{q}$, D) like bag of words, human
engineered features etc. Let the set of features be represented by $M$.

Rank-Efficiency:

The sum total of feature attributions for a particular ordering should be
equal to the difference between the NDCG value of the ordering obtained when
all features are included and the NDCG value of the ordering obtained when no
features are included.

$\displaystyle\sum_{i=1}^{m}\phi_{R}(f_{R},\vec{x},i)=\textit{NDCG}(f_{R}(\vec{x}))-\textit{NDCG}(f_{R}(\emptyset))$

Rank-Missingness:

If a feature $i$ does not change the NDCG score for any ranking computed based
on a coalition $S\subseteq M\setminus\\{i\\}$, it has to be assigned a value
of $0$. More formally, if for all $S\subseteq M\setminus\\{i\\}$ we have
$\textit{NDCG}(f_{R}(S\cup i,\vec{x}))=\textit{NDCG}(f_{R}(S,\vec{x}))$, then
$\phi_{R}(f_{R},\vec{x},i)=0$.

The Rank-missingness axiom naturally extends the null/dummy axiom.

Rank-Symmetry:

two features $i,j\in M$ are called _symmetric_ if they produce equal NDCG
scores when added individually to every possible feature coalition $S\subseteq
M\setminus\\{i,j\\}$. The Rank-symmetry axiom requires that we assign equal
values to symmetric features, i.e. If for all $S\subseteq M\setminus\\{i,j\\}$
we have $\textit{NDCG}(f_{R}(S\cup i,\vec{x}))=\textit{NDCG}(f_{R}(S\cup
j,\vec{x}))$, then $\phi_{R}(f_{R},\vec{x},i)=\phi_{R}(f_{R},\vec{x},j)$.

Rank-Monotonicity:

For input $\vec{x}$ and two ranking models $f_{R}$ and $f_{R}^{\prime}$, if
for all coalitions $S$ where $S\subseteq M\setminus\\{i\\}$, the marginal
contribution of adding $i$ to $S$ in $f_{R}$, is greater than or equal to the
marginal contribution of adding $i$ to $S$ in $f_{R}^{\prime}$, the
attribution of $i$ in $f_{R}$ has to be greater in $f_{R}$ as compared to
$f_{R}^{\prime}$. If for all $S\subseteq M\setminus\\{i\\}$ we have
$\textit{NDCG}(f_{R}(S\cup
i,\vec{x}))-\textit{NDCG}(f_{R}(S,\vec{x}))\geq\textit{NDCG}(f_{R}^{\prime}(S\cup
i,\vec{x}))-\textit{NDCG}(f_{R}^{\prime}(S,\vec{x}))$ then
$\phi_{R}(f_{R},\vec{x},i)\geq\phi_{R}(f_{R}^{\prime},\vec{x},i)$.

Theorem 3.1 is a direct corollary of Young (1985).

######  Theorem 3.1.

There exists a unique ranking feature attribution method that jointly
satisfies the properties of Rank-Efficiency, Rank-Missingness, Rank-Symmetry
and Rank-Monotonicity. It is given by $\phi_{R}$.

$\displaystyle
m_{S}=\textit{NDCG}(f_{R}(S\cup\\{i\\},\vec{x}))-\textit{NDCG}(f_{R}(S,\vec{x}))$
$\displaystyle\phi_{R}(f_{R},\vec{x},i)=\sum_{S\subseteq
M\setminus\\{i\\}}\frac{|S|!\times(m-|S|-1)!}{m!}\times m_{S}$

Hence we present RankSHAP, a valid feature attribution for our ranking
function $f_{R}$ that satisfies several natural properties. Computing the
Shapley value is NP-Hard (Shapley, 1953) in the form presented above, as it is
exponential to the number of features. NDCG computation for each subset
further adds to the computation as a function of the number of documents by a
factor of $O(nlogn)$. Moreover, the NDCG metric is dependent on relevance
labels, which are at times not available. In the next section, we discuss how
we approximate the RankSHAP solution expression to be able to use it in day to
day applications.

###  3.2 Approximating RankSHAP

We look into approximations of RankSHAP and discuss viable implementations.

Inferring NDCG:

In the RankSHAP formulation, we compute the value of every coalition by
calculating the NDCG value of the generated ordering. The NDCG score relies on
knowing relevance scores for each document, which can ideally be obtained from
ground truth labels of a human annotated dataset. However, ground truth labels
are not always available. Several neural rankers (e.g.: BERT (Nogueira & Cho,
2019)) assign a similarity score to each query-document pair, which are then
aggregated to generate an ordering. In situations where a similarity score is
available per document, this score can be used to infer relevance for that
document and ultimately compute NDCG. In situations where a similarity score
is not available, there are methods to indirectly infer relevance scores.
Common methods in web search and product recommendation include click-through
rate, time spent on page, bounce rate, scroll depth, bookmarking, times shared
on social media etc (Zoeter et al., 2008).

Kernel-RankSHAP:

While the Shapley value itself is NP-hard to compute, KernelSHAP (Lundberg &
Lee, 2017) leverages a kernel based model to faithfully approximate the
Shapley value. It is a linear approximation of the ranking model $f_{R}$ at
instance $\vec{x}$, which extends the LIME optimization function with
following kernel and loss function. Below we customize it for ranking
attributions and name it Kernel-RankSHAP. Let $G$ be the class of all linear
additive attributions.

$\displaystyle\phi_{R}(f_{R},\vec{x},i)$ $\displaystyle=\arg\min_{g\in
G}L(f_{R},g,\pi_{\vec{x}})$ $\displaystyle L(f_{R},g,\pi_{\vec{x}})$
$\displaystyle=\sum_{\vec{z}\in
Z}[\textit{NDCG}(f_{R}(\vec{z}))-\textit{NDCG}(g(\vec{z}))]^{2}\pi_{\vec{x}}(\vec{z})$
$\displaystyle\pi_{\vec{x}}(\vec{z})$
$\displaystyle=\frac{(m-1)}{{m\choose{|\vec{z}|}}|\vec{z}|(m-|\vec{z}|)}$

Table 1:  Comparing performance between different feature attribution methods (Random, EXS, RankLIME and RankSHAP) for word-based explanations of textual rankers (BM25, BERT and T5) based on model fidelity (Fidelity) and weighted-fidelity (wFidelity). We generate feature attributions for each ranking models with the 10, 20 and 100 most relevant documents for a particular query. Ranker $\downarrow$ |  | Top-10 | Top-20 | Top-100
---|---|---|---|---
Metric $\rightarrow$ |  | Fidelity | wFidelity | Fidelity | wFidelity | Fidelity | wFidelity
BM25 | Random | 0.08 | 0.04 | -0.05 | 0.01 | -0.09 | 0.00
EXS | 0.39 | 0.34 | 0.31 | 0.29 | 0.24 | 0.17
RankLIME | 0.45 | 0.37 | 0.37 | 0.29 | 0.26 | 0.20
RankSHAP | 0.63 | 0.48 | 0.54 | 0.40 | 0.47 | 0.33
BERT | Random | 0.12 | 0.02 | -0.09 | 0.01 | -0.12 | 0.00
EXS | 0.36 | 0.21 | 0.28 | 0.18 | 0.11 | 0.08
RankLIME | 0.41 | 0.24 | 0.32 | 0.24 | 0.23 | 0.18
RankSHAP | 0.59 | 0.41 | 0.45 | 0.38 | 0.39 | 0.29
T5 | Random | 0.05 | 0.03 | -0.11 | 0.01 | -0.04 | 0.00
EXS | 0.36 | 0.21 | 0.30 | 0.17 | 0.11 | 0.06
RankLIME | 0.35 | 0.24 | 0.35 | 0.20 | 0.27 | 0.17
RankSHAP | 0.56 | 0.39 | 0.43 | 0.37 | 0.36 | 0.29

##  4 Computational and Human Evaluation Studies

In this section, we first describe our experimental design. Next we evaluate
the performance of our attribution algorithm based on various metrics. Lastly
we conduct a user study to investigate how well different attribution
algorithms align with human intuition for the ranking task.

###  4.1 Experimental Settings

Dataset:

We test our hypothesis on the Microsoft Machine Reading Comprehension (MS
Marco, (msm, 2016)) dataset. MS MARCO is a large-scale dataset aggregated from
anonymized Bing search queries containing $>8M$ passages from diverse text
sources. Within MS MARCO, we focus on the passage re-ranking task, which gives
a query and a set of relevant documents, retrieved using some rudimentary
algorithm. The goal is to re-rank this set of documents to enhance the
relevance of the top results for a given query. The average length of a
passage in the MS MARCO dataset is 1131 words.

Ranking models:

While the documents in our approach can be phrased in any representation, we
choose bag of words for our experiments based on its widespread prevalence.
For our ranking models, we use the classic BM25 ranker (Robertson & Spärck
Jones, 1994) as well as the BERT (Devlin et al., 2019) and T5 (Raffel et al.,
2020) neural rankers. We in particular use the MS MARCO fine-tuned version
made available by Nogueira & Cho (2019) for the neural ranking models. All the
above ranking approaches are score based rankers, i.e they return a relevance
score corresponding to each document. This makes relevance score estimation
straightforward for NDCG computation, as discussed previously, in Section 2.

Task:

To suit our case, we sample 200 queries from the test set of the MS MARCO
passage re-ranking dataset, and retrieve the 100 documents most relevant to
each query. We sample a variety of queries, including both factoid and non-
factoid queries, as well as queries with high and low relevant document
presence in the corpus. We use these query-document sets as input to the
ranking models described before and obtain an ordering of the documents with
respect to the query. We then generate feature attributions for these query
and ordered-document sets, generating attributions for top-10, top-20 as well
as top-100 documents.

Competing Systems:

We generate feature attributions for each of these query and ordered-document
sets using EXS (Singh & Anand, 2019), RankLIME (Chowdhury et al., 2023) and
RankSHAP (our approach). We also add comparison with a randomly generated
attribution system. For EXS, we use the rank-based setting as it is directly
comparable to RankSHAP. For RankLIME, we generate attributions using the
NeuralNDCG loss function in the single perturbation setting. All three of EXS,
RankLIME and Kernel-RankSHAP are perturbation based algorithms. They are all
assigned 5000 neighbourhood samples per decision, to create an even playing
field between the competing systems. We only process the 7 most significant
features generated by each algorithm for evaluation, in order to stay within
human comprehension limits (Chowdhury et al., 2023).

Evaluation Metrics:

Apart from human evaluations, feature attributions can be evaluated using
either fidelity-based metrics or via reference-based metrics. Since ground
truth attributions are not available in this case, we stick to fidelity-based
metrics, which measure the degree to which an attribution algorithm is
successful in recreating the underlying model decision (Anand et al., 2022).
For this, we use the metrics Fidelity and its weighted version wFidelity.
Given a feature attribution and a set of documents, we reconstruct the
attribution ordering by sorting linearly-combined document features multiplied
with their attributions. For Fidelity, we report the simplified unweighted
Kendall’s Tau between the original ranking model’s ordering and the
attribution-reconstructed ordering. For wFidelity, we penalize each swap in
passages found in the reconstructed ordering by the difference in their
positions. More formally, let $o_{f_{R}}=f_{R}(\vec{x})$ be the ranking output
by $f_{R}$, and let $o_{\phi}$ be the ranking of documents by their total
importance according to some feature attribution method $\phi$; that is, each
document in $\vec{x}$ is given a score of
$\sum_{i=1}^{M}\phi_{R}(f_{R},\vec{x},i)\times x_{i}$. Given two rankings $a$
and $b$, let $r_{a}[i]$ be the rank of the $i$-th element in $a$. We set
$\tau(a,b)$ to be the Kendall’s tau distance between $a$ and $b$, i.e.

$\displaystyle\frac{2}{n(n-1)}\sum_{i<j}\text{sgn}(r_{a}[i]-r_{a}[j])\cdot\text{sgn}(r_{b}[i]-r_{b}[j])$

and let $\tau_{w}(a,b)$ be the weighted Kendall’s Tau distance:

$\displaystyle\frac{\sum_{i<j}w_{ij}\cdot(\text{sgn}(r_{a}[i]-r_{a}[j])\cdot\text{sgn}(r_{b}[i]-r_{b}[j]))}{\sum_{i<j}w_{ij}}$

Here, $w_{ij}$ is $|r_{a}[i]-r_{a}[j]|$. Fidelity is then simply defined as
$\tau(o_{\phi},o_{f_{R}})$, and wFidelity is set to be
$\tau_{w}(o_{\phi},o_{f_{R}})$. The range of both Fidelity and wFidelity is
$[-1,1]$.

###  4.2 Experimental Results

The results of the above experiments can be seen in Table 1. We discuss them
below.

Competing Systems:

On an average, RankSHAP outperformed the best competing system by $30.78\%$ on
Fidelity and $23.68\%$ on weighted Fidelity. This shows that RankSHAP is quite
successful at recovering the original ordering. All metrics, other than
Random, show positive correlation between the original and reconstructed
ranking.

Across Metrics:

RankSHAP (and all the competing algorithms) perform worse on weighted-Fidelity
as compared to Fidelity, although both metrics have equal range. RankSHAP on
an average performs $17\%$ worse on weighted-Fidelity as compared to Fidelity,
showing the presence of inversions between documents that have significant
rank differences.

Number of documents:

The value of both metrics fall as the number of documents in the ordered-set
we generate attribution for increases. There is $20\%$ average drop in
performance between attributions for $10$ documents and attributions for $20$
documents and $14.6\%$ average Fidelity drop between $20\%$ documents and a
$100$ documents for RankSHAP, hinting towards a performance drop proportional
to the logarithm of documents. We witness that while EXS and RankLIME have
comparable Fidelity on generating attributions with $10$ and $20$ documents,
EXS attributions perform considerably worse in the $100$ document setting.

Ranking Model Complexity:

BM25 is a heuriristic-based ordering model while BERT and T5 are complex
neural ranking models. BERT and T5 show similar performance, which is on and
average $13.3\%$ and $14.7\%$ lower than BM25 on Fidelity for RankSHAP. As a
result, we see that Fidelity drops, with an increase in model complexity.

P1 So while raw honey and sugar both contain glucose and fructose, our raw
liquid gold is a nutrient dense food. Just as good grey sea salt is far
healthier that iodized table salt, raw honey beats the buzz off of regular
white sugar in both nutritional content and effects on the body.
---
P2 Agave nectar has a taste and appearance similar to honey, making it a
popular substitute for strict vegans and others who avoid honey. There are two
kinds of agave nectar: dark agave nectar and light. Though the two can be used
interchangeably, dark agave nectar has a stronger taste.
P3 Clover honey is named so because it’s produced from the nectar of clover
blossoms. It has a mild flavor and sweetness and is favored in the United
States for its abundance and versatility. When honey is labeled “pure,” it
means it has no additives such as sugar, corn syrup, or flavorings.
P4 Energy Source. According to the USDA, honey contains about 64 calories per
tablespoon. Therefore, it is used by many people as a source of energy. On the
other hand, one tablespoon of sugar will give you about 15 calories.
P5 1 Honey crystallizes because it is a supersaturated solution. 2 This
supersaturated state occurs because there is so much sugar in honey (more than
70%) relative to the water content (often less than 20%).
Actual Query: Is honey as a substitute for sugar healthier?

Figure 1: A user study instance, depicting 5 passages and a corresponding
feature attribution. Words that influenced the ranking decision positively are
shown in green, whereas those that influenced negatively were shown in red.
Participants were asked to re-rank documents based on the feature attribution
to try and estimate the query.

###  4.3 User Study

The goal behind developing a feature attribution method like RankSHAP, is to
aid a human user in understanding, why a complex model decided to order items
in a certain way. It is necessary to evaluate how well RankSHAP and other
ranking attribution methods fare amongst human subjects. We recruit $30$
participants for this study. Participants are shown sets of $5$ short passages
and a corresponding feature attribution generated by one of the competing
systems, displayed as a bar chart. An example set from the user study is shown
in Figure 1. Based on the feature attribution shown, the participants are
asked to (i) Re-order the $5$ documents from most important to least, and (ii)
Estimate what the query might have been, based on their intuition. The goal of
the first part is to figure which attribution algorithm is able to regenerate
the original ordering better and the goal for the second part is to see if the
attributions are indeed explaining the task in hand or something else.

We sampled query-document sets from the MS Marco dataset and BERT was used to
rank the documents. Each participant was shown $10$ such passage sets, with
attributions randomly selected from one of Random/EXS/RankLIME and RankSHAP.
At the end, the Kendall’s Tau score (Fidelity) between the original-ordering
and the predicted-ordering was used to identify the algorithm better in
reconstructing the underlying model decision. A text based semantic similarity
measuring LLM was used to measure the quality of the queries estimated by the
participants. The similarity metric returned scores in [0,1], 1 being when the
estimated query was most similar to the original query. The result of this
study can be witnessed in Table 2.

Table 2: Table demonstrating the user study results for (i)the passage reordering task, measured using Kendall’s Tau (Fidelity) (ii) the query estimation task, measured using GPT-4 similarity ($\mu$) | Random | EXS | RankLIME | RankSHAP
---|---|---|---|---
Q1 ($\tau$) | 0.23 | 0.43 | 0.47 | 0.56
Q2 ($\mu$) | 0.3 | 0.48 | 0.52 | 0.69

Results. At first, we observe that the randomly generated feature attributions
show a concordance of 0.23 on Task 1, which is higher than their metric
evaluated value (Table 1, Top-10 Fidelity - Random). This might imply that
participants also use pre-conceived notions on the topics, while participating
in the study. There was an inter-annotator agreement of 0.65 in Task 1.
RankSHAP fared $19.1\%$ better than the strongest baseline, following similar
trends as in Table 1. In the second task, certain query estimates by all
algorithms, were off-topic. We also witnessed a considerable variance in the
queries estimated by different participants for the very same set of documents
and feature attribution. However, amongst the compared systems, participants
when shown RankSHAP attributions fared at least $30.9\%$ better than when they
were shown other attributions. This shows that using RankSHAP, for ranking
feature attributions in practical use cases, would considerably aid human
understanding.

Table 3: Analyzing which feature attribution algorithms satisfy the ranking attributions axioms with NDCG value function, introduced in Section 3. Algorithm | Rank-Efficiency | Rank-Missingness | Rank-Symmetry | Rank-Monotonicity
---|---|---|---|---
EXS | ✗ | ✓ | ✗ | ✗
RankLIME | ✗ | ✓ | ✓ | ✗
RankSHAP | ✓ | ✓ | ✓ | ✓

##  5 Discussion

We first discuss alternate value functions for ranking attribution tasks.
Next, we deep dive in popular existing ranking attribution methods and conduct
a theoretical analysis to see if they satisfy the Shapley axioms for alterate
value functions.

###  5.1 Alternate Value Functions

The following functions could have been used in place of NDCG to measure the
effectiveness of an ordered-list.

Precision:

Measures proportion of retrieved documents that are more relevant than some
threshold.

Recall:

Measures proportion of documents, more relevant than some threshold, that have
been successfully retrieved.

F1-score:

The harmonic mean of Precision and Recall.

Top-k results:

Measures how often a desired correct results appear in the top $k$ results.

Mean Reciprocal Rate (MRR):

Reports the rank of the first relevant passage in ordered set.

Cumulative Gain (CG):

Reports the sum of all relevant documents in the list, without any
discounting.

However, NDCG is a better choice for the value function because of numerous
reasons. NDCG discounts the relevance of a document based on its position
logarithmically - which is known to mimic the human attention span as they
scan through a list. Unlike metrics that depend on binary relevance
(precision, recall, F1), NDCG accommodates a smooth relevance function. NDCG
is also more effective in measuring the quality of a large ranked list.
Lastly, a large number of product search on e-commerce and recommendation
systems (which we might be trying to explain) already use NDCG to evaluate the
quality of their orderings.

###  5.2 Axiomatic Analysis of EXS and RankLIME

Neither EXS nor RankLIME offer an axiomatic analysis, arguing that their
methods satisfy natural desiderata. In this subsection, we analyse if the
simplified algorithms behind EXS and RankLIME, satisfy the axioms defined
before. We also explore alternate value functions, that would potentially
allow these algorithms to follow the Shapley framework. We display the results
of this analysis in Table 3.

####  5.2.1 EXS

In EXS, Singh & Anand first use LIME (Ribeiro et al., 2016) to generate
feature attributions for each query-document pair individually. In order to
compute LIME attributions, they perturb documents and assign a binary
relevance label to each perturbation, assigning a positive label if the rank
of the perturbed document is higher than $k$. The authors then add the LIME
attributions of top-$k$ documents to assign a feature attribution to the
entire query-document set.

In our axiomatic analysis, we assume an implementation of EXS using KernelSHAP
(Lundberg & Lee, 2017) instead of LIME, as KernelSHAP is known to satisfy the
axioms of Efficiency, Missingness, Symmetry and Monotonicity for
regression/classification.

Rank-Efficiency :

The EXS algorithm described above does not satisfy the NDCG based Rank-
Efficiency axiom defined in Section 3 as-is. For a simple example, assume we
have an ordering with 5 documents with relevance scores $[3,2,3,1,2]$
respectively and the relevance limit $k$ of EXS set to $3$. The EXS
attributions would not align with the NDCG value function as, EXS would be
limited to the contribution of top-3 documents, would ignore the discount
factor based on the position, and would not be normalized to belong to a range
(like [0,1]). However, if in our Shapley axioms, instead of
$\textit{NDCG}_{n}$, we use a value function that sums up the contribution of
only the top-$k$ documents, without applying any discount, like Cumulative
Gain at $k$($\textit{CG}_{k}$), EXS would inherently satisfy the Rank-
Efficiency axiom.

Rank-Missingness :

Theoretically, EXS does satisfy the Rank-Missingness axiom, as is. We know
that EXS is a sum of $k$ KernelSHAP attribution and we know that KernelSHAP
satisfies the missingness axiom. Using the additivity axiom (another axiom
that KernelSHAP satisfies), we can conclude that EXS satisfies Missingness.
This is because the aggregate effect of non-contributing features across the
top-$k$ documents would also be non-contributory to the overall ranking
explanation provided by EXS.

Rank-Symmetry :

EXS does not satisfy rank-symmetry with the NDCG value function. Let us assume
two documents $d_{i}$ and $d_{j}$ and two feature $i$ and $j$. The attribution
of $i$ in $d_{i}$ is the same as the attribution of $j$ in $d_{j}$. Both
features do not have any attribution in the other document. Now assume we
assign $k=1$ for EXS. This will lead to computation of
$\textit{NDCG}_{2}([d_{i},d_{j}])$ for $\phi_{i}$ and $NDCG_{2}([d_{j},d_{i}]$
for $\phi_{j}$. These two will not be equal due to the potentially different
relevance scores of $d_{i}$ and $d_{j}$. Precision$@k$ is an alternate value
function for which EXS will satisfy Rank-Symmetry.

Rank-Monotonicity :

EXS does not satisfy monotonicity, as increases in a feature’s values does not
directly correlate with an increase in NDCG scores. A value function like
$CG_{k}$ would make EXS satisfy the Rank-Monotonicity axiom.

Singh & Anand propose other methods to compute the binary label of the
perturbed document as well, which we discuss in Appendix.

####  5.2.2 Rank-LIME

Chowdhury et al. (2023) extend LIME in a listwise manner, replacing the
$L_{2}$ loss function in LIME with a differentiable ordering distance metric.
This enables them to directly compute LIME feature attributions from ranked
lists. Two of the loss functions used by them, which lead to the best results,
are ApproxNDCG (Qin et al., 2010) and NeuralNDCG (Pobrotyn & Bialobrzeski,
2021), both differentiable alternatives to NDCG, with all relevance scores set
to 1. For the purpose of a theoretical analysis, we assume that both these
metrics successfully replicate the true NDCG value. To simplify the model for
analysis, we view each document independently, and frame the prediction of
their discounted relevance score ($\frac{1}{\log p}$ for the document with
rank $p$) as a regression task. We then add up the attributions over all
documents. We assume KernelSHAP has been used to compute each document’s
attribution since it satisfies desirable properties. Below we analyze if this
simplified version of RankLIME satisfies the ranking axioms.

Rank-Efficiency :

RankLIME in its simplified setting does not always satisfy Rank-efficiency
with an NDCG value function. It is also difficult to find an alternate value
function that satisfied Rank-Efficiency in Rank-LIME.

Rank-Missingness :

If a feature does not affect the ranking of any document within the list, and
consequently does not influence the NDCG score, KernelSHAP would attribute
zero importance to this feature. As a result, Rank-LIME satisfies Rank-
Missingness with the NDCG value function.

Rank-Symmetry :

In this simplified setting, Rank-LIME would satisfy symmetry. The aggregated
attribution would be the sum of multiple KernelSHAP attributions which would
individually satisfy summetry.

Rank-Monotonicity :

RankLIME does not satisfy monotonicity, as increases in a feature’s values
does not necessarily correlate with an increase in NDCG scores. A value
function like $CG_{k}$ would make RankLIME satisfy the Rank-Monotonicity
axiom.

##  6 Related Work

There have been a number of works in the information retrieval community
towards various kinds of explanations for ranking models. Most relevant to us
are EXS (Singh & Anand, 2019) and RankLIME (Chowdhury et al., 2023), model
agnostic, feature attribution generating methods for ranked lists, discussed
in detail in the previous sections. Fernando et al. (2019) extend the DeepSHAP
algorithm (Lundberg & Lee, 2017; Shrikumar et al., 2017) from
classification/regression to propose a set of model-intrinsic feature
attributions for neural ranking models. They observe significant variance
between feature attributions generated by LIME based methods as compared to
their approach on the very same feature space. Some authors have compute
attributions for ranking, but on a non-word feature space. Völske et al.
(2021) attempt to explain the performance of neural ranking models using well
understood axioms from statistical IR like term frequency and document length
(Fang et al., 2004), lower bounding term frequency (Lv & Zhai, 2011), query
aspects (Wu & Fang, 2012) etc. Others (Singh et al., 2020) introduce the
axioms of validity and completeness as a means to measure “ goodness” or
ranking explanations and propose an algorithm to suggest a small set of
features which are sufficient to explain a ranking decision. LRIME (Verma &
Ganguly, 2019) propose a set of heuristics to make regression/classification
based feature attribution methods work better for information retrieval, like
choosing a diverse set of perturbation samples and parameters. However, none
of these works have theoretical grounding or proper evaluation metrics, due to
the complexity involved in evaluating feature attribution methods. While the
above discussed methods to find attributions for learning to rank models were
post-hoc and model-agnostic in nature, there have also been efforts to
generate intrinsically interpretable learning to rank models. Ranking GAMs
(Zhuang et al., 2020) in an inherently interpretable structure that can be
distilled into a set of compact piece-wise linear functions with little
accuracy loss. Anand et al. (2022) summarize explanation methods proposed for
the ranking task comprehensively, also delving into free-text explanations
(Rahimi et al., 2021), adversarial example based explanations (Raval & Verma,
2020; Wu et al., 2023), probing the representation space (Choi et al., 2022),
rationale based explanations (Zhang et al., 2021; Wojtas & Chen, 2020) etc.

##  7 Conclusion and Future Work

While there has been considerable work towards axiomatically valid feature
attributions in regression/classification, there seemed to exist no standard
attribution when it came to ranking. In this work we proposed RankSHAP: a
valid feature attribution framework for ranking problems. We first introduced
axioms that we believe a valid feature attribution algorithms must satisfy. We
then provided with a solution that satisfied the axioms. We also presented a
computationally efficient approximation of RankSHAP and conducted a user study
to show the degree to which our proposed solution aligned with human
intuition. We also did an analysis of existing ranking feature attribution
algorithms EXS and RankLIME and discussed if they followed the axioms we
defined initially. We also proposed modifications to our value function, in
order to make existing algorithms Shapley compliant. All the properties above,
show that our solution makes a good candidate to be used as a gold standard
for generating ranking feature attributions.

##  8 Ethics

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none
which we feel must be specifically highlighted here.

## References

  * msm (2016) MS MARCO: Microsoft Machine Reading Comprehension.  https://microsoft.github.io/msmarco/, 2016.  Accessed: [3.5.2023].
  * Anand et al. (2022) Anand, A., Lyu, L., Idahl, M., Wang, Y., Wallat, J., and Zhang, Z.  Explainable information retrieval: A survey.  _arXiv preprint arXiv:2211.02405_ , 2022.
  * Choi et al. (2022) Choi, J., Jung, E., Lim, S., and Rhee, W.  Finding inverse document frequency information in bert.  _arXiv preprint arXiv:2202.12191_ , 2022.
  * Chowdhury et al. (2023) Chowdhury, T., Rahimi, R., and Allan, J.  Rank-LIME: Local model-agnostic feature attribution for learning to rank.  In _Proceedings of the 46th International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pp. 33–37, 2023\.
  * Datta et al. (2016) Datta, A., Sen, S., and Zick, Y.  Algorithmic transparency via quantitative input influence: Theory and experiments with learning systems.  In _Proceedings of the 37th IEEE Conference on Security and Privacy (Oakland)_ , pp. 598–617, 2016.
  * Devlin et al. (2019) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Annual Conference of the North American Chapter of the Association for Computational Linguistics (NAACL)_ , pp. 4171–4186, 2019.
  * Fang et al. (2004) Fang, H., Tao, T., and Zhai, C.  A formal study of information retrieval heuristics.  In _Proceedings of the 27th International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pp. 49–56, 2004\.
  * Fernando et al. (2019) Fernando, Z. T., Singh, J., and Anand, A.  A study on the interpretability of neural retrieval models using DeepSHAP.  In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pp. 1005–1008, 2019.
  * Järvelin & Kekäläinen (2002) Järvelin, K. and Kekäläinen, J.  Cumulated gain-based evaluation of ir techniques.  _ACM Transactions on Information Systems_ , 20(4):422–446, 2002.
  * Lundberg & Lee (2017) Lundberg, S. M. and Lee, S.-I.  A unified approach to interpreting model predictions.  In _Proceedings of the 30th Annual Conference on Neural Information Processing Systems (NeurIPS)_ , pp. 4768–4777, 2017.
  * Lv & Zhai (2011) Lv, Y. and Zhai, C.  Lower-bounding term frequency normalization.  In _Proceedings of the 20th ACM international conference on Information and Knowledge Management (CIKM)_ , pp. 7–16, 2011.
  * Nogueira & Cho (2019) Nogueira, R. and Cho, K.  Passage re-ranking with bert.  _arXiv preprint arXiv:1901.04085_ , 2019.
  * Pobrotyn & Bialobrzeski (2021) Pobrotyn, P. and Bialobrzeski, R.  Neuralndcg: Direct optimisation of a ranking metric via differentiable relaxation of sorting.  _arXiv preprint arXiv:2102.07831_ , 2021.
  * Qin et al. (2010) Qin, T., Liu, T.-Y., and Li, H.  A general approximation framework for direct optimization of information retrieval measures.  _Information retrieval_ , 13(4):375–397, 2010\.
  * Raffel et al. (2020) Raffel, C., Shazeer, N., Roberts, A., Lee, K., Narang, S., Matena, M., Zhou, Y., Li, W., and Liu, P. J.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _Journal of Machine Learning Research_ , 21(140):1–67, 2020.
  * Rahimi et al. (2021) Rahimi, R., Kim, Y., Zamani, H., and Allan, J.  Explaining documents’ relevance to search queries.  _arXiv preprint arXiv:2111.01314_ , 2021.
  * Raval & Verma (2020) Raval, N. and Verma, M.  One word at a time: adversarial attacks on retrieval models.  _arXiv preprint arXiv:2008.02197_ , 2020.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  “why should i trust you?” explaining the predictions of any classifier.  In _Proceedings of the 22nd Special Interest Group on Knowledge Discovery in Data (SIGKDD)_ , pp. 1135–1144, 2016.
  * Robertson & Spärck Jones (1994) Robertson, S. E. and Spärck Jones, K.  Simple, proven approaches to text retrieval.  Technical report, University of Cambridge, Computer Laboratory, 1994.
  * Shapley (1953) Shapley, L. S.  A value for n-person games.  _Contributions to the Theory of Games_ , 2(28):307–317, 1953.
  * Shrikumar et al. (2017) Shrikumar, A., Greenside, P., and Kundaje, A.  Learning important features through propagating activation differences.  In _Proceedings of the 34th International Conference on Machine Learning (ICML)_ , pp. 3145–3153, 2017.
  * Singh & Anand (2019) Singh, J. and Anand, A.  Exs: Explainable search using local model agnostic interpretability.  In _Proceedings of the 12th ACM International Conference on Web Search and Data Mining (WSDM)_ , pp. 770–773, 2019.
  * Singh et al. (2020) Singh, J., Khosla, M., and Anand, A.  Valid explanations for learning to rank models.  _arXiv preprint arXiv:2004.13972_ , 2020.
  * Verma & Ganguly (2019) Verma, M. and Ganguly, D.  Lirme: locally interpretable ranking model explanation.  In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR)_ , pp. 1281–1284, 2019.
  * Völske et al. (2021) Völske, M., Bondarenko, A., Fröbe, M., Stein, B., Singh, J., Hagen, M., and Anand, A.  Towards axiomatic explanations for neural ranking models.  In _Proceedings of the 2021 ACM SIGIR International Conference on the Theory of Information Retrieval (ICTIR)_ , pp. 13–22, 2021.
  * Wojtas & Chen (2020) Wojtas, M. and Chen, K.  Feature importance ranking for deep learning.  _Advances in Neural Information Processing Systems_ , 33:5105–5114, 2020.
  * Wu et al. (2023) Wu, C., Zhang, R., Guo, J., De Rijke, M., Fan, Y., and Cheng, X.  Prada: practical black-box adversarial attacks against neural ranking models.  _ACM Transactions on Information Systems_ , 41(4):1–27, 2023.
  * Wu & Fang (2012) Wu, H. and Fang, H.  Relation based term weighting regularization.  In _Proceedings of the 34th European Conference on Information Retrieval (ECIR)_ , pp. 109–120, 2012.
  * Young (1985) Young, H. P.  Monotonic solutions of cooperative games.  _International Journal of Game Theory_ , 14(2):65–72, 1985.
  * Zhang et al. (2021) Zhang, Z., Rudra, K., and Anand, A.  Explain and predict, and then predict again.  In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining (WSDM)_ , pp. 418–426, 2021.
  * Zhuang et al. (2020) Zhuang, H., Wang, X., Bendersky, M., Grushetsky, A., Wu, Y., Mitrichev, P., Sterling, E., Bell, N., Ravina, W., and Qian, H.  Interpretable learning-to-rank with generalized additive models.  _arXiv preprint arXiv:2005.02553_ , 2020.
  * Zoeter et al. (2008) Zoeter, O., Taylor, M., Snelson, E., Guiver, J., Craswell, N., and Szummer, M.  A decision theoretic framework for ranking using implicit feedback.  In _SIGIR 2008 workshop on Learning to Rank for Information Retrieval_ , 2008.


# A Parallel Monte-Carlo Tree Search-Based Metaheuristic For Optimal Fleet
Composition Considering Vehicle Routing Using Branch & Bound

T.M.J.T. Baltussen1, M. Goutham1, M. Menon2, S.G. Garrow2, M. Santillo2, S.
Stockar1

###### Abstract

In this paper, a Monte-Carlo Tree Search (MCTS)-based metaheuristic is
developed that guides a Branch & Bound (B&B) algorithm to find the globally
optimal solution to the Fleet Size and Mix Vehicle Routing Problem with Time
Windows (FSMVRPTW). The metaheuristic and exact algorithms are implemented in
a parallel hybrid optimization algorithm where the metaheuristic rapidly finds
feasible solutions that provide candidate upper bounds for the B&B algorithm
which runs simultaneously. The MCTS additionally provides a candidate fleet
composition to initiate the B&B search. Experiments show that the proposed
approach results in significant improvements in computation time and
convergence to the optimal solution.

Keywords: Fleet composition, Vehicle Routing, Branch & Bound, Monte-Carlo Tree
Search, Metaheuristic

11footnotetext: The Ohio State University22footnotetext: Ford Motor
Company††footnotetext: This work has been submitted to the IEEE for possible
publication. Copyright may be transferred without notice, after which this
version may no longer be accessible. Manuscript submitted 1 Feb 2023 to the
IEEE Intelligent Vehicles Symposium 2023

##  I Introduction

In the industrial sector, Reconfigurable Manufacturing Systems (RMS) are
increasingly being adopted because of their ability to scale and diversify
production by supporting the adaptability of process controls, functions and
operations [1, 2]. A key enabler is the added production flexibility provided
by the adoption of fleets of autonomous mobile robots (AMRs) that move
material within an RMS plant [3]. In particular, multi-load AMRs enhance
efficiency by picking up and dropping off multiple items in a single mission
[4]. The design of such a fleet is a strategic problem and involves
considerable capital investment. Therefore, all costs related to the
acquisition and operation should be considered [5]. Despite this, the combined
fleet composition and routing problem has received insufficient attention [5].
In this paper, the Vehicle Routing Problem with Time Windows (VRPTW) and
capacity constraints on the cargo mass, volume and vehicle range is used to
obtain the operational cost. The combined VRPTW with the heterogeneous fleet
composition problem, is called the Fleet Size and Mix Vehicle Routing Problem
with Time Windows (FSMVRPTW). This problem accommodates a heterogeneous fleet
and takes vehicle costs into consideration in addition to operational costs
[5].

Fleet composition optimization problems are typically posed as a capacitated
VRPTW where the fleet size can be varied [6]. Exact algorithms that guarantee
optimality for this $\mathcal{NP}$-hard combinatorial optimization problem
structure the optimization problem as a tree exploration problem and are
solved using Branch & Bound (B&B) methods [6]. Due to the $\mathcal{NP}$-hard
nature of the problem, the application of exact algorithms is restricted to
small problem instances [7]. Real-life Vehicle Routing Problem applications
are considerably larger in scale [7] and finding the optimal solution to such
a problem is computationally expensive. Therefore, most VRPTWs are solved
using metaheuristic methods due to their ability to find near-optimal solution
in a limited time[7, 6]. However, such approximate methods do not provide
guarantees on the optimality of the solution [6, 8].

Hybrid optimization methods can improve performance and efficiency of the
optimizer by combining the strengths of metaheuristics and exact algorithms.
Successful metaheuristics provide a balance between exploration and
exploitation of the search space [9]. As such, Monte-Carlo Tree Search (MCTS)
is a reinforcement learning algorithm that balances this exploration and
exploitation and it is well suited to large-scale combinatorial optimization
problems [6, 10, 11]. In fact, MCTS has already been used in literature as a
metaheuristic that guides a CPLEX solver toward the optimal solution [12].
Moreover, it is frequently hybridized with other optimization algorithms [11].
MCTS has been found to obtain state-of-the-art results in resource allocation
problems (RAP) [13] and in single vehicle instances of the VRPTW, called
travelling salesperson problems with time windows (TSPTW) [14]. It has also
been used to solve VRP problems with variable fleet sizes[13, 14, 15].
However, MCTS has not yet been used to solve FSMVRPTWs that permit different
types of vehicles.

The first contribution of this paper is the development of an exact Fleet
Composition Optimizer (FCO) algorithm for the FSMVRPTW. This algorithm employs
a divide and conquer approach where the VRPTW is partitioned to a resource
allocation problem (RAP) that first assigns tasks to each robot using a
parallel B&B algorithm, and then finds the optimal sequence in which the
assigned tasks are completed by solving a nested TSPTW, using another B&B
algorithm. The second contribution is a hybrid MCTS-based metaheuristic, that
uses the Upper Confidence bounds applied to Trees (UCT) algorithm [16]. This
UCT-based metaheuristic (UCT-MH) uses the UCT algorithm in fleet composition
levels to guide its search and solves the nested TSPTW using a B&B algorithm.
The third novelty presented in this paper is the hybrid optimization framework
where the UCT-MH guides the FCO to find the optimal solution to the FSMVRPTW.
When possible, the FCO is initialized with a fleet composition identified by
the rapid search space exploration enabled by the UCT-MH. Additionally, the
best solutions found by the UCT-MH update the upper bound used by the FCO,
which allows sub-optimal solutions to be pruned earlier. The performance of
the proposed method is verified on various real-life case studies. Results
show a significant reduction in computation time when the FCO algorithm is
guided by the proposed UCT-MH, especially for large problems sizes.

##  II Problem Formulation & Methodology

Consider a manufacturing plant with a known layout that comprises various
spatial constraints, and a set of material handling tasks $\mathcal{T}$. Each
task involves picking-up certain cargo items at inventory locations and
dropping them off at their designated drop-off locations within defined time
windows. The objective of the optimization is to find the optimal fleet of
multi-load capacitated AMRs that completes all the defined tasks $\mathcal{T}$
while minimizing fixed and operational costs.

Let the set $\mathcal{H}:=\\{1,2,...,h\\}$ identify $h$ different AMR types
available, each with specific traveling speeds, energy efficiency, cargo
capacity, driving range, charge-time etc. Let $k_{i}\leq
k_{i}^{max}:i\in\mathcal{H}$ denote the number of each type of AMR that forms
a fleet so that any fleet composition can be fully defined by a vector
$\mathbf{k}\in\mathbb{N}_{0}^{h}$. This fleet is associated with a fixed cost
$J^{f}(\mathbf{k})$ composed of purchase costs, depreciation, etc., that can
be captured by $J^{f}(\mathbf{k})=\mathbf{c}^{\top}\mathbf{k}$ for some
$\mathbf{c}\in\mathbb{R}^{h}$. For completing all the tasks in $\mathcal{T}$,
the operational cost $J^{o}(\mathbf{k})$ can be any combination of relevant
metrics to be minimized such as energy, slack time, number of turns, asset
depreciation, etc. [17, 18, 19]. The total cost to be minimized is:

$\displaystyle\min_{\mathbf{k}}J=\mathbf{c}^{\top}\mathbf{k}+J^{o}(\mathbf{k})$
(1)

The fleet operational cost $J^{o}(\mathbf{k})$ is posed as an RAP that finds
the optimal partition of tasks to be assigned to AMRs that minimizes total
operational cost. If the total number of robots in the heterogeneous fleet
$\mathbf{k}$ is given by $m=\sum_{i=1}^{h}k_{i}$, every robot in this fleet
can be identified by $r\in\mathcal{R}_{k}:=\\{1,2,...,m\\}$. Let the set
$\mathcal{T}_{r}\subseteq\mathcal{T}$ denote the tasks assigned to robot $r$
by the partitioning of $\mathcal{T}$, denoted by
$\mathfrak{T}:=\\{\mathcal{T}_{r}:r\in\mathcal{R}_{k}\\}$, meaning
$\bigcup_{r\in\mathcal{R}_{k}}\mathcal{T}_{r}=\mathcal{T}$ and $\forall
r,s\in\mathcal{R}_{k}:r\neq
s,\mathcal{T}_{r}\bigcap\mathcal{T}_{s}=\varnothing$. The optimal partition of
task set $\mathcal{T}$ minimizes $J^{o}(\mathbf{k})$ in Eq (2).

$\displaystyle
J^{o}(\mathbf{k})=\min_{\mathfrak{T}}\sum_{r\in\mathcal{R}_{k}}J^{r}(r,\mathcal{T}_{r})$
(2)

The minimum operational cost $J^{r}(r,\mathcal{T}_{r})$ for each robot in
fleet $\mathbf{k}$ is dependent on the AMR type, and is also affected by the
sequence with which task locations are visited, as it is possible for the
robot to pickup multiple items before dropping them off so long as each pickup
is visited before the corresponding drop-off. The objective function and
constraints that yield $J^{r}(r,\mathcal{T}_{r})$ are defined in Eq. (3-15).

Let robot $r$ of the fleet be assigned $n_{r}=|\mathcal{T}_{r}|$ tasks. The
set of pickup and drop-off locations are defined as
$\mathcal{V_{P}}:=\\{1,2,...,n_{r}\\}$ and
$\mathcal{V_{D}}:=\\{n_{r}+1,n_{r}+2,...,2n_{r}\\}$ respectively, so that an
item picked up at location $i$ must be dropped off at location $n_{r}+i$. The
origin and final destination locations of the robot are identified by
$\\{0,2n_{r}+1\\}$. Let
$\mathcal{V}:=\\{\mathcal{V_{P}}\cup\mathcal{V_{D}}\cup\\{0,2n_{r}+1\\}\\}$ be
the set of all locations in a graph representation
$\mathcal{G}:=(\mathcal{V},\mathcal{A})$ where
$\mathcal{A}:=\\{(i,j)\in\mathcal{V}\times\mathcal{V}\\}$ is the arc set.
Between every pair of nodes $(i,j)\in\mathcal{A}$, the operational costs
$D_{ij}\in\mathbb{R}^{+}$, energy consumed $\delta e_{ij}\in\mathbb{R}^{+}$
and travel time $\delta t_{ij}\in\mathbb{R}^{+}$ are pre-computed before
initializing the optimization by solving a path planning problem between every
two locations $i,j\in\mathcal{V}:(i,j)\in\mathcal{A}$.

$\displaystyle J^{r}(r,\mathcal{T}_{r})=\min_{\begin{subarray}{c}x_{ij}\\\
\forall ij\in\mathcal{A}\end{subarray}}\leavevmode\nobreak\
\sum_{i=0}^{2n_{r}}\leavevmode\nobreak\ \sum_{j=1}^{2n_{r}+1}D_{ij}x_{ij}$ (3)
s.t. $\displaystyle x_{ij}\in\\{0,1\\}\leavevmode\nobreak\
\forall(i,j)\in\mathcal{A}$ (4) $\displaystyle\sum_{j=1}^{n_{r}}x_{0j}=1$ (5)
$\displaystyle\sum_{i=0}^{2n_{r}}$ $\displaystyle
x_{il}=\sum_{j=1}^{2n_{r}+1}x_{lj}=1,\quad\forall
l\in\\{\mathcal{V}\setminus\\{0,2n_{r}+1\\}\\}$ (6)
$\displaystyle\sum_{i=n_{r}+1}^{2n_{r}}x_{i,2n_{r}+1}=1$ (7) $\displaystyle
z_{j}=$ $\displaystyle\begin{cases}z_{i}-\delta e_{ij}\\\ \leavevmode\nobreak\
\;\text{if }x_{ij}=1\wedge z_{i}-\delta e_{ij}>0\leavevmode\nobreak\
\leavevmode\nobreak\ \forall(i,j)\in\mathcal{A}\\\ 1-\delta e_{0j}\\\
\leavevmode\nobreak\ \;\text{if }x_{ij}=1\wedge z_{i}-\delta e_{ij}\leq
0\leavevmode\nobreak\ \leavevmode\nobreak\ \forall(i,j)\in\mathcal{A}\\\
\end{cases}$ (8) $\displaystyle\quad z_{0}=$ $\displaystyle
1;\leavevmode\nobreak\ 0\leq z_{i}\leq 1\leavevmode\nobreak\ \forall
i\in\mathcal{V}$ (9) $\displaystyle T_{ij}=$ $\displaystyle\begin{cases}\delta
t_{ij}\\\ \quad\leavevmode\nobreak\ \text{if }x_{ij}=1\wedge z_{i}-\delta
e_{ij}>0\leavevmode\nobreak\ \leavevmode\nobreak\
\forall(i,j)\in\mathcal{A}\\\ \delta t_{0i}+(1-z_{i}-\delta
e_{i0})p^{-1}+\delta t_{0j}\\\ \quad\leavevmode\nobreak\ \text{if
}x_{ij}=1\wedge z_{i}-\delta e_{ij}\leq 0\leavevmode\nobreak\
\leavevmode\nobreak\ \forall(i,j)\in\mathcal{A}\end{cases}$ (10)
$\displaystyle x_{ij}=1\rightarrow t_{i}+s_{i}+T_{ij}\leq
t_{j}\leavevmode\nobreak\ \forall(i,j)\in\mathcal{A}$ (11) $\displaystyle\quad
t_{i}+s_{i}+T_{i,n+i}\leq t_{n+i}\leavevmode\nobreak\ \forall i\in\mathcal{V}$
(12) $\displaystyle\quad e_{i}\leq t_{i}\leq l_{i}\leavevmode\nobreak\ \forall
i\in\mathcal{V}$ (13) $\displaystyle x_{ij}=1\rightarrow
y_{j}=y_{i}+q_{j}\leavevmode\nobreak\ \forall(i,j)\in\mathcal{A}$ (14)
$\displaystyle\quad y_{0}=0;\leavevmode\nobreak\ 0\leq y_{i}\leq
Q_{r}\leavevmode\nobreak\ \forall i\in\mathcal{V}$ (15)

The binary flow variable $x_{ij}=1$ signifies that the robot uses directed arc
$(i,j)\in\mathcal{A}$. Constraints related to the robot starting from the
depot $0$, visiting every location once and terminating the sequence at
$2n_{r}+1$ are enforced by Eq. (4-7).

The battery states of charge $z_{j}$ in Eq. (8-10) are updated as the robot
goes about its mission. Whenever the battery is depleted, the robot heads to
the depot where it is fully charged up with a constant recharging rate $p$.
The variable $T_{ij}$ in Eq. (10-12) updates the travel time between locations
$i$ and $j$ based on whether a recharge is required between the two locations.

Time variables $t_{i}$ in Eq. (11-13) denote the arrival time of the robot at
location $i\in\mathcal{V}$. Each location is associated with a time $s_{i}$
for material handling and a time window $[e_{i},l_{i}]$ which represents the
earliest and latest time at which material handling can start.

Cargo constraints are captured in Eq. (14, 15) where payload variables $y_{i}$
capture the cargo mass being carried by the robot as it leaves location
$i\in\mathcal{V}$. Each robot $r$ has a cargo capacity limitation of $Q_{r}$
and each location $i\in\mathcal{V}$ is associated with a cargo load
$q_{i}\in\mathbb{R}$ such that $q_{i}+q_{n+i}=0$. Volumetric constraints are
modelled similarly.

###  II-A Exact Algorithm: Fleet Composition Optimizer

The FCO systematically partitions the search space into subsets that are
arranged in a tree structure. The root of the tree is the original problem and
the leaves of the tree are its individual candidate solutions. Between the
root and the leaves are intermediate nodes that represent subproblems obtained
by recursively partitioning the original problem by a process called
branching. B&B algorithms are used to solve these sub-problems. The order
according to which these subproblems are examined is determined by a best-
first selection criteria, i.e. exploitation, that first explores the problem
with cheapest cost.

For minimization problems, the upper bound is the incumbent solution defined
as the cheapest candidate solution to the original problem found at the leaf
node. The upper bound is continuously updated as the tree is explored, and is
used to prune sub-optimal branches without recursively evaluating their
solutions up to the leaf node. Thus, as the algorithm searches from the root
to the leaves, branching is conducted only if the cost at the node is lower
than the incumbent solution, and branching can potentially find a better
solution than the incumbent solution. Following this process, the B&B
algorithm recursively decomposes the original problem until further branching
is futile when the solution cannot be improved, or until the original problem
has been solved when every feasible branch has been evaluated.

Figure 1: Parallel implementation of the RAP B&B algorithm for the Resource
Allocation Problem formulation.

The $\mathcal{NP}$-hard RAP problem described by Eq. 2 is solved by the FCO
using a B&B algorithm implemented in a parallel framework that uses $p$
processing cores, as shown in Fig. 1, where robots in the fleet are identified
by subscript $r\in\mathcal{R}:=\\{1,2,...,m\\}$. Thus, by splitting the
arborescence at some task assignment level and assigning the emanating sub-
trees to the available processors, several subproblems are explored
simultaneously. During each processor’s exploration, updated incumbent
solutions are instantaneously made available to every processor in an
asynchronous information sharing method using a shared work pool. For each
processor, this RAP B&B algorithm is implemented by a recursive function to
minimize memory and computational requirements as the tree is explored.
Further, since the computation time of B&B algorithms increases with the
number of feasible branches at each node, the fleet is initiated with a
smaller candidate fleet $\mathbf{f}^{1}\in\mathbb{N}_{0}^{h}$ than the maximal
fleet $\mathbf{k}^{max}\in\mathbb{N}_{0}^{h}$. After evaluating the total cost
of this candidate fleet, the number of robots is incrementally raised until
further increments do not reduce the total cost or additional robots remain
idle. For each fleet increment, only RAP subproblems that include at least one
of the newly added robots are evaluated, since other solutions are guaranteed
to have been evaluated already.

For $h$ different AMR types available, the fleet is initiated with a candidate
fleet $\mathbf{f}^{1}\leq\mathbf{k}^{max}$, which is chosen based on problem
parameters and prior experience so that feasible solutions exist. The RAP of
fleet $\mathbf{f}^{1}$ is then solved using the described parallel B&B
algorithm, and its minimum total cost $J^{1}$ is found, which utilizes robots
$\mathbf{k}^{1}\in\mathbb{N}_{0}^{h}:\mathbf{k}^{1}\leq\mathbf{f}^{1}$. In the
increment step, only robot types $i$ that satisfy
$\mathbf{k}^{1}_{i}=\mathbf{f}^{1}_{i}$ are incremented by 1 for the next
candidate fleet $\mathbf{f}^{2}$. These increments are conducted so long as
both $J^{i+1}\leq J^{i}$ and $\exists\leavevmode\nobreak\
i:\mathbf{k}^{1}_{i}=\mathbf{f}^{1}_{i}$. The optimal fleet that completes all
tasks while minimizing total cost is then $k^{i}$ when $J^{i+1}>J^{i}$ or when
$\nexists\leavevmode\nobreak\ i:\mathbf{k}^{i+1}_{i}=\mathbf{f}^{i+1}_{i}$,
i.e, the additional robots were idle.

At each instance that the RAP subproblem is solved at a node in the
arborescence shown in Fig. 1, the TSPTW problem defined by Eq. (3-15) is
solved to find the cost at that node. This TSPTW problem is solved using
recursive Algorithm 1 that employs another B&B to find the optimal sequence of
task completion for each robot. In summary, the FCO incrementally increases
the fleet size while minimizing total cost $J$ of Eq. 1, which includes
operational cost $J^{o}(\mathbf{k})$ of Eq. 2 found using the RAP B&B
algorithm, which includes the cost $J^{r}(r,\mathcal{T}_{r})$ of Eq. (3-15)
found using the TSPTW B&B algorithm.

Algorithm 1 Recursive TSPTW B&B
Cost = _B &B_(State, taskList, current Location)

1: Find feasible next locations based on completed pickups, time, cargo,
battery constraints

2: Sort feasible next locations by operational cost of branching to that
location (Best First Search)

3: for  _i_ in feasible locations do

4:  _branchCost_ = tourCost + operational cost(i)

5: if  _branchCost_ $\geq$ State.bestCost then

6: continue { skip to next _location i+_}

7: else if  _branchCost_ $<$ State.bestCost then

8:  _State+_ = Update State of Time, SOC, final position, remaining locations

9: if remainingLocations $>0$ then

10: Cost = _B &B_(State+, taskList, _location_(i))

11: else

12: State.bestCost = Cost

13: end if

14: end if

15: end for

###  II-B Metaheuristic: Monte-Carlo Tree Search

Each iteration of MCTS involves four steps [20]: Selection: at every node $v$
in the arborescence, the tree policy selects the next node $v^{\prime}$. This
node selection is initiated at the root node $v_{0}$ and is used for
navigation until the leaf node $v_{l}$ is reached. Expansion: at the leaf node
$v_{l}$, a random action is taken to expand the tree. Simulation: a Monte-
Carlo simulation is performed starting from the expansion node to complete the
solution. Backpropagation: the cost/reward of the expansion and simulation is
propagated back to the root node $v_{0}$. The Upper Confidence bounds applied
to Trees (UCT) algorithm was the first variant and formal introduction of MCTS
[16].

The proposed UCT-based metaheuristic (UCT-MH), is a hybrid MCTS algorithm
which is intended to guide the exact FCO algorithm to the optimal solution.
While in typical VRPTWs the fleet size is a free variable [6], the proposed
metaheuristic selects a fleet size $m$ and a composition $\mathbf{k}$ in the
fleet sizing and composition stages and tries to solve the VRPTW optimally,
fully utilizing that composition. By doing so, the algorithm finds an estimate
of the expected total cost associated with a particular fleet size and
composition. This estimate serves as a measure for the quality of that branch
and can be used by the MCTS to navigate the search. MCTS is most effective as
a heuristic at the early stages of the decision problem [12]. Moreover, for
smaller problem instances, B&B algorithms are often more suitable than MCTS
[14]. As such, the proposed hybrid MCTS algorithm is aimed to utilize the
strengths of the different algorithms and combine them into an effective
hybrid MCTS-based metaheuristic.

Figure 2: Overview of the multi-stage design problem, with the FSMVRPTW (red)
and the nested VRPTW (blue), and the proposed UCT-MH Algorithm.

Although MCTS was originally designed to solve Markov Decision Processes
(MDP), without loss of generality, MCTS can be used to solve a design problem
by formulating it as a deterministic MDP [11]. The optimization problem is
modeled as a 3-tuple $\langle S,A,g\rangle$, where $S$ is a set of states, $A$
is a set of actions and $g(s,a):S\times A\rightarrow[0,g_{max}]$ is a scalar
cost function for taking action $a$ at state $s$. The state $s(v)$ contains
the parameters that follow from the decisions up to node $v$. At the root node
$v_{0}$, the fleet size $m$ is determined by action $a_{0}$, where
$g_{1}(s_{0}(v_{0}),a_{0}):=0$, for the fleet cost is determined by its
composition. Subsequently, the fleet composition $\mathbf{k}$ is determined by
$a_{1}\in A_{1}(m)$, with fixed cost
$g_{2}(s_{1}(m),a_{1})=J^{f}(\mathbf{k})$. Fig. 2 provides a schematic
overview of the problem and the proposed metaheuristic.

At the fleet sizing and composition stages, the UCT-MH utilizes the UCB1 tree
policy [16] for the selection step at node $v$ of the decision tree:

$\text{UCB1}(v)=\operatorname*{arg\,max}_{v^{\prime}\in\text{children of
}v}\frac{Q(v^{\prime})}{N(v^{\prime})}+\sqrt{\frac{2\ln N(v)}{N(v^{\prime})}}$
(16)

Here, $Q(v^{\prime})$ is the total reward of all plays through child node
$v^{\prime}$, $N(v^{\prime})$ denotes the number of visits of child node
$v^{\prime}$, and $N(v)$ is the number of visits of the parent node $v$. The
policy function is dependent on the quality of the node being considered as
well as the number of evaluations of that node, balancing the exploration and
exploitation of the search space[20]. In order to apply the UCB1 policy and
have a proper balance between exploration and exploitation, the problem is
transformed such that the stage reward $R_{i}(v)\in[0,1]$ [16]:

$R_{i}(v^{\prime})=1-\frac{g_{i}(v^{\prime})}{g_{max}}$ (17)

where $R_{i}(v^{\prime})$ is the reward of the transition from state
$s_{i-1}(v)$ to state $s_{i}(v^{\prime})$ and $v^{\prime}\in\text{children of
}v$. It follows that $Q(v^{\prime})$ is the sum of all rewards of all
$N(v^{\prime})$ plays through node $v^{\prime}$ back to the root node $v_{0}$:

$Q(v^{\prime})=\sum_{i=1}^{N(v^{\prime})}R_{i}(v^{\prime})+R_{i}(v)+...+R_{i}(v_{0})$
(18)

Considering that the number of permutations of the RAP is exponential with the
number of tasks, it is deemed sufficient to determine the task assignment by a
random rollout $\left(\xi_{1},...,\xi_{n}\right)$. In order to prevent any
bias toward another fleet size, it is ensured that the full fleet size is
utilized, i.e. each AMR in the fleet will have at least one assignment. The
assigned tasks do not have any associated costs/rewards.

Since many of the TSPTW instances encountered are small problem instances, it
is advantageous to use the same recursive B&B algorithm for TSPTW as described
in Section II-A to find the optimal sequence in which the assigned tasks are
completed by each robot. Each TSPTW B&B is terminated after a one second time
cap since the metaheuristic is not aimed at local convergence. Considering the
best first order of exploration, this still finds reasonably good estimates
for the operational cost $\tilde{J}^{o}(k)$. The cost that is obtained through
the rollout of the RAP and the TSPTW, is backpropagated through the tree and
are assigned to $Q(v)$ at node $v$ that is associated with a particular fleet
size or composition. This is in turn used by the UCB1 policy function to
determine the decisions in the next iteration. As a result, at the root node,
the term $\frac{Q(v^{\prime})}{N(v^{\prime})}$ in (16) is proportional to the
total mean cost-to-go for a given fleet size or composition at node
$v^{\prime}$. As the total number of plays at the root node $N(v_{0})$ grows
to infinity, the UCB1 function converges to the expected value of the total
cost for a given fleet size.

###  II-C Hybrid Optimization: Guiding B&B with the UCT-MH

The hybrid optimization framework utilizes the search results of the UCT-MH to
guide the exact FCO. Multiple processors are allocated to the FCO algorithm
that systematically navigates the tree to solve the problem exactly.
Meanwhile, one processor is dedicated to run the UCT-MH which efficiently
samples the entire design space to get an estimate of the associated costs.
Considering the parallelization overhead of the exact FCO algorithm, it can be
expected that the UCT-MH already finds a fleet composition candidate
$\mathbf{\hat{k}}$ by the time the FCO is initiated. If such a composition is
available, then it is used as the candidate fleet
$\mathbf{f}^{1}=\mathbf{\hat{k}}$ that initializes the FCO algorithm.
Moreover, whenever the guiding UCT-MH finds a new best solution, it provides
this solution with its associated cost to the guided FCO by adding it to the
pooled best cost shown in Fig. 1. This information is used to preemptively
prune sub-optimal branches whenever possible and helps guide the FCO towards
the optimal fleet size and composition, reducing the search space and
computation time.

##  III Results

(a) The quality of candidate fleet sizes as determined by the UCT-MH

(b) The UCT-MH balances exploration and exploitation of candidate solutions.

(c) Performance of UCT-MH and B&B algorithm and their parallelization.

Figure 3: Case Study 1: The UCT-MH and B&B algorithm, number of tasks $n=10$,
maximum number of AMRs: $m_{max}=6$, $\mathbf{k}_{max}^{\top}=[2,2,2]$.

(a) The quality of candidate fleet sizes as determined by the UCT-MH, first 60
minutes of simulation.

(b) Performance of UCT-MH and B&B algorithm and their parallelization.

Figure 4: Case Study 1: The UCT-MH and B&B algorithm, number of tasks $n=20$,
maximum number of AMRs: $m_{max}=21$, $\mathbf{k}_{max}^{\top}=[7,7,7]$.

###  III-A Computational Experiments

To study the performance of the proposed hybrid algorithm, the guiding UCT-MH
and the guided FCO are compared against the standalone FCO. Four real-life
case studies are conducted in MATLAB 2022a at the Ohio Supercomputer Center
[21]. For each experiment, a set of $n$ tasks is defined, each consisting of a
unique item of known mass, pickup and drop off locations and respective time
windows. The fleet size is limited to $m_{max}$, equally distributed over
$h=3$ different AMR types. Each algorithm is allowed to run for a limited time
after which the incumbent solutions are compared. Firstly, two smaller
problems are studied in detail to illustrate the behavior of the UCT-MH in
Fig. 4-4. The best-known costs by each algorithm are summarized for all case
studies in Table I.

###  III-B Case Studies

####  III-B1 $n=10$ and $m_{max}=6$

Figure 3a shows the UCT-MH exploration of the various fleet sizes, where the
mean of the cost-to-go starts to converge and the algorithm gains more
confidence in particular solutions as the number of evaluations increases. The
guiding UCT-MH finds that $m=6$ is the best candidate and dedicates more
visits to these branches as shown in Fig. 3b. As a result, the guided FCO
quickly focuses on local convergence (Fig. 3c). As the entire search space is
explored, this solution is the guaranteed global optimum.

####  III-B2 $n=20$ and $m_{max}=21$

In Fig. 3a several patterns are observed. While small fleet sizes yield
infeasible solutions, larger fleet sizes initially show a transient behavior
due to the stochastic exploration. The largest fleet sizes always yield
feasible solutions, irrespective of the lower level decisions. Here, an
increase in fleet size results in an incremental increase of the mean cost-to-
go which is associated with the fleet cost. Remarkably, Fig. 3b shows that the
standalone FCO is initially faster, however, as the guided FCO already starts
from a good candidate branch, the underlying TSPTW is expected to be more
difficult to solve. Consequently, the guided FCO discards suboptimal fleets
and focuses on local convergence thereby reducing the overall computation time
of the guided FCO.

###  III-C Discussion

The time taken to initialize the parallel FCO algorithm is sufficient for the
guiding UCT-MH to find a strong candidate fleet that warm starts the guided
FCO. The UCT-MH provides a reduction of computation time ranging from $38.3\%$
up to $86.5\%$. The performance of the UCT-MH is dependent on the problem size
due to the time cap imposed at the TSPTW level. Thus, for higher number of
tasks where the TSPTW is larger, the gap with the incumbent solution is
greater. However, the guided FCO is able to close this gap since it conducts
local searches systematically. Further, for the case with 100 tasks, the
standalone FCO was unable to find any feasible solution in 24 hours while the
UCT-MH provided multiple solutions by its efficient stochastic exploration of
the design space.

Table I: Experimental Results - B&B and UCT-MH

Case Study | Standalone FCO | Guided FCO | Guiding UCT-MH |  |
---|---|---|---|---|---
$\mathbf{n}$ | $\mathbf{m}_{\text{max}}$ | t${}_{\text{max}}$ [min] | Cost | t${}_{\text{found}}$ [min] | Cost | Rel. Gap | t${}_{\text{found}}$ [min] | Reduction | Cost | Rel. Gap | t${}_{\text{found}}$ [min] |  |
10 | 6 | 2 | 9633.5∗ | 16 | 9633.5∗ | 0.00% | 2.17 | 86.5% | 10190.8 | 5.79% | 78.52 |  |
20 | 21 | 4 | 24,797.0 | 229.02 | 24,761 | \- 0.15% | 141.35 | 38.3% | 24,891.5 | 0.53% | 27.88 |  |
50 | 30 | 12 | 38,047.3 | 420.43 | 40,638.3 | \+ 6.81% | 212.30 | 49.5% | 53,161.3 | 39.72% | 2.883 |  |
100 | 60 | 24 | N/A† | 1440 | 74,250.5 | $-$ | 358.33 | $-$ | 103,193.0 | 38.98% | 30.05 |  |

∗Globally optimal solution.
†No solution was found after 24 hours.

##  IV Conclusions

In this paper, a hybrid optimization algorithm was developed that uses a
Monte-Carlo Tree Search-based metaheuristic (UCT-MH) to guide an exact fleet
composition optimization (FCO) algorithm. The FCO is a parallelized exact
Branch & Bound algorithm which solves a real-life Fleet Size and Mix Vehicle
Routing Problem with Time Windows. The UCT-MH yields a significant improvement
in the computation time and convergence of the FCO by constantly sharing the
expected optimal fleet composition as well as the upper bound on the cost.
Although in this study MCTS was only employed at the fleet sizing and
composition level, future research needs to determine to what depth MCTS can
be effective. Moreover, modifications to the selection policy as well as a bi-
directional communication between the UCT-MH and the FCO algorithm could
further improve computations times.

## Acknowledgments

This research was supported by the Ford Motor Company as part of the Ford-OSU
Alliance Program.

## References

  * [1] J. Morgan, M. Halton, Y. Qiao, and J. G. Breslin, “Industry 4.0 smart reconfigurable manufacturing machines,” pp. 481–506, 4 2021.
  * [2] M. G. Mehrabi, A. G. Ulsoy, and Y. Koren, “Reconfigurable manufacturing systems: Key to future manufacturing,” _Journal of Intelligent Manufacturing_ , vol. 11, pp. 403–419, 2000.
  * [3] Z. Ghelichi and S. Kilaru, “Analytical models for collaborative autonomous mobile robot solutions in fulfillment centers,” _Applied Mathematical Modelling_ , vol. 91, pp. 438–457, 3 2021.
  * [4] R. Yan, L. Jackson, and S. Dunnett, “A study for further exploring the advantages of using multi-load automated guided vehicles,” _Journal of Manufacturing Systems_ , vol. 57, pp. 19–30, 10 2020.
  * [5] A. Hoff, H. Andersson, M. Christiansen, G. Hasle, and A. Løkketangen, “Industrial aspects and literature survey: Fleet composition and routing,” _Computers and Operations Research_ , vol. 37, no. 12, pp. 2041–2061, 12 2010.
  * [6] G. Desaulniers, O. B. G. Madsen, and S. Ropke, “The Vehicle Routing Problem with Time Windows,” in _Vehicle Routing Problems, Methods, and Applications._ , 2nd ed., 2014, pp. 119–159.
  * [7] R. Elshaer and H. Awad, “A taxonomic review of metaheuristic algorithms for solving the vehicle routing problem and its variants,” _Computers and Industrial Engineering_ , vol. 140, 2 2020.
  * [8] M. N. A. Wahab, S. Nefti-Meziani, and A. Atyabi, “A comparative review on mobile robot path planning: Classical or meta-heuristic methods?” _Annual Reviews in Control_ , vol. 50, pp. 233–252, 1 2020.
  * [9] I. Boussaïd, J. Lepagnot, and P. Siarry, “A survey on optimization metaheuristics,” _Information Sciences_ , vol. 237, pp. 82–117, 7 2013\.
  * [10] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, “Mastering the game of Go with deep neural networks and tree search,” _Nature_ , vol. 529, no. 7587, pp. 484–489, 1 2016.
  * [11] M. Świechowski, K. Godlewski, B. Sawicki, and J. Mańdziuk, “Monte Carlo Tree Search: a review of recent modifications and applications,” _Artificial Intelligence Review_ , 2022.
  * [12] A. Sabharwal, H. Samulowitz, and C. Reddy, “Guiding Combinatorial Optimization with UCT,” in _International Conference on Integration of Artificial Intelligence (AI) and Operations Research (OR) Techniques in Constraint Programming_. Springer, 6 2012, pp. 356–361.
  * [13] B. Kartal, E. Nunes, J. Godoy, and M. Gini, “Monte Carlo Tree Search for Multi-Robot Task Allocation,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2016, pp. 4222–4223.
  * [14] S. Edelkamp, M. Gath, C. Greulich, M. Humann, O. Herzog, and M. Lawo, “Monte-Carlo Tree Search for Logistics,” in _Lecture Notes in Logistics_. Springer Cham, 2015, pp. 427–440.
  * [15] C. Barletta, W. Garn, C. Turner, and S. Fallah, “Hybrid fleet capacitated vehicle routing problem with flexible Monte–Carlo Tree search,” _International Journal of Systems Science: Operations and Logistics_ , 2022\.
  * [16] L. Kocsis and C. Szepesvári, “Bandit based Monte-Carlo Planning,” in _European Conference on Machine Learning_ , 9 2006, pp. 282–293.
  * [17] J. Lu, Y. Chen, J.-K. Hao, and R. He, “The time-dependent electric vehicle routing problem: Model and solution,” _Expert Systems with Applications_ , vol. 161, p. 113593, 2020.
  * [18] I. Kucukoglu, R. Dewil, and D. Cattrysse, “The electric vehicle routing problem and its variations: A literature review,” _Computers & Industrial Engineering_, vol. 161, p. 107650, 2021.
  * [19] M. Goutham, S. Boyle, M. Menon, S. Mohan, S. Garrow, and S. Stockar, “Optimal path planning through a sequence of waypoints,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 8566–8573, 2022.
  * [20] C. B. Browne, E. Powley, D. Whitehouse, S. M. Lucas, P. I. Cowling, P. Rohlfshagen, S. Tavener, D. Perez, S. Samothrakis, and S. Colton, “A survey of Monte Carlo tree search methods,” pp. 1–43, 3 2012.
  * [21] O. S. Center, “Ohio supercomputer center,” 1987. [Online]. Available: http://osc.edu/ark:/19495/f5s1ph73