# Heuristics for Multi-Vehicle Routing Problem Considering Human-Robot
Interactions

Venkata Sirimuvva Chirala, Kaarthik Sundar, Saravanan Venkatachalam, Jonathon
M. Smereka, Sam Kassoumeh Venkata Sirimuvva Chirala and Saravanan
Venkatachalam are with the Department of Industrial and Systems Engineering,
Wayne State University, Detroit, MI, 48202 USA (e-mail<EMAIL_ADDRESS>fz1533@wayne.edu).Kaarthik Sundar is with Los Alamos National Labaratory, Los
Alamos, New Mexico, 87545 USA (e-mail: kaarthik@lanl.gov).Jonathon M. Smereka
and Sam Kassoumeh are with US Army CCDC, Ground Vehicles Systems Center, MI,
48397 USA (e-mail<EMAIL_ADDRESS>sam.a.kassoumeh.civ@army.mil).

###### Abstract

Unmanned ground vehicles (UGVs) are being used extensively in civilian and
military applications for applications such as underground mining, nuclear
plant operations, planetary exploration, intelligence, surveillance and
reconnaissance (ISR) missions and manned-unmanned teaming. We consider a
multi-objective, multiple-vehicle routing problem in which teams of manned
ground vehicles (MGVs) and UGVs are deployed respectively in a leader-follower
framework to execute missions with differing requirements for MGVs and UGVs
while considering human-robot interactions (HRI). HRI studies highlight the
costs of managing a team of follower UGVs by a leader MGV. This paper aims to
compute feasible paths, replenishments, team compositions and number of MGV-
UGV teams deployed such that the requirements for MGVs and UGVs for the
missions are met and the path, replenishment, HRI and team deployment costs
are at minimum. The problem is first modeled as a a mixed-integer linear
program (MILP) that can be solved to optimality by off-the-shelf commercial
solvers for small-sized instances. For larger instances, a variable
neighborhood search algorithm is offered to compute near optimal solutions and
address the challenges that arise when solving the combinatorial multi-
objective routing optimization problem. Finally, computational experiments
that corroborate the effectiveness of the proposed algorithms are presented.

###### Index Terms:

unmanned autonomous vehicles, human-robot interaction, vehicle routing,
variable neighborhood search, multi-objective optimization

##  I Introduction

Industry 4.0 [1] concepts and technologies such as artificial intelligence
(AI), machine learning, cloud computing are establishing their presence across
innumerable fields in the recent years. Worldwide artificial intelligence (AI)
software revenue is forecast to total $62.5 billion in 2022, an increase of
21.3% from 2021 [2]. Advancement in AI technology is causing a paradigm shift
in the way autonomous vehicle technology is applied across various civilian
and military applications. UGVs and unmanned aerial vehicles (UAVs),
previously controlled manually or remotely, with the incorporation of AI are
able to operate with increasing autonomy. The use of UGVs and UAVs to assist
humans can be seen across various industries such as agricultural spraying and
harvesting, nuclear plant operations, crowd control, environmental monitoring,
and ISR missions such as explosive ordnance disposal (EOD), firefighting,
search and rescue.

Figure 1: A piece-wise linear discontinuous function for the team’s
performance based on a HRI study in [3]

In 2017, the United States Army published its strategy for the development and
deployment of robotic and autonomous systems [4]. This strategy emphasizes
primarily on human-machine collaboration by integrating army with robotic and
autonomous systems solutions. Teaming humans/MGVs with robots (UGVS and UAVs)
and AI will provide significant advantage to military forces by increasing
combat effectiveness, lightening the soldiers’ cognitive and physical
workloads, facilitating maneuvering to difficult terrains and reducing harm to
human personnel. MGV-UGV teams in a leader-follower framework are deployed in
semi-autonomous convoy operations thereby providing logistics formations the
ability to rest drivers for critical tasks only humans/MGVs can perform [4]. A
convoy is a group of vehicles travelling together accompanied by armed troops
or other vehicles. Within the leader-follower framework, the leader refers to
an MGV. The MGV leader lays down “virtual train tracks” which a set of
follower UGVs use to follow the leader. Determining ideal team compositions,
i.e., the ratio of human operators to robots is crucial for missions involving
MGV-UGV teams [5], [6]. The team compositions used in this effort are based on
HRI studies. HRI is a rapidly evolving field. The ideal team compositions for
MGV-UGV teams are derived from various characteristics such as neglect
tolerance and activity time. The autonomy of a robot is measured by its
neglect time. A human operator can ignore a robot for the duration of its
neglect time. Activity time measures the amount of time a human operator is
involved in instructing a robot [7]. These HRI studies provide a model for the
human costs of managing a team of UGVs by MGVs [8]. Authors in [3] conducted
an urban search and rescue mission (USAR) to evaluate the task performance
while a human operator was carrying varying number of robots. Fig. 1 explains
the HRI/team composition costs based on [3]. Authors in [7] also measure task
effectiveness as a function of the number of robots carried by a human
operator.

In this research, we consider a multiple vehicle routing problem (MVRP) in
which MGV-UGV teams are deployed to implement missions requiring human-robot
collaboration. Given a set of point of interests (POIs) with differing
requirements for MGVs and UGVs, the objective is to find ideal number of
teams, paths and replenishments for each team, and team compositions such that
the following conditions are satisfied:

  1. 1.

The demand for the MGVs and UGVs at the POIs are satisfied. The demand for the
UGVs is met either by the MGV-UGV team travelling to the POIs or by a
combination of the former and additional replenishments supplied to the POI
from the depot that then join the MGV-UGV team at the POI.

  2. 2.

The MGV-UGV teams follow a leader-follower framework where follower UGVs
follow a leader MGV based on HRI studies with an exception for the
replenishment UGVs supplied from the depot until they join the MGV-UGV team at
a POI.

  3. 3.

The sum of travel (path and replenishment) costs, HRI (team composition) costs
and costs for the number of teams deployed to complete the missions are
minimum.

(a)

(b)

Figure 2:  (a) shows the location of depot, POIs and the demand for UGVs at
the POIs. (b) shows feasible routes such that all POIs are visited, and
feasible assignments such that the demand for the UGVs at the POI is
satisfied: replenishment UGVs from depot (blue arrows) and total UGVs
traversing the edges (green arrows), given two MGV-UGV teams.

Given two vehicles, seven POIs and the demand for UGVs at POIs indicated in
Fig 2a, a feasible solution for the MVRP instance is shown in Fig 2b. In the
next section, we present a detailed overview of the literature that concerns
this problem.

##  II Related Work

The multiple vehicle, MGV-UGV team routing problem described in this paper is
a combination of the capacitated vehicle routing problem (CVRP) [9] and the
generalized assignment problem (GAP) [10]. The CVRP problem determines the
optimal routes for the MGV-UGV teams and the GAP problem is used to determine
the optimal team compositions for each team; in this case, the costs for the
team compositions are based on HRI studies.

The CVRP is a variant of the vehicle routing problem (VRP) [11] which
determines the optimal routes for a fleet of capacitated vehicles initially
located at one or more depots to serve a set of customers [12]. There has been
extensive research conducted in literature to apply exact methods [13, 14, 15]
that find the optimal solution to the CVRP using methods such as branch-and-
cut. Heuristics [16, 17] that solve the problem faster and efficiently,
however compromising on optimality, have also been explored. Given $n$ agents
and $m$ tasks, the GAP determines optimal allocation of agents to tasks, such
that each agent executes only one task subject to capacity restrictions on the
agents [10]. Similar to the CVRP, there is extensive research on exact methods
[18, 19, 20] and heuristics [21, 22] for GAP.

In our earlier work presented in [23], we propose an optimization model to
solve the MVRP problem. This was a first attempt at proposing a solution to
multi-vehicle route planning problems involving MGV-UGV teams taking human-
robot interaction into consideration. The optimization model is presented in
this paper in Sec. III to keep the presentation self-contained. The CVRP
problem and the GAP problem, the two key sub-problems of the MVRP, are NP-Hard
[24]. Therefore, scaling the optimization model to solve large scale instances
is very challenging, time consuming and computationally expensive. Hence, the
use of heuristic methods becomes inexorable.

VNS, proposed by Mladenovic and Hansen [25], is a very well-known
metaheuristic used to solve such NP-hard optimization problems. VNS utilizes
systematic change of neighborhoods within a local search algorithm exploring
increasingly distant neighborhoods of the current incumbent solution and
moving to a new solution only if an improvement is made [26]. Authors in [27]
use the VNS with seven neighborhoods, three intra-route and four inter-route
neighborhoods, alongside the guided local search metaheuristic to solve
capacitated vehicle routing problems with a single depot and $2,400$-$20,000$
customers. They show that the proposed methodology is able to produce fast,
high-quality feasible solutions for numerous large-scale benchmark instances
within reasonable time. Authors in [28] apply VNS algorithm to multi-depot
vehicle routing problem with time windows for instances containing up to $300$
customers. They compare the performance of VNS to Tabu-search heuristic [29]
and conclude that VNS is able to perform better with solution quality and
scaling. VNS has also been an efficient heuristic for assignment problems [30,
31]. Adaptation of VNS to multi-objective combinatorial optimization problems
is gaining pace in the research community. Authors in [32] suggest an adapted
VNS approach to solve two combinatorial bi-objective optimization problems.
They report that the suggested approach outperforms both Non-dominated Sorting
Genetic Algorithm (NGSA) and the Strength Pareto Evolutionary Algorithm, two
standard and popular algorithms used multi-objective optimization. Authors in
[33] apply the multi-objective VNS to a path planning problem in mobile
robotics by optimizing three objectives, namely path safety to avoid
obstacles, path length to control robot operation time and path smoothness to
optimize energy consumption to obtain accurate and efficient paths. They apply
VNS to eight different realistic scenarios and conclude that the multi
objective VNS approach is superior in comparison to the NSGA. Given the
efficacy of VNS in solving various classes of vehicle routing problems
including problem with multiple vehicles, multiple depots, time windows,
capacity restrictions and multi objective combinatorial optimization problems,
in this research, we propose the variable neighborhood search (VNS) algorithm
approach to the MVRP problem and address the challenge of scaling the approach
to large test instances.

The remainder of the paper is organized as follows. The optimization model
used to solve the MVRP problem is presented in section III. The variable
neighborhood search algorithm and approach are presented in section IV.
Computational results for randomly generated instances are presented in
section V. Finally, the paper concludes in section VI.

##  III Optimization Model

The optimization model presented in our previous work in [23] is included in
this section with additional explanation for completeness.

###  III-A Problem Description

Let $P=\\{p_{1},\dots,p_{l}\\}$ denotes the set of POIs. POIs are the
locations at which MGV-UGV teams are required to execute missions requiring
human-robot collaboration. Let $d_{0}$ denote the depot where $K$ different
categories of heterogeneous MGVs and UGVs represented as $M_{k}$, $m\in
M_{k}$, are initially stationed. $V=P\cup d_{0}$ denotes the set of vertices
and $E$ denotes the set of edges joining any pair of vertices. For each edge
$(i,j)\in E$, $c_{ij}$ represents the travel cost incurred while traversing
the edge $(i,j)$. Set $N$ (indexed by $n$) is used to represent the number of
UGVs in a team, and the parameter $h_{n}$ represents the HRI costs of managing
‘$n$’ UGVs. Parameter $T$ represents the cost incurred for each MGV-UGV team
deployed to execute the missions. For any set $S\subset V$,
$\delta^{+}(S)=\\{(i,j)\in E:i\in S,j\notin S\\}$ and
$\delta^{-}(S)=\\{(i,j)\in E:i\notin S,j\in S\\}$ denote the set of outgoing
and incoming edges from the set $S$. When $S=\\{i\\}$, we shall simply write
$\delta^{+}(i)$ and $\delta^{-}(i)$ instead of $\delta^{+}(\\{i\\})$ and
$\delta^{-}(\\{i\\})$, respectively. Also, for any $A\subseteq E$, let
$x(A)=\sum_{(i,j)\in A}x_{ij}$.

The decision variables for the MVRP are as follows: each edge $(i,j)\in E$ is
associated with a variable $x_{ij}^{m}$ that equals $1$ if the edge $(i,j)$ is
traversed by the ‘$m$th’ MGV, and $0$ otherwise. Variable $y_{j}^{m}$
represents the number of UGVs visiting the POI $j\in P$ along with the MGV
$m$. Variable $\bar{y}_{nj}^{m}$ represents variable $y_{j}^{m}$ as a series
of binary variables. Variable $z_{j}^{m}$ represents the number of UGV’s
provided as reinforcement to the POI $j\in P$ from the base depot $d_{0}$ for
the MGV $m$. Variable $r^{m}_{nj}$ represents variable $z_{j}^{m}$ as a series
of binary variables. Also, variable $w_{ij}^{m}$ represents the number of UGVs
that visit a POI $i$ then a POI $j$, where $i,j\in P$ along with the MGV $m$.

###  III-B Objective Function

The objective of the problem is to minimize the sum of (1) routing costs, (2)
replenishment costs, (3) HRI costs and finally (4) the cost incurred by the
teams to the mission. The expressions for each of these costs functions are as
follows:

Routing costs:
$\displaystyle\mathcal{R}^{1}\triangleq\sum\limits_{\begin{subarray}{c}(i,j)\in
E,n\in N\\\ m\in M_{k},k\in K\end{subarray}}c_{ij}x_{ij}^{m}$ (1a)
Replenishement costs:
$\displaystyle\mathcal{R}^{2}\triangleq\sum\limits_{\begin{subarray}{c}(i,j)\in
E,n\in N\\\ m\in M_{k},k\in K\end{subarray}}c_{d_{0}j}{r}^{m}_{nj}$ (1b) HRI
costs: $\displaystyle\mathcal{H}\triangleq\sum\limits_{\begin{subarray}{c}j\in
V,n\in N\\\ m\in M_{k},k\in K\end{subarray}}h_{n}\bar{y}_{nj}^{m}$ (1c) Team
costs: $\displaystyle\mathcal{T}\triangleq\sum\limits_{\begin{subarray}{c}j\in
V,\\\ m\in M_{k},k\in K\end{subarray}}T\cdot x_{0j}^{m}$ (1d)

The objective function is the sum of the of above costs weighted by user
provided parameters $\alpha$, $\beta$, and $\gamma$ as follows:

$\displaystyle\min~{}\alpha\cdot(\mathcal{R}^{1}+\mathcal{R}^{2})+\beta\cdot\mathcal{H}+\gamma\cdot\mathcal{T}$
(2)

The weighted sum method is a classical approach used to solve multi-objective
optimization problems. The optimal solution to the problem depends on the
weights assigned to the objectives. Weights are assigned to the objective in
proportion to the preference factor assigned to the objective [34].
Alternatively, a problem instance can be evaluated for different weights to
give decision makers multiple perspectives.

###  III-C Routing Constraints

Routing constraints impose the requirements that all POIs $i\in P$ should be
visited and that the demand for the number of MGVs and UGVs at the POIs is
satisfied. Constraint (3a) ensures that a feasible solution is connected and
are referred to as sub tour elimination constraints. The optimization problem
is initially solved by relaxing these constraints. Once a feasible solution is
obtained, the callback feature offered by solvers such as CPLEX, Gurobi etc.
are used to add the constraints to the model if the constraints in (3a) are
violated. A constraint is said to be in violation if the tour does not contain
the depot $d_{0}$. Constraint (3b) states that the number of MGV-UGV teams
leaving the depot must come back to the depot after completing the missions at
the visited POIs. Constraint (3c) states that and MGV-UGV team entering a POI
must leave the POI. For each POI $i$, the pair of constraints in (3d) and (3e)
state that the demand for MGVs and UGVs at POI $i$ are satisfied.

$\displaystyle\begin{split}&x^{m}(S)\leqslant|S|-1\quad\forall(i,j)\in S,\\\
&S\subset V\setminus\\{d_{0}\\}:S\neq\emptyset,\quad\forall m\in M_{k},k\in
K,\end{split}$ (3a) $\displaystyle
x^{m}(\delta^{+}(d_{0}))=x^{m}(\delta^{-}(d_{0}))\quad\forall m\in M_{k},k\in
K,$ (3b) $\displaystyle x^{m}(\delta^{+}(i))=x^{m}(\delta^{-}(i))\quad\forall
i\in P,m\in M_{k},k\in K,$ (3c) $\displaystyle\begin{split}&\sum_{m\in
M_{k}}x^{m}(\delta^{+}(i))\geq a_{ik}\text{ and }\\\ &\sum_{m\in
M_{k}}x^{m}(\delta^{-}(i))\geq a_{ik}\quad\forall i\in P,k\in K,\end{split}$
(3d) $\displaystyle\sum_{m\in M_{k}}y^{m}_{j}\geq b_{jk}\quad\forall j\in
P,k\in K,$ (3e)

###  III-D HRI/Teaming Constraints

The teaming constraints decide the team sizes for the MGV-UGV teams. For each
POI, constraint (4a) states that the number of UGV’s along side MGV $m$
visiting the POI should be equal to the the incoming UGV’s from the previously
visited POI and the reinforcement provided to the current POI from the base
depot. Constraints in (4b)-(4d) state that UGVs’ travel alongside MGV $m$ if
and only if the edge $x_{ij}^{m}$ is traversed. These constraints are used to
convert the non-linear relationship ($\sum_{m\in M_{k}}w^{m}_{ij}=\sum_{m\in
M_{k}}y^{m}_{i}$ only if $x_{ij}^{m}$ = 1) to a linear one. Constraint (4e)
state the dependence of UGVs on an MGV following a leader-follower framework,
and constraints (4f) define the capacity for UGVs. Constraints (4g) state that
team size does not change for each manned GV ‘$m$’, and the constraints (4h)
are the special order sets of constraints (type 1) stating that only one
binary variable can be chosen. Constraint (4i) is used to account for
replenishment costs at each POI. Finally, constraint (4j) imposes restrictions
on the decision variables.

$\displaystyle\sum_{m\in M_{k},i\in P,d_{0}}w^{m}_{ij}+z^{m}_{j}=\sum_{m\in
M_{k}}y^{m}_{j}\quad\forall j\in(P,d_{0})$ (4a)
$\displaystyle\begin{split}&\sum_{m\in M_{k}}w^{m}_{ij}\leq\sum_{m\in
M_{k}}B*x^{m}_{ij}\\\ &\quad\forall i\in(P,d_{0}),\forall j\in(P,d_{0}),k\in
K,\end{split}$ (4b) $\displaystyle\begin{split}&\sum_{m\in
M_{k}}w^{m}_{ij}\leq\sum_{m\in M_{k}}y^{m}_{i}\\\ &\quad\forall
j\in(P,d_{0}),\forall i\in(P,d_{0}),k\in K,\end{split}$ (4c)
$\displaystyle\begin{split}&\sum_{m\in M_{k}}w^{m}_{ij}\geq\sum_{m\in
M_{k}}(y^{m}_{i}-(1-x^{m}_{ij})*B)\\\ &\quad\forall i\in(P,d_{0}),\forall
j\in(P,d_{0}),k\in K,\end{split}$ (4d) $\displaystyle 0\leqslant
y^{m}_{i}\leqslant\ell_{k}x_{ji}^{m}\quad\forall(i,j)\in E,m\in M_{k},k\in K,$
(4e) $\displaystyle 0\leqslant\sum_{m\in
M_{k}}y^{m}_{i}\leqslant\ell_{k}\quad\forall k\in K,$ (4f) $\displaystyle
y^{m}_{j}=\sum_{n\in N}n_{n}\bar{y}_{nj}^{m}\quad\forall j\in P,m\in
M_{k},k\in K,$ (4g) $\displaystyle\sum_{n\in N}\bar{y}^{m}_{nj}\leq
1\quad\forall m\in M_{k},k\in K,$ (4h) $\displaystyle z^{m}_{j}=\sum_{n\in
N}n_{n}{r}_{nj}^{m}\quad\forall j\in P,m\in M_{k},k\in K,$ (4i)
$\displaystyle\begin{split}&x_{ij}^{m}\in\\{0,1\\},y^{m}_{i}\in\mathbb{Z}^{+},w^{m}_{ij}\in\mathbb{Z}^{+},z^{m}_{i}\in\mathbb{Z}^{+},\\\
&\bar{y}^{m}_{nj},\bar{r}^{m}_{nj}\in\\{0,1\\},0\leqslant s_{d}^{m}\leqslant
1,\\\ &\quad\forall(i,j)\in E,n\in N,m\in M_{k},k\in K.\end{split}$ (4j)

The optimization model proposed to solve the MVRP contains the CVRP and GAP,
both of which individually are combinatorial optimization problems. They are
NP-Hard [24]. Therefore, the optimization model also is combinatorial. The
model can find optimal solution to small instances, the results of which are
presented in [23]. Scaling the optimization model to larger instances is
computationally challenging. In the next section, we present the variable
neighborhood search algorithm to address this challenge and compute high
quality sub-optimal solutions for the MVRP.

##  IV Variable Neighborhood Search

The skewed variable neighborhood search algorithm (SVNS) is a variant of the
VNS meta-heuristic that allows for the exploration of valleys far from the
incumbent solution [25]. SVNS is particularly advantageous when there may be
near optimal solutions in separated and distant valleys. The exploration of
far away valleys is done by allowing for ‘recentering’ of the search. A
solution outside of the incumbent solution is accepted if it lies far away
from the incumbent solution. Recentering requires a metric to define the
distance between the incumbent and current solutions. In this research, we use
the SVNS, a variant of VNS.

###  IV-A Construction Heuristic

SVNS starts with an initial feasible solution. Hence, we develop an initial
feasible solution to the MVRP. Construction heuristic start with an empty
solution to a given problem and extend the current solution repeatedly until a
complete solution is found [35]. The initial feasible solution to the MVRP
consists of two components: feasible routes i.e., the sequence of visits to
the POIs by MGV-UGV teams such that all POIs are visited and feasible
assignments i.e., the team compositions for MGV-UGV teams for each route.
Feasible assignments ensure that the demand for the UGVs at the POI is met
either by the incoming team to the POI from the previously visited POI or by
the replenishment UGVs supplied from the base depot. We first construct the
initial routes to the multiple travelling salesman problem (mTSP) by using an
extension of the Lin-Kernighan-Helsgaun heuristic [36], LKH-3 heuristic solver
[37]. We then construct feasible assignments to the MVRP by comparing the
following two assignment strategies and implementing the strategy that
provides a better solution:

  1. 1.

asgn-rule1: For each MGV-UGV team visiting a subset of POIs, consider the
maximum UGV demand from the subset and provide this UGV demand from the depot
allowing for no replenishments from the depot.

  2. 2.

asgn-rule2: For each MGV-UGV team visiting a subset of POIs, provide required
UGV demand and replenish POIs as and when required based on the UGV demand at
POIs.

(a)  (b)

(c)  (d)

Figure 3:  (a) shows the location of depot, POIs and the demand for UGVs at
the POIs for MGV-UGV team missions, (b) shows a feasible routes with 2 MGV-UGV
teams, (c) shows asgn-rule1 implementation, and (d) shows asgn-rule2
implementation.

The choice of assignment strategy implemented in the solution is dependent on
the weights $\alpha$ and $\beta$ used in the objective function for the MVRP
problem. asgn-rule1 strategy would provide a better solution when weight
$\alpha$ is prioritized over $\beta$, since no replenishments are allowed
leading to a replenishment cost of 0, however, increasing the HRI cost due to
the increase in team sizes. Alternatively, asgn-rule2 increases the
replenishment costs thereby reducing the team sizes and the HRI costs proving
effective when weight $\beta$ is prioritized over $\alpha$. Figure 3 explains
the assignment strategies. Finally, we reverse the routes generated by the
LHK-3 solver and construct feasible assignments once again as mentioned above.
The best solution overall is used as the initial feasible solution for the
SVNS. The costs for the construction heuristic solution may be high since the
routes obtained from the LKH-3 heuristic may not be the best with respect to
the assignments of the MGV-UGV teams. A different sequence of routes may yield
a better suited solution for assignments. In the next section, we define the
neighborhoods used in SVNS.

###  IV-B Neighborhoods

(a)  (b)

(c)  (d)

Figure 4:  (a) shows the locations of the depot and POIs and path for a MGV-
UGV team before applying intra neighborhood moves, (b) shows path after
applying reverse neighborhood, (c) shows path after applying POI-swap-intra
neighborhood, and (d) shows path after applying 2-opt-intra neighborhood.

(a)  (b)

(c)  (d)

Figure 5:  (a) shows the locations of the depot and POIs and paths for two
MGV-UGV teams before applying inter neighborhood moves, (b) shows paths after
applying POI-remove-insert neighborhood, (c) shows paths after applying POI-
swap-inter neighborhood, and (d) shows paths after applying POI-
SequenceExchange-inter neighborhood.

We define seven neighborhoods for the MVRP, three intra route neighborhoods,
namely (i) reverse, (ii) POI-swap-intra, (iii) 2-opt-intra and three inter
route neighborhoods, namely (iv) POI-remove-insert, (v) POI-swap-inter and
(vi) POI-SequenceExchange-inter. Given routes for the MGV-UGV teams, the
reverse neighborhood reverses the sequence of visits to POIs for each route.
The POI-swap-intra neighborhood selects randomly, two POIs from a route and
interchanges them with each other. 2-opt was first introduced by Croes[38].
The 2-opt-intra neighborhood makes use of the 2-opt swap mechanism where,
given a route and a segment of the route, the segment is reversed and
reinserted back in its place. The intra neighborhoods are explained in Fig. 4.
Given two routes, the POI-remove-insert takes a POI randomly from one route
and inserts it into a random position in another route. The POI-swap-inter
neighborhood takes a POI at random from two routes and exchanges the POIs.
Finally, the POI-SequenceExchange-inter neighborhood takes one segment each,
from two routes and and exchanges the segments. The inter route neighborhoods
are explained in Fig. 5. For assignments, all intra route and inter route
neighborhoods are evaluated with the two assignment strategies mentioned in
section IV-A and the best solution is used.

Additionally, we define an assignment neighborhood named (vii) asgn-rule3 as
follows: Given a MGV-UGV team visiting a set of POIs, the asgn-rule3
neighborhood meets the UGV demand at the POIs by allowing replenishments only
at a few selected POIs. For POIs where no replenishments are allowed the UGV
demand is carried from the MGV-UGV team coming into the POI. This neighborhood
allows for a balance of the replenishments and team sizes (replenishment and
HRI costs) for the MVRP, in contrast to the two assignment strategy asgn-rule1
and asgn-rule2 in section IV-A.

###  IV-C Algorithm

The construction heuristic solution obtained in section IV-A is given as the
initial feasible solution to SVNS. The SVNS algorithm consists of three steps:
(i) shaking, (ii) local search, and (iii) recentering. The POI-swap-intra
neighborhood presented in Sec. IV-B is used in the shaking phase of the
algorithm. The shaking phase is used to escape from the incumbent solution.
The 2-opt-intra neighborhood, all three inter route neighborhoods and the
asgn-rule3 neighborhood are used in local search phase of the algorithm. The
recentering step of the algorithm performs the following: if the solution from
the shaking and local search phase $f^{{}^{\prime\prime}}$ is better than the
incumbent solution computed thus far $f$, the search is recentered i.e.,
$f^{{}^{\prime\prime}}$ is assigned to $f$ and the new solution becomes the
incumbent solution. This is done in line 15 of the SVNS algorithm.
Additionally in the SVNS, line 23 of the algorithm, the search is also
recentered when $f^{{}^{\prime\prime}}$ is worse than $f$, provided the
distance between $f^{{}^{\prime\prime}}$ and $f$ is large enough [39]. Given
two feasible solutions where a solution represents the total cost for the
MVRP, a solution $f^{{}^{\prime\prime}}$ is said to be better than $f$ if the
total cost $f^{{}^{\prime\prime}}$ is strictly less than $f$ since the MVRP is
a cost minimization problem. Solution $f^{{}^{\prime\prime}}$ is said to be
worse than $f$ if the total cost $f^{{}^{\prime\prime}}$ is strictly greater
than $f$. The distance metric for the MVRP is defined as the relative gap
$(\%)$ in the total cost of the two solutions. Given two feasible solutions,
$f^{{}^{\prime\prime}}$ and $f$, the distance between them is said to be large
enough if the relative gap is greater than or equal to 20% and less than or
equal to 50%. If the search is recentered, a parameter $k$ is initialized to
1, otherwise, the parameter $k$ is increased by 1 unit. Only when the total
objective value decreases, the number of iterations with no improvement
(UnImproved) is initialized to 0. In all the other cases, UnImproved is
increased by 1. The pseudo-code of the SVNS is presented in Algorithm 1.

Algorithm 1 Skewed Variable Neighborhood Search

Input: $N_{j}$ for Shaking phase: POI-swap-intra
$N_{j}$ for Local Search Phase: 2-opt-intra, POI-remove-insert, POI-swap-
inter, POI-SequenceExchange-inter, asgn-rule3

1:f $\leftarrow$ construction heuristic solution.

2:while UnImproved $<$ UnImproved_max do

3: $k$ $\leftarrow$ 0

4: while $k$ $<$ $k\\_max$ do

5: $f^{{}^{\prime}}$ $\leftarrow$ Shake($f$) $\triangleright$ Shake

6: $j$ $\leftarrow$ 0

7: while j $<$ 5 do $\triangleright$ Local Search

8: Find best neighbor in $N_{j}(f^{{}^{\prime}})$ $\rightarrow$
$f^{{}^{\prime\prime}}$

9: if $f^{{}^{\prime\prime}}$ is better than $f^{{}^{\prime}}$ then

10: $f^{{}^{\prime}}$ $=$ $f^{{}^{\prime\prime}}$

11: $j$ $\leftarrow$ 0

12: else

13: $j$ $\leftarrow$ $j+$ 1

14: $f^{{}^{\prime\prime}}$ $=$ $f^{{}^{\prime}}$

15: if $f^{{}^{\prime\prime}}$ is better than $f$ then $\triangleright$
Recenter

16: $f$ $=$ $f^{{}^{\prime\prime}}$

17: $k$ $\leftarrow$ 0

18: UnImproved $\leftarrow$ 0

19: if $f$ is better than IncumbentSol then

20: IncumbentSol $\leftarrow$ $x$

21: else

22: UnImproved$\leftarrow$UnImproved$+$1

23: if $f^{{}^{\prime\prime}}$ is worse than and far from $f$ then

24: $f$ = $f^{{}^{\prime\prime}}$

25: $k$ $\leftarrow$ 0

26: else

27: $k$ $\leftarrow$ $k+$ 1

Output: IncumbentSol

##  V Computational Results

The SVNS algorithm for the MVRP is implemented using Python. All the
computational experiments were run on a Intel(R) Core(TM) i5-6200U PC running
at 2.30GHz and 4 GB RAM.

###  V-A Performance of SVNS algorithm

To show the overall performance and effectiveness of the algorithm, we
generate three sets of instances: 80 small scale instances with five POIs and
two vehicles; 80 medium scale instances with 20 POIs, three and five vehicles;
and 80 large scale instances with 40 POIs, six and eight vehicles. For each of
these instances, the path traversal costs and the replenishment costs are
measured as the Euclidean distance costs for the edges. The edge for
replenishment cost is the edge connecting the depot to the POI where
replenishment is supplied. The HRI (team compositions) are based on HRI
studies. We choose HRI costs based on Fig 1. The costs for the number of teams
deployed linearly increases as the number of vehicles deployed for the
missions increase. The demand for the UGVs is uniformly distributed between
one and 12 at the POIs. The total capacity for UGVs that can be carried by an
MGV is 12. The weights for the objective function $\alpha$, $\beta$, and
$\gamma$ can take values 0, 0.1, 0.3, 0.6, 1 such that the that the sum of
weights adds up to one.

(a)

(b)

(c)

Figure 6: Improvement of MVRP cost from construction heuristic (initial
feasible solution) cost using SVNS for (a) small scale instances, (b) medium
scale instances, (c) large scale instances.

The small scale instance runs are made for both the optimization model
presented in section III and the SVNS heuristic to evaluate the performance of
SVNS. The average run time for small scale instances for the optimization
model and SVNS are 4.28 and 2.02 seconds, respectively. SVNS results in
optimal solutions for all 80 small scale instances. The initial feasible
solution given by the construction heuristic in section IV-A alone produces
optimal solutions for 38 instances of the 80 small scale instances. For the
remaining 42 instances, the SVNS improves the initial feasible solution
provided by the construction heuristic by a maximum of 36.35% and an average
improvement of 11.17%. For the remainder of the paper, given 2 MVRP objective
function costs, $f^{{}^{\prime}}$ and $f^{{}^{\prime\prime}}$, improvement is
defined as the decrease from $f^{{}^{\prime}}$ to $f^{{}^{\prime\prime}}$ as a
percentage of $f^{{}^{\prime}}$. Fig 6a summarizes the performance of the SVNS
algorithm for the remaining 42 small scale instances.

The optimization model due to being combinatorial in nature did not scale to
find optimal solutions for medium and large scale instances. However with
SVNS, medium scale instances take an average time of 11.6 seconds and 1.32
minutes for instances with three and five vehicles, respectively, and large
scale instances take an average time of 5.79 and 13.21 minutes for instances
with six and eight vehicles, respectively, to obtain sub-optimal solutions.
SVNS improves the quality of the initial feasible solution by a maximum of
59.84% and an average of 25.51% for 56 of the 80 medium scale instances and a
maximum of 57.34% and an average of 25.48% for 35 of the 80 large scale
instances. Fig 6b and 6c summarize the performance of the SVNS algorithm for
SVNS improved medium and large scale instances.

Figure 7: Improvement of MVRP cost from construction heuristic (initial
feasible solution) cost using SVNS by varying parameters UnImproved_max and
$k\\_max$ in Algorithm 1

###  V-B Parametric Analysis

Parameters, UnImproved_max and $k\\_max$ are used in line 2 and 4 respectively
in the SVNS algorithm presented in Algorithm 1. We generate 30 random
instances with 10 small, 10 medium and 10 large scale instances to test the
SVNS algorithm for different values of parameters, UnImproved_max and
$k\\_max$. Figure 7 summarizes the average improvement in the MVRP cost from
the construction heuristic using SVNS for differing values of the two
parameters. The average improvement improves as the values of the parameters
are increased. The lowest average improvement of 35.55% occurs when both
parameters UnImproved_max and $k\\_max$ are set to a value of 10. The average
improvements are comparable when UnImproved_max and $k\\_max$ are set to 40
and 30 with 40.45% (case 1) and 40 and 60 with 40.76% (case 2) respectively.
The time taken by SVNS increases by 4.01% for the 10 large scale instances
between case 1 and case 2. Based on the results, we set the parameter values
of UnImproved_max to 40 and $k\\_max$ to 30. Increasing the values may yield
better solutions, however, compromising on the time taken to solve the
instances.

Figure 8: Instances resulting in optimal solution (indicated in blue color) from 80 small scale instances (indicated in orange color) when SVNS algorithm is used with one neighborhood in the local search phase at a time. TABLE I: SVNS improvement of MVRP cost from construction heuristic cost with individual neighborhoods in local search phase for small, medium and large scale instances  Instance set | Neighborhood | Number of instances improved | Maximum improvement | Average improvement
---|---|---|---|---
|  | (Total instances: 80) | (%) | (%)
Small scale | 2-opt-intra | 13 | 28.56 | 10.79
| POI-remove-insert | 37 | 31.81 | 11.73
| POI-swap-inter | 37 | 36.35 | 11.27
| POI-SequenceExchange-inter | 37 | 36.35 | 11.27
| asgn-rule3 | 14 | 28.56 | 9.77
Medium scale | 2-opt-intra | 45 | 51.51 | 21.53
| POI-remove-insert | 48 | 59.84 | 26.22
| POI-swap-inter | 47 | 59.84 | 27.28
| POI-SequenceExchange-inter | 50 | 61.35 | 24.39
| asgn-rule3 | 46 | 49.22 | 19.53
Large scale | 2-opt-intra | 31 | 43.74 | 18.77
| POI-remove-insert | 30 | 54.15 | 25.14
| POI-swap-inter | 31 | 55.66 | 24.44
| POI-SequenceExchange-inter | 30 | 54.89 | 24.78
| asgn-rule3 | 32 | 41.65 | 18.21

###  V-C Neighborhood Study

Figure 8 summarizes the number of instances that result in optimal solution
for small scale instances when SVNS algorithm is used with one neighborhood in
the local search phase at a time. SVNS with POI-swap-inter and POI-
SequenceExchange-inter neighborhoods result in optimal solutions for 65 of the
80 small scale instances. The POI-remove-insert results in optimal solutions
for 64 instances. All 3 neighborhoods are inter route neighborhoods and
perform comparatively. They are also better in comparison to the 2-opt-intra,
intra route neighborhood.

The performance of the SVNS algorithm for small, medium and large scale
instances for each neighborhood operator used in the local search phase of the
algorithm is summarized in Table I. The POI-remove-insert neighborhood gives
the highest average improvement for small and large scale instances. For
medium scale instances, the highest average improvement is observed with the
POI-swap-inter neighborhood. The performance of all three inter route
neighborhoods: POI-remove-insert, and POI-SequenceExchange-Inter used in the
local search of the SVNS are nearly identical. All the inter route
neighborhoods perform better than the 2-opt-intra route neighborhood with
small, medium and large scale instances. Overall, results suggest that the
inter route neighborhoods in the local search phase of the SVNS show superior
performance when compared to intra route neighborhood. The performance of the
asgn-rule3 neighborhood improves as the instance size increases. However, the
average improvement seen in asgn-rule3 neighborhood remains low as the
neighborhood does not perform any route modifications.

###  V-D Sensitivity of weights on multi-objective function

Figure 9: MVRP travel vs HRI cost plots for differing weights $\alpha$ and $\beta$ for five small scale instances (ss1, ss2, ss3, ss4 and ss5) to observe the effect of weights on the overall MVRP solution. TABLE II: Solutions for problem instances with differing weights ($\alpha$, $\beta$, $\gamma$) Instance | User Parameters | Travel Costs | HRI Cost | Total Cost | Run Time (sec)
---|---|---|---|---|---
| alpha | beta | gamma | Path Cost | Replenishment Cost |  |  | SVNS | opt
ss1 | 1 | 1 | 1 | 521.89 | 0 | 250.00 | 871.89 | 1.43 | 6
| 0.6 | 0.3 | 0.1 | 521.89 | 0 | 250.00 | 871.89 | 1.43 | 6
| 0.6 | 0.1 | 0.3 | 454.15 | 0 | 410.00 | 964.15 | 1.84 | 5
| 0.3 | 0.6 | 0.1 | 521.89 | 0 | 250.00 | 871.89 | 1.29 | 6
| 0.1 | 0.6 | 0.3 | 533.59 | 196.34 | 200.03 | 1029.96 | 1.39 | 5
| 0.1 | 0.3 | 0.6 | 510.47 | 89.56 | 220.03 | 920.06 | 1.26 | 6
| 0.3 | 0.1 | 0.6 | 454.15 | 0 | 410.00 | 964.15 | 1.65 | 5
ss2 | 1 | 1 | 1 | 226.60 | 36.50 | 190.02 | 553.14 | 0.95 | 7
| 0.6 | 0.3 | 0.1 | 226.60 | 36.50 | 190.02 | 553.14 | 0.63 | 7
| 0.6 | 0.1 | 0.3 | 220.80 | 0 | 290.00 | 610.82 | 0.54 | 6
| 0.3 | 0.6 | 0.1 | 264.80 | 36.50 | 160.03 | 561.36 | 0.54 | 5
| 0.1 | 0.6 | 0.3 | 275.30 | 77.20 | 140.03 | 592.59 | 0.61 | 4
| 0.1 | 0.3 | 0.6 | 275.30 | 77.20 | 140.03 | 592.59 | 0.49 | 6
| 0.3 | 0.1 | 0.6 | 220.80 | 0 | 290.00 | 610.82 | 0.54 | 6
ss3 | 1 | 1 | 1 | 422.71 | 83.34 | 160.03 | 766.08 | 1.47 | 6
| 0.6 | 0.3 | 0.1 | 420.82 | 0 | 250.00 | 770.82 | 1.27 | 5
| 0.6 | 0.1 | 0.3 | 420.82 | 0 | 250.00 | 770.82 | 2.47 | 5
| 0.3 | 0.6 | 0.1 | 422.71 | 83.34 | 160.03 | 766.08 | 1.31 | 4
| 0.1 | 0.6 | 0.3 | 464.32 | 83.34 | 150.04 | 797.70 | 2.22 | 3
| 0.1 | 0.3 | 0.6 | 422.71 | 83.34 | 160.03 | 766.08 | 1.24 | 6
| 0.3 | 0.1 | 0.6 | 420.82 | 0 | 250.00 | 770.82 | 1.28 | 5

The weights for the objective function to the MVRP, given by $\alpha$,
$\beta$, and $\gamma$ can either be predetermined by the decision makers to
find the best solution based on the priority given to each of the objectives
or pareto solutions obtained by differing the weights can be offered to the
decision makers. In multi objective optimization (MOO), a solution is said to
be non dominated, pareto optimal if none of the objective functions can be
improved in value without degrading some of the other objective values [40].
The weighted sum method is a classical A priori method used to solve MOO
problems. The set of objectives are scalarized to give a single objective by
adding weights to each of the objectives. The weights $\alpha$, $\beta$, and
$\gamma$ are parameters for scalarization.

Fig 9 plots the travel (path + replenishment) costs vs the HRI costs for
differing weights $\alpha$, $\beta$, and $\gamma$ for a small-scale problem
instance (ss1, ss2, ss3, ss4 and ss5) with five POIs and two vehicles. Table
II details the complete solutions ss1, ss2 and ss3. A cost of 200 for each
instance is added to the Total Cost in Table II to account for the team cost
for deploying two vehicles. It can be observed that the travel costs (path and
replenishments costs) are lesser when weight $\alpha$ is prioritized over
$\beta$. This can be particularly advantageous when MGV-UGV teams are
maneuvering through unfamiliar and difficult terrains and providing
replenishments in such a scenario would be unsuitable and expensive. However,
this increases the team sizes for the MGV-UGV teams, thereby increasing
HRI/team composition costs. Alternatively, weight $\beta$ can be prioritized
to reduce load on the MGVs and reduce HRI costs. However, such scenarios will
have higher travel costs as observed. The best course of action for a given
problem instance can be implemented based on the interests of the decision
makers by studying the solution characteristics for the instance.

##  VI Conclusion

In this research, we propose the SVNS algorithm to the multi-objective routing
problem with multiple vehicles considering human robot interactions.
Specifically, MGV-UGV teams are deployed in a leader-follower (MGV leader -
UGV followers) framework to execute missions at POIs with differing
requirements for MGVs and UGVs. To the best of our knowledge, our previous
work in [23] was a first attempt to model this problem. The model developed is
a combinatorial optimization problem and is NP-Hard. SVNS is a variant of the
VNS meta-heuristic that enhances the exploration of far away valleys. The SVNS
addresses the challenges faced in scaling the previously modelled
combinatorial optimization problem for the MVRP in [23]. We generate three
sets of instances: small, medium and large scale to evaluate the performance
of SVNS. The SVNS produces optimal solutions for all small scale instances.
For medium and large scale instances, SVNS generates high quality sub-optimal
solutions for the optimization problem for more realistic and practical
scenarios in reasonable time. Our future work involves focusing on extending
the model to include uncertainties in the availability of UGVs, travel time,
etc. Additional functionalities such as vehicle refueling can be incorporated
into the model. Machine learning algorithms can be used to find the
appropriate values for the weights, $\alpha$, $\beta$, and $\gamma$ used in
the objective function.

## References

  * [1] C. Bai, P. Dallasega, G. Orzes, and J. Sarkis, “Industry 4.0 technologies assessment: A sustainability perspective,” _International journal of production economics_ , vol. 229, p. 107776, 2020.
  * [2] G. Inc. (2021) Gartner forecasts worldwide artificial intelligence software market to reach $62 billion in 2022. [Online]. Available: www.gartner.com/en/newsroom/press-releases/2021-11-22-gartner-forecasts-worldwide-artificial-intelligence-software-market-to-reach-62-billion-in-2022
  * [3] P. Velagapudi, P. Scerri, K. Sycara, H. Wang, M. Lewis, and J. Wang, “Scaling effects in multi-robot control,” in _2008 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2008, pp. 2121–2126.
  * [4] U.S. Army Training and Doctrine Command and Army Capabilities Integration Center, _The U.S. Army Robotic and Autonomous Systems Strategy_. U.S. Army, 2017.
  * [5] H. A. Yanco and J. L. Drury, “A taxonomy for human-robot interaction,” in _Proceedings of the AAAI fall symposium on human-robot interaction_ , 2002, pp. 111–119.
  * [6] Y. Liu and G. Nejat, “Robotic urban search and rescue: A survey from the control perspective,” _Journal of Intelligent & Robotic Systems_, vol. 72, no. 2, pp. 147–165, 2013.
  * [7] D. R. Olsen Jr and S. B. Wood, “Fan-out: Measuring human control of multiple robots,” in _Proceedings of the SIGCHI conference on Human factors in computing systems_ , 2004, pp. 231–238.
  * [8] M. A. Goodrich and A. C. Schultz, _Human-robot interaction: a survey_. Now Publishers Inc, 2008.
  * [9] T. K. Ralphs, L. Kopman, W. R. Pulleyblank, and L. E. Trotter, “On the capacitated vehicle routing problem,” _Mathematical programming_ , vol. 94, no. 2, pp. 343–359, 2003.
  * [10] D. G. Cattrysse and L. N. Van Wassenhove, “A survey of algorithms for the generalized assignment problem,” _European journal of operational research_ , vol. 60, no. 3, pp. 260–272, 1992.
  * [11] P. Toth and D. Vigo, _The vehicle routing problem_. SIAM, 2002.
  * [12] T. Paolo and V. Daniele, “Models, relaxations and exact approaches for the capacitated vehicle routing problem,” _Discrete Applied Mathematics_ , vol. 123, no. 1-3, pp. 487–512, 2002.
  * [13] G. Laporte and Y. Nobert, “Exact algorithms for the vehicle routing problem,” in _Surveys in Combinatorial Optimization_ , ser. North-Holland Mathematics Studies, S. Martello, G. Laporte, M. Minoux, and C. Ribeiro, Eds. North-Holland, 1987, vol. 132, pp. 147–184. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0304020808732353
  * [14] A. N. Letchford, J. Lysgaard, and R. W. Eglese, “A branch-and-cut algorithm for the capacitated open vehicle routing problem,” _Journal of the Operational Research Society_ , vol. 58, no. 12, pp. 1642–1651, 2007. [Online]. Available: https://doi.org/10.1057/palgrave.jors.2602345
  * [15] R. Fukasawa, H. Longo, J. Lysgaard, M. P. d. Aragão, M. Reis, E. Uchoa, and R. F. Werneck, “Robust branch-and-cut-and-price for the capacitated vehicle routing problem,” _Mathematical programming_ , vol. 106, no. 3, pp. 491–511, 2006.
  * [16] W. Y. Szeto, Y. Wu, and S. C. Ho, “An artificial bee colony algorithm for the capacitated vehicle routing problem,” _European Journal of Operational Research_ , vol. 215, no. 1, pp. 126–135, 2011.
  * [17] S.-W. Lin, Z.-J. Lee, K.-C. Ying, and C.-Y. Lee, “Applying hybrid meta-heuristics for capacitated vehicle routing problem,” _Expert Systems with Applications_ , vol. 36, no. 2, pp. 1505–1512, 2009.
  * [18] R. M. Nauss, “Solving the generalized assignment problem: An optimizing and heuristic approach,” _INFORMS Journal on Computing_ , vol. 15, no. 3, pp. 249–266, 2003.
  * [19] J. S. Park, B. H. Lim, and Y. Lee, “A lagrangian dual-based branch-and-bound algorithm for the generalized multi-assignment problem,” _Management Science_ , vol. 44, no. 12-part-2, pp. S271–S282, 1998.
  * [20] M. Savelsbergh, “A branch-and-price algorithm for the generalized assignment problem,” _Operations research_ , vol. 45, no. 6, pp. 831–841, 1997.
  * [21] M. Yagiura, T. Ibaraki, and F. Glover, “An ejection chain approach for the generalized assignment problem,” _INFORMS Journal on computing_ , vol. 16, no. 2, pp. 133–151, 2004.
  * [22] J. A. Dıaz and E. Fernández, “A tabu search heuristic for the generalized assignment problem,” _European Journal of Operational Research_ , vol. 132, no. 1, pp. 22–38, 2001.
  * [23] V. S. Chirala, S. Venkatachalam, J. M. Smereka, and S. Kassoumeh, “A multi-objective optimization approach for multi-vehicle path planning problems considering human–robot interactions,” _Journal of Autonomous Vehicles and Systems_ , vol. 1, no. 4, p. 041002, 2021.
  * [24] J. K. Lenstra and A. R. Kan, “Complexity of vehicle routing and scheduling problems,” _Networks_ , vol. 11, no. 2, pp. 221–227, 1981.
  * [25] P. Hansen and N. Mladenovic, “A tutorial on variable neighborhood search,” _Les Cahiers du GERAD ISSN_ , vol. 711, p. 2440, 2003.
  * [26] N. Mladenović and P. Hansen, “Variable neighborhood search,” _Computers & operations research_, vol. 24, no. 11, pp. 1097–1100, 1997\.
  * [27] J. Kytöjoki, T. Nuortio, O. Bräysy, and M. Gendreau, “An efficient variable neighborhood search heuristic for very large scale vehicle routing problems,” _Computers & Operations Research_, vol. 34, no. 9, pp. 2743–2757, 2007. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0305054805003394
  * [28] M. Polacek, R. F. Hartl, K. Doerner, and M. Reimann, “A variable neighborhood search for the multi depot vehicle routing problem with time windows,” _Journal of heuristics_ , vol. 10, no. 6, pp. 613–627, 2004.
  * [29] M. Gendreau, A. Hertz, and G. Laporte, “A tabu search heuristic for the vehicle routing problem,” _Management science_ , vol. 40, no. 10, pp. 1276–1290, 1994.
  * [30] S. Mitrović-Minić and A. P. Punnen, “Local search intensified: Very large-scale variable neighborhood search for the multi-resource generalized assignment problem,” _Discrete Optimization_ , vol. 6, no. 4, pp. 370–377, 2009.
  * [31] H. Gavranović, M. Buljubašić, and E. Demirović, “Variable neighborhood search for google machine reassignment problem,” _Electronic Notes in Discrete Mathematics_ , vol. 39, pp. 209–216, 2012.
  * [32] A. Duarte, J. J. Pantrigo, E. G. Pardo, and N. Mladenovic, “Multi-objective variable neighborhood search: an application to combinatorial optimization problems,” _Journal of Global Optimization_ , vol. 63, no. 3, pp. 515–536, 2015.
  * [33] A. Hidalgo-Paniagua, M. A. Vega-Rodríguez, and J. Ferruz, “Applying the movns (multi-objective variable neighborhood search) algorithm to solve the path planning problem in mobile robotics,” _Expert Systems with Applications_ , vol. 58, pp. 20–35, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0957417416301336
  * [34] K. Deb, “Multi-objective optimization,” in _Search methodologies_. Springer, 2014, pp. 403–449.
  * [35] Wikipedia. (2021) Constructive heuristic. [Online]. Available: https://en.wikipedia.org/wiki/Constructive˙heuristic
  * [36] K. Helsgaun, “An effective implementation of the lin–kernighan traveling salesman heuristic,” _European journal of operational research_ , vol. 126, no. 1, pp. 106–130, 2000.
  * [37] H. Keld. (2021) Lkh-3 version 3.0.7. [Online]. Available: http://webhotel4.ruc.dk/~keld/research/LKH-3/
  * [38] G. A. Croes, “A method for solving traveling-salesman problems,” _Operations research_ , vol. 6, no. 6, pp. 791–812, 1958.
  * [39] K. Sundar, S. G. Manyam, P. Sujit, and D. Casbeer, “Coordinated air-ground vehicle routing with timing constraints,” in _2019 Sixth Indian Control Conference (ICC)_ , 2019, pp. 116–121.
  * [40] Wikipedia. (2022) Multi-objective optimization. [Online]. Available: https://en.wikipedia.org/wiki/Multi-objective˙optimization


No relation for Wigner’s friend

Leonardo Castellani

Dipartimento di Scienze e Innovazione Tecnologica
Università del Piemonte Orientale, viale T. Michel 11, 15121 Alessandria,
Italy
INFN, Sezione di Torino, via P. Giuria 1, 10125 Torino, Italy
Arnold-Regge Center, via P. Giuria 1, 10125 Torino, Italy

We argue that Wigner’s friend thought experiment does not support observer
dependence of quantum states. An analysis in terms of history vectors suggests
that quantum collapse is to be understood as collapse of histories rather than
collapse of states.

<EMAIL_ADDRESS>
The core of subjective or relational interpretations of QM (see for ex. [1, 2,
3, 4]) is the alleged observer dependence of quantum states, meaning that
different observers give different descriptions of the same quantum system.
This is claimed by some authors to be a consequence of the usual quantum
mechanical rules, when these are applied to “third person” situations, as in
Wigner’s friend thought experiment [5]. In the present note we argue that the
usual textbook QM description of the experiment does not lead to an essential
observer dependence.

Let us briefly recall the setting: an observer $F$ and a system $S$ are
contained in a perfectly isolated laboratory. A second observer $W$ stays
outside of the laboratory. The initial state of $S$ at $t=t_{0}$ is prepared
in the superposition

$\alpha|0\rangle+\beta|1\rangle$ (1)

At $t=t_{1}$ the observer $F$ effects a measurement on $S$ in the
$|0\rangle,|1\rangle$ basis 111i.e. a measurement of an observable having
$|0\rangle$ and $|1\rangle$ as eigenvectors.. According to the usual QM rules,
$F$ obtains the result 0 with probability $|\alpha|^{2}$ and the result 1 with
probability $|\beta|^{2}$. Suppose for example that $F$ obtains 1: “for him”
the initial state $|\psi\rangle$ collapses into the state $|1\rangle$:

$\alpha|0\rangle+\beta|1\rangle\longrightarrow|1\rangle$ (2)

How does the external observer $W$ describe the evolution, from $t_{0}$ to
$t_{1}$, of the composite system $S+F$ ? According to QM, in absence of any
interaction with the environment, the system $S+F$ evolves unitarily. The
initial state of the composite system is the tensor product
$(\alpha|0\rangle+\beta|1\rangle)\otimes|init\rangle$, where $|init\rangle$ is
the initial state of the observer $F$. A measurement by $F$ on $S$ is modeled
by the $S+F$ (unitary) evolution

$\displaystyle|0\rangle\otimes|init\rangle\longrightarrow|0\rangle\otimes|0_{F}\rangle$
(3)
$\displaystyle|1\rangle\otimes|init\rangle\longrightarrow|1\rangle\otimes|1_{F}\rangle$
(4)

where $|0_{F}\rangle$ (resp. $|1_{F}\rangle$) is the state of $F$ after he has
obtained 0 (resp. 1), or equivalently the state of a notebook where $F$ writes
down the result of the measurement.

The crucial point is that the state of $S$ becomes correlated to the state of
$F$. Thus if the initial state of $S+F$ is
$(\alpha|0\rangle+\beta|1\rangle)\otimes|init\rangle$, by linearity the system
evolves from $t_{0}$ to $t_{1}$ as

$(\alpha|0\rangle+\beta|1\rangle)\otimes|init\rangle\longrightarrow\alpha|0\rangle\otimes|0_{F}\rangle+\beta|1\rangle\otimes|1_{F}\rangle$
(5)

i.e. in an entangled state of the $S+F$ system. All predictions on future
measurements by $W$ must be computed using this entangled state. In particular
$W$ could observe interference effects in $S+F$, since the state at $t_{1}$ is
a superposition.

The observer $F$, on the contrary, ”sees” no superposition. In fact he cannot
observe $S+F$, since he is part of the composite system. Only $S$ is
accessible to him, and after his measurement at $t=t_{1}$, $F$ describes the
state of $S$ as one of the two states $|0\rangle$, $|1\rangle$.

So far no difficulty arises. $W$ and $F$ are describing different physical
systems: no wonder that their descriptions differ.

But suppose that $W$ focuses his attention on the subsystem $S$. He can make
measurements on $S$, since $S$ is a subsystem of $S+F$ and therefore
accessible to him. How does $W$ describe the state of $S$ after $t_{1}$,
before making any measurement himself? Is this description different from the
one made by $F$ ? Here the question acquires relevance, since both observers
refer to the same system $S$. In fact it lies at the heart of the measurement
problem in QM, highlighting the clash between unitary evolution of isolated
systems and nonunitary projection due to measurement.

The answer of textbook QM involves the reduced density operator $\rho^{(S)}$
for system $S$. The density operator $\rho^{(S+F)}$ for $S+F$ at $t=t_{1}$ is
given by:

$\rho^{(S+F)}=(\alpha|0\rangle\otimes|0_{F}\rangle+\beta|1\rangle\otimes|1_{F}\rangle)(\alpha^{\star}\langle
0|\otimes\langle 0_{F}|+\beta^{\star}\langle 1|\otimes\langle 1_{F}|)$ (6)

The reduced density operator $\rho^{(S)}$ is defined by tracing on the
subsystem $F$. Assuming $|0_{F}\rangle$ and $|1_{F}\rangle$ to be
orthogonal222since they are eigenvectors of a “position of the pointer”
observable. we find

$\rho^{(S)}=Tr_{F}(\rho^{(S+F)})=|\alpha|^{2}|0\rangle\langle
0|+|\beta|^{2}|1\rangle\langle 1|$ (7)

Thus $W$ describes $S$ to be in a mixed state, and precisely in the state
$|0\rangle$ with probability $|\alpha|^{2}$ and in the state $|1\rangle$ with
probability $|\beta|^{2}$. Note that this mixed state is radically different
from the superposed (pure) state $\alpha|0\rangle+\beta|1\rangle$. It is the
description of system $S$ given by $W$, after a measurement has been performed
by $F$, but without acquisition of the result (since $S+F$ is isolated, $F$
cannot communicate with $W$).

In summary, the $F$ and $W$ descriptions of the state of $S$ at $t_{0}$ and
$t_{1}$ can be represented as follows:

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}t_{0}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}t_{1}$
$\displaystyle
F:~{}~{}~{}\alpha|0\rangle+\beta|1\rangle\longrightarrow~{}|0\rangle~{}or~{}|1\rangle~{}~{}with~{}probabilities~{}|\alpha|^{2}~{}or~{}|\beta|^{2}$
$\displaystyle
W:~{}~{}~{}\alpha|0\rangle+\beta|1\rangle\longrightarrow~{}mixed~{}state~{}(statistical~{}ensemble~{}of~{}|0\rangle~{}and~{}|1\rangle)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}with~{}probabilities~{}|\alpha|^{2}~{}and~{}|\beta|^{2}$

The two descriptions do not coincide, because $F$ knows the result of his
measurement, whereas $W$ does not. This is the only reason of the difference,
due to ignorance of $W$, and reflected in the statistical ensemble he must use
to describe $S$ at $t_{1}$. But both $F$ and $W$ do agree that the initial
state of $S$ has collapsed in one of the two basis states $|0\rangle$ or
$|1\rangle$, even if only $F$ knows which one. This difference has no profound
significance: $F$ would describe $S$ in exactly the same way as $W$, if he
chose to measure $S$ (producing collapse) without registering the result.
Moreover the difference is not due to quantum mechanical effects, but only to
incomplete information of one of the two observers. The same difference arises
in a classical world, where $S$ is for example a coin under a cup, and $F$
“measures” it by removing the cup and registering whether $S$ shows head or
tails. If $W$ cannot communicate with $F$ ($S+F$ being isolated), the
information on the result of $F$ measurement is not available, and he can only
describe the state of $S$ statistically, with a 50-50 probability of heads or
tails. But this could hardly be a motivation to develop a relational
interpretation of classical mechanics333besides the obvious system of
reference dependence of position, velocities etc..

To acquire the missing information, $W$ must interact with the $S+F$ system,
for example with a measurement on $S$, or on $F$, or on the whole $S+F$. If he
interacts with the whole of $S+F$, he will be able to detect that the
composite system is in the superposition state (5). But if he limits his
measurements to the subsystems $S$ or $F$, these are not in superposed (pure)
states, but in mixed states, and no interference effects can be detected.

Suppose that $W$ measures $S$ and obtains the result 0. The $S+F$ system
collapses:

$\alpha|0\rangle\otimes|0_{F}\rangle+\beta|1\rangle\otimes|1_{F}\rangle\longrightarrow|0\rangle\otimes|0_{F}\rangle$
(8)

and by a subsequent measure on $F$, $W$ can verify that 0 was indeed the
result registered by $F$. Thus no contradiction can arise between the results
of $F$ and $W$.

But suppose that $F$ has obtained 1 in his measurement on $S$. How is it
possible that $W$ has still a probability $|\alpha|^{2}$ of obtaining 0 on the
same system, after the measurement of $F$ ? This “paradox” is resolved
operationally, since Wigner’s measurement will never reveal any disagreement
with his friend’s measurement.

In fact this situation illustrates the “history collapse” effect of
measurements, discussed in [6, 7]. Consider the circuit version of Wigner’s
friend experiment in its simplest version, where $S$ is a single qubit system
and $F$ is a single qubit “observer”. The evolutions (3), (4) can be realized
just with a CNOT gate444The CNOT is a two-qubit gate acting on the
computational basis as as $|0\rangle|0\rangle\rightarrow|0\rangle|0\rangle$,
$|0\rangle|1\rangle\rightarrow|0\rangle|1\rangle$,
$|1\rangle|0\rangle\rightarrow|1\rangle|1\rangle$,
$|1\rangle|1\rangle\rightarrow|1\rangle|0\rangle$, the first qubit being the
control and the second being the target. where $S$ is the control and $F$ is
the target qubit, and choosing $|init\rangle=|0\rangle$:

Fig. 1 Circuit describing the measurement of $F$ on $S$, as seen by $W$, and
the corresponding history diagram.

As discussed in [7], a quantum system can be described by a history vector
$|\Psi\rangle$ that encodes its whole time evolution:

$|\Psi\rangle=\sum_{\gamma}A(\gamma)|\gamma\rangle$ (9)

where $\gamma$ are histories with nonvanishing amplitudes. A short summary
follows, where we recall the definition of histories and corresponding
amplitudes.

By history we mean a sequence of possible measurement results, at discrete
times $t_{1},t_{2},...$ (time discretization is assumed for simplicity). The
measured observables can be different at each $t_{i}$. With usual Born rules
we can compute the probability $p(\psi,\gamma)$ of obtaining a particular
sequence $\gamma=\gamma_{1},\gamma_{2},...$ of results at times
$t_{1},t_{2},...$, if the system starts in state $|\psi\rangle$. We find

$p(\psi,\gamma)=Tr(C_{\psi,\gamma}C_{\psi,\gamma}^{\dagger})$ (10)

where the chain operator $C_{\psi,\gamma}$ is defined as

$C_{\psi,\gamma}=P_{\gamma_{n}}U(t_{n},t_{n-1})~{}P_{\gamma_{n-1}}~{}U(t_{n-1},t_{n-2})\cdots
P_{\gamma_{1}}~{}U(t_{1},t_{0})P_{\psi}$ (11)

with $P_{\psi}=|\psi\rangle\langle\psi|$; $P_{\gamma_{i}}$ are projectors on
eigensubspaces of observables corresponding to the eigenvalues $\gamma_{i}$
and $U(t_{i+1},t_{i})$ is the unitary evolution operator between times $t_{i}$
and $t_{i+1}$. If $\gamma_{n}$ is nondegenerate, the chain operator can be
written as:

$C_{\psi,\gamma}=|\gamma_{n}\rangle A(\psi,\gamma)\langle\psi|$ (12)

the amplitude $A(\psi,\gamma)$ being defined by

$A(\psi,\gamma)=\langle\gamma_{n}|P_{\gamma_{n-1}}~{}U(t_{n-1},t_{n-2})\cdots
P_{\gamma_{1}}~{}U(t_{1},t_{0})|\psi\rangle$ (13)

and the probability (10) is just the square modulus of the amplitude (13). The
generalization to a degenerate eigenvalue $\gamma_{n}$ is straightforward (see
[7]). The vectors $|\gamma\rangle$ in (9) are in 1-1 correspondence with all
possible histories $\gamma_{1},...\gamma_{n}$, and assumed to form an
orthonormal basis. Probabilities for (any sequences of) measurements can be
easily computed using appropriate projections of the history vector (9) [7].

We apply now this formalism to the circuit that models Wigner’s friend
experiment, given in Fig.1. Note that it represents the system $S+F$
accessible only to $W$: the measurement of $F$ on the subsystem $S$ is
represented by the unitary gate $CNOT$, and not by a projection. Starting from
an initial state $|\psi\rangle|0\rangle$ the history vector contains only two
histories, since the only nonvanishing amplitudes are

$\displaystyle A(00,00)=\langle 00|CNOT|00\rangle\langle
00|\psi,0\rangle=\alpha$ (14) $\displaystyle A(10,11)=\langle
11|CNOT|10\rangle\langle 10|\psi,0\rangle=\beta$ (15)

Suppose now that $W$ at time $t_{2}$ measures $F$ and obtains 1. This result
is compatible only with the history $(10,11)$. Its probability is given the
square modulus of the corresponding amplitude, i.e. $|\beta|^{2}$. The history
vector collapses to

$|\Psi\rangle=|10,11\rangle$ (16)

Even if the measurement is performed at time $t_{2}$, the history involves
also values at time $t_{1}<t_{2}$: in this sense the collapse “modifies the
past”. Even if at $t_{1}$ $F$ had obtained 0 in his measurement of $S$, the
measurement of $W$ “undoes” this event and realigns it to agree with his
outcome at $t=t_{2}$. This “undoing” of past events has been discussed by
various authors (see for ex. [8]), and is really a matter of interpretation,
since there is no possibility of “checking” the undoing. However note that
history collapse is not a matter of interpretation, if we require simultaneity
of collapse in entangled spacelike separated systems to be valid in every
reference frame [6].

In conclusion, at first sight it may appear that Wigner ($W$) and his friend
($F$) give different descriptions of $F$’s measurement of $S$. In particular
$F$’s description involves a collapse, whereas $W$ describes a unitary
evolution. $W$ could detect interference effects in measurements on $S+F$,
whereas $F$ sees no such effects after his measurement on $S$. But this
difference is only due to $W$ and $F$ considering different physical systems:
$W$ describes $S+F$, whereas $F$ can only describe $S$. If we compare instead
their descriptions of the same system $S$, the difference essentially vanishes
as discussed above. The only residual difference is due to lack of
information, not to quantum mechanical effects, and in our opinion does not
motivate any “relational” or subjective interpretation of quantum states.

Note:  in recent times Gedanken experiments [9, 10, 11] extending Wigner’s
friend setup have been proposed to probe the issue of observer independence. A
real experiment has been carried out [12] in 2019, showing violation of Bell-
type inequalities derived by assuming free-choice, locality and observer
independence. However extreme care is necessary to uncover all assumptions
made in real or Gedanken experiments. For example the analysis in ref. [9] has
been criticized by a number of authors [13, 14, 15], and inconsistences have
been pointed out. It is therefore not clear to us that these experiments could
directly support observer independence of quantum states or events.

## Acknowledgements

We thank Carlo Rovelli for correspondence on recent experiments. This work is
supported by the research funds of the Eastern Piedmont University and INFN -
Torino Section.

## References

  * [1] C. A. Fuchs, R. Schack, “Quantum-Bayesian coherence,” Rev. Mod. Phys. 85, 1693 (2013).
  * [2] N. D. Mermin, “Physics: QBism puts the scientist back into science,” Nature 507, 421–423 (2014).
  * [3] H. Everett, “Relative state formulation of quantum mechanics,” Rev. Mod. Phys. 29, 454–462 (1957).
  * [4] C. Rovelli, “Relational Quantum Mechanics,” International Journal of Theoretical Physics 35 (1996) 1637; arXiv:quant-ph/9609002.
  * [5] E. P. Wigner, in The Scientist Speculates (ed. Good, I. J.) 284–302 (Heinemann, 1961).
  * [6] L. Castellani, “History operators in quantum mechanics,” Int. J. Quant. Inf. 17, no. 08, 1941001 (2019) doi:10.1142/S0219749919410016 [arXiv:1810.03624 [quant-ph]].
  * [7] L. Castellani, “History entanglement entropy,” arXiv:2009.02331 [quant-ph].
  * [8] C. Brukner, “On the quantum measurement problem”. In R. Bertlmann and A. Zeilinger, editors, Quantum [Un]Speakables II: Half a Century of Bell’s Theorem, pages 95–117. Springer, 2017; arXiv:1507.05255.
  * [9] D. Frauchiger and R. Renner, “Quantum theory cannot consistently describe the use of itself,” Nature Communications 9, 3711 (2018).
  * [10] C. Brukner, “A no-go theorem for observer-independent facts, ” Entropy 20, 350, (2018), arXiv:1804.00749
  * [11] K. Bong, A. Utreras-Alarcón, F. Ghafari et al., “A strong no-go theorem on the Wigner’s friend paradox,” Nat. Phys. (2020). https://doi.org/10.1038/s41567-020-0990-x
  * [12] M. Proietti et al., “Experimental test of local observer independence,” Science Advances 20 Sep 2019: Vol. 5, no. 9, DOI:10.1126/sciadv.aaw9832, arXiv:1902.05080 [quant-ph].
  * [13] M. Araújo, “The flaw in Frauchiger and Renner’s Argument,” (2018) http://mateusaraujo.info/2018/10/24/the-flaw-in-frauchiger-and-renners-argument/
  * [14] A. Sudbery, “The hidden assumptions of Frauchiger and Renner,” Int. J. Quant. Foundations 5, 98 (2019), arXiv:1905.13248..
  * [15] R. E. Kastner, “Unitary-Only Quantum Theory Cannot Consistently Describe the Use of Itself: On the Frauchiger–Renner Paradox,” Found. Phys. 50, no. 5, 441 (2020) doi:10.1007/s10701-020-00336-6 [arXiv:2002.01456 [quant-ph]].