© 2016 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

–

DOI:10.1109/TASE.2015.2432746

# An Iterative Approach for Collision Free Routing and Scheduling in
Multirobot Stations

Domenico Spensieri, Johan S. Carlson, Fredrik Ekstedt, and Robert Bohlin All
authors are with the Geometry and Motion Planning Group, at the Fraunhofer-
Chalmers Research Centre for Industrial Mathematics, Gothenburg, 41288,
Sweden, e-mail<EMAIL_ADDRESS>
###### Abstract

This work is inspired by the problem of planning sequences of operations, as
welding, in car manufacturing stations where multiple industrial robots
cooperate. The goal is to minimize the station cycle time, _i.e._ the time it
takes for the last robot to finish its cycle. This is done by dispatching the
tasks among the robots, and by routing and scheduling the robots in a
collision-free way, such that they perform all predefined tasks. We propose an
iterative and decoupled approach in order to cope with the high complexity of
the problem. First, collisions among robots are neglected, leading to a min-
max Multiple Generalized Traveling Salesman Problem (MGTSP). Then, when the
sets of robot loads have been obtained and fixed, we sequence and schedule
their tasks, with the aim to avoid conflicts. The first problem (min-max
MGTSP) is solved by an exact branch and bound method, where different lower
bounds are presented by combining the solutions of a min-max set partitioning
problem and of a Generalized Traveling Salesman Problem (GTSP). The second
problem is approached by assuming that robots move synchronously: a novel
transformation of this synchronous problem into a GTSP is presented.
Eventually, in order to provide complete robot solutions, we include path
planning functionalities, allowing the robots to avoid collisions with the
static environment and among themselves. These steps are iterated until a
satisfying solution is obtained. Experimental results are shown for both
problems and for their combination. We even show the results of the iterative
method, applied to an industrial test case adapted from a stud welding station
in a car manufacturing line.

_Note to Practitioners_ —This article is motivated by the problem of planning
robot operations in welding applications in the automotive industry. Here, a
number of welding tasks have been introduced along the car body: the goal is
to let the robots perform such tasks while minimizing the cycle time (or
_makespan_). The main difficulties, from the manufacturing engineer
perspective, lie in assigning the tasks to the robots, deciding the order and
the timing of the operations, avoiding collisions between the robots and the
environment, and among the robots themselves. We present in this work an
iterative approach, consisting of two steps: first, sequences for the robot
operations are computed in order to minimize the cycle time, while neglecting
collisions among robots; then, given the assignment of tasks to robots, the
operations are reordered and scheduled while avoiding conflicts among robots.
Robot motions are also automatically computed to avoid collisions with the
static environment. We show an optimal algorithm, for the first part, based on
implicit enumeration (branch and bound) and introduce a novel suboptimal
algorithm, for the second part, to synchronize the robots. These algorithms
are iterated while fetching information about the problem that are hard to
compute, thus following a lazy approach. Tests on problems adapted from the
literature and from the automotive industry, show clear improvements over more
sequential approaches and good running times. The algorithms are especially
suited for cases with up to 40 tasks and 4 robots, as for typical _geometry
stations_. In future works, we will further investigate efficient heuristic
optimization approaches in order to handle larger problems, consisting of more
than 100 tasks, as for typical _assembly stations_.

###### Index Terms:

Multirobot systems, Robot Scheduling, Path planning, Computer aided
manufacturing

##  I Introduction

In the last decades, for the automotive manufacturers short production times,
on one hand, and large volumes, on the other, have become crucial. In order to
achieve that, efficient equipment utilization is needed. Nevertheless,
resource optimization plays even an important role when considering
sustainable production systems, both economically and for the ecological
aspects, in terms of less energy consumption and working space, see [1].

In this work, we restrict the focus to the automotive manufacturing process.
The car body, also called ”Body in White” (BiW), is usually assembled together
on a line where several stations are placed serially and in parallel. At each
station multiple robots perform operations like stud, spot welding, and
sealing, and often share the same workspace, see Fig. 1.

In this framework we are mainly concerned with the problem of maximizing the
number of products assembled in a line. This can be directly translated into
minimizing the time needed at each station to perform predefined tasks in a
collision-free way. If a manufacturing line reduces the cycle time by $r$
($0<r<1$), then it would potentially increase the number of assembled products
by $r/(1-r)$. For a car manufacturer, this strong relation means that reducing
the cycle time, for example, by 33% would directly translate into 50% more
cars produced. The immediate impact is tremendous.

A very powerful way to shorten ramp up times and increase the volumes is to
simulate as soon as possible all the upcoming processes, including the design
and the manufacturing ones. Simulation software is therefore strongly needed
and provides the main key to virtual manufacturing. Unfortunately, nowadays,
there is still much work done manually inside the software tools, which is
both time consuming and error prone. On the other side, few or sometimes no
automatic software tools at all are available, due to the difficulties in
modeling the process and to the complexity of synthesizing solutions.

Here, we describe an algorithm and an automatic tool, aiming at planning robot
motions in order to achieve near-optimal programs, performing all predefined
tasks in the station, avoiding collisions between them and the environment,
and among the robots themselves.

Figure 1: A car body assembly line equipped with stud welding robots, modeled
in the simulation software IPS, Industrial Path Solutions, [2].

##  II Problem description and related literature

The challenge is to concurrently:

  * •

distribute the tasks among robots such that each of them is assigned to one
robot (set partitioning, dispatching),

  * •

decide in which way each robot should perform a task among several
alternatives,

  * •

find a sequence of tasks for each robot (routing),

  * •

compute robot paths that are collision-free w.r.t the static environment (path
planning),

  * •

schedule the sequences of paths such that no collision occurs among the robots
(scheduling).

The problem in its entire complexity has received attention only in the last
five years, at the best of the authors’ knowledge.

A reason for that might be that it is very interdisciplinary, integrating
together combinatorial optimization and path planning issues. Moreover, only
in the last years, with the development of computers and software tools, it
has been possible to face large problems in a reasonable time. Anyway, a wide
amount of literature has been produced in each of these two fields. There are
several works regarding the Vehicle Routing Problem (VRP), which is a more
general case of the Multiple Traveling Salesman Problem (MTSP); for a good
review of the VRP and its variations see [3]. Furthermore, several articles
have dealt with path planning for mobile and industrial robots; for a good
treatment of the subject see [4] and [5]. Here, we refer to the works that
most closely relate to ours, combining together the sub-problems described
above.

From a combinatorial optimization perspective, if we disregard collisions
among robots, the problem presented in this work can be modeled as a min-max
Generalized Multiple Traveling Salesman Problem (min-max MGTSP), or as an
uncapacitated min-max Generalized Vehicle Routing Problem, (min-max GVRP),
with multidepots. In contrast with the classical VRP, in the generalized one,
GVRP, customers are clustered in groups, and it is enough to visit one
customer in each group. The other difference is that the objective to minimize
is the length of the longest tour (min-max), not the sum of the tours lengths
(min-sum), as in the more classical formulation.

Little work has been produced in this area, compared to the wide literature
referring to the min-sum MTSP and VRP.

An exact algorithm for the min-max VRP, able to prove the optimality of a
solution provided to the problem presented at the Whizzkids ’96 competition,
has been proposed by Applegate _et al._ in [6]. The algorithm is based on
branch-and-cut and requires a very lare computing time on a computer network.
Some heuristic approaches have also been investigated, see [7]. In this work,
the author also adapts into No Depot min-max MTSP, and solves some instances
from the TSPLIB, see [8], thus giving results to publicly available benchmark
instances. In [9] a load balancing algorithm for min-max VRP is provided. It
runs by iterating the solution of a linear programming sub-problem and of a
TSP, to achieve the minimization of the longest tour, without explicitly
considering the min-max as objective function. There, theoretical lower bounds
are provided for problem instances in the 2D Euclidean space, based on
geometrical space partitions. The GVRP with min-sum objective has also been
studied. Ghiani _et al._ in [10] introduced an efficient transformation into a
Capacitated Arc Routing Problem (CARP). Anyway, we are mostly interested in
the min-max GVRP, where there is a lack of work.

All the articles so far cited do not deal at all with resource allocation
problems such as conflicts between the tours. By conflict we mean that parts
of the tours may not occur at the same time. These types of conflict are often
introduced as a way of modeling physical collisions among moving objects,
_e.g._ industrial/mobile robots or ground/air vehicles. Only in the last years
the problem has been treated. In [11] a min-max MTSP with conflicts is solved
by using a genetic algorithm and local search heuristics. In [12] the problem
is generalized to a min-max MGTSP with conflicts and solved exploiting genetic
algorithm as well, by using the solution of a GTSP as local search. Two works
trying to exploit the geometrical properties of the problems are [13] and
[14]. Here, the assignment of welding tasks to robots exploits their
geometrical distribution, in order to achieve a good trade-off between:
separating the tasks among robots to avoid collisions vs. partitioning them to
balance the robot loads and minimize the longest tour. The method is applied
to a complete line composed of three stations with a total of ten robots. In
[15], the Welding Cell Problem (WCP) is introduced, which is a variant of the
min-max MTSP with conflicts. Here, the problem is to find routes for robots
such that the tour time for each robot is within a predefined maximum time.
The authors discuss also the min-max MTSP. The approach uses a branch and
bound (B&B) method together with column generation. A related problem, the
Laser Source Problem, is studied in [16]. It is similar to the min-max MTSP
with conflicts, but with an additional requirement, which is to find the
minimum number of laser sources. The motivation is the same, a multirobot
station for car body manufacturing, and their solving approach is to use NP-
hard sub-problems in a B&B framework. Conflicts are handled as well and
results for problem instances with up to 40 jobs are illustrated. High
computing times are however needed for large instances and no path planning
capability to avoid conflicts is included.

Another similar problem is the one of routing of automated guided vehicles
(AGVs) and supervision of automated manufacturing systems (AMSs), where a
number of unmanned vehicles operate in a factory, or warehouse environment
transporting goods between various places along a network of uni- or bi-
directional lanes. A good overview of concepts and methods for AGV routing and
scheduling is provided by [17] and [18]. The main differences between our
problem and the AGV scheduling and routing, in general, are the absence, in
our settings, of a predefined roadmap network and the presence of several
alternatives that perform the tasks.

A large amount of literature concerning AGVs and AMSs is related to discrete
event modeling and different extensions of Petri nets are used, see [19, 20,
21, 22, 23, 24, 25, 26, 27, 28], for example. In these works, often, the main
focus is on modeling and on synthesizing a maximally permissive supervisor or
to enforce and guarantee _liveness_ of the system. Liveness is a property
implying that no _deadlock_ occurs. Deadlock is a very separated concept from
the one of _conflict_ used here. A conflict occurs when two agents or jobs
want to use the same resource or machine: it is explicitly specified within
the problem definition. A deadlock is a state requiring four conditions:
mutual exclusion, no preemption, hold while wait and circular wait. It is not
necessarily identified at the problem definition level, since its detection
often requires deep analysis of the system. Note that there exist problems
such as the classical job shop scheduling problem (JSSP) that

  * •

do contain conflicts, defined by the fact that two operations may not be done
at the same time on the same machine, and

  * •

do not generate any deadlock in their open-loop behavior since the ”hold while
wait” condition is not present as stated in [30].

The papers [31, 32] also take a look at alternative graph-based approaches.

In [33, 34, 35], various optimization is present, but with different
objectives and/or conditions than in our case. In the first case, the problem
consists in routing AGVs with only one fixed start and one fixed goal in a
predefined network. In the second case, robot energy is minimized with
collision avoidance as a side-constraints, and in the third case total
transport time is minimized.

The conflict-free routing problem, given a fixed scheduling, is addressed in
[36] utilizing a Dijkstra-like procedure on a time-window graph. In [37], a
column generation method is used to minimize the makespan given a fixed set of
assigned jobs. In [38] a more dynamic approach for conflict resolution is
suggested, computing AGV routes incrementally and thereby enabling conflict
avoidance under changing circumstances using a modified Banker’s algorithm. In
[39] the problem is decomposed into a master problem for the task assignment
(dispatching) and for routing, and a sub-problem for vehicle routing. It is
solved by subdividing the time horizon into time slots and applying MIP
algorithms. The network on which vehicles can move is precomputed. This
characteristic is also present in [40] and [41] where Constraint Programming
and MIP techniques are used. In [41], the assignment of operations to machines
is already determined together with the order in which operations will be done
by several jobs. The problem is very close to job shop scheduling. In [40],
results for different types of problems show that for _compact or spread out_
requests sets there is a natural decoupling between the scheduling and routing
problems.

That property, together with the geometrical characteristics about the
problems addressed here, motivate the strategy to consider the two extremes of
a scale. On the one hand, we try to generate lower bounds by assuming that no
collision between the robot paths is present. To do that, we use a B&B
approach, naturally branching on the tasks to be assigned, as in [16], but we
generalize it to solve the min-max MGTSP. Moreover, we present two different
ways to improve the lower bounds, thus decreasing the number of expanded nodes
and running times. On the other hand, we want to generate feasible solutions
even for cluttered environments, where a lot of collisions may occur, by
preferring as quality measure the feasibility of the solution to the cycle
time minimization. To achieve that, we introduce here a novel algorithm that
solves a specialized synchronous routing and scheduling problem, in which the
assignment of tasks to robots is fixed and robots are assumed to move
synchronously. The algorithm proposed is based on a new transformation of the
problem into a GTSP. This may be interpreted as an attempt to redesign robot
paths by means of combinatorial optimization. In the rest of this article we
will use the terms _robots_ and _agents_ indifferently.

##  III Graph formulation and notation

The problem can be represented in graph terms. Given $N_{A}$ agents
$a_{1},a_{2},...,a_{N_{A}}$ and a set
$G=\\{g_{1},\ldots,g_{N_{A}},\ldots,g_{N_{G}}\\}$ of $N_{G}$ number of tasks,
all possible ways that agent $a_{k}$ can perform the tasks in $G$ may be
represented by a set $V^{k}$. Each agent has also a home position, which may
be achieved by several configurations. We have included these home tasks into
the set of tasks $G$ and denoted them $g_{1},\ldots,g_{N_{A}}$, without loss
of generality. Thus, each $V^{k}$ can be partitioned into groups (sometimes
called clusters in the literature) and the information about which elements of
$V^{k}$ perform a specific task is given by the function
$\gamma^{k}:V^{k}\rightarrow G$, which maps each element into a task in $G$.
In the rest of the paper we will denote with $D^{k}_{i}$ the set of vertices
in $V^{k}$ that can perform task $g_{i}$, _i.e._ $D^{k}_{i}=\\{v\in
V^{k}:\gamma^{k}(v)=g_{i}\\}$. A set of edges $E^{k}\subseteq V^{k}\times
V^{k}$ can be naturally added to the set of vertices, thus forming the graph
$(V^{k},E^{k})$. This graph can be extended into a weighted one by the help of
a non-negative weight (or cost) function
$c^{k}:E^{k}\rightarrow\mathbb{R^{+}}$. Since we assume this cost function to
be symmetric, we will talk about arcs, instead of edges. The weight function
models the time needed for agent $a_{k}$ to move from a start configuration
$u\in V^{k}$ to an end configuration $v\in V^{k}$, with $u\neq v$, by
$c^{k}(u,v)$. As a special case the time needed to perform task
$\gamma^{k}(v)$ by configuration $v$ is modeled as $c^{k}(v,v)$: we will, in
the rest of the paper, indicate this processing time with $c^{k}(v)$ and omit
the $k$ since the $c^{k}(u,v)$ is completely determined by the vertices $u,v$.
At this point we may define a tour $T^{k}$ for agent $a_{k}$ as a sequence of
vertices, starting and ending at a vertex belonging to the corresponding home
task $g_{k}$: $T^{k}=\\{u^{k}_{1},\ldots,u^{k}_{N_{T^{k}}},u^{k}_{1}\\}$,
where $u^{k}_{i}\in V^{k},\forall i=1,\ldots,N_{T^{k}}$, and such that
$\gamma^{k}(u^{k}_{i})\neq\gamma^{k}(u^{k}_{j}),\forall i\not=j$ and
$\gamma^{k}(u^{k}_{1})=g_{k}$. The cost associated to a tour $T^{k}$ is given
by

$c(T^{k})=\sum\limits_{i=1}^{N_{T^{k}}-1}\left[c(u^{k}_{i})+c(u^{k}_{i},u^{k}_{i+1})\right]+c(u^{k}_{N_{T}^{k}})+c(u_{N_{T^{k}}},u^{k}_{1})$
(1)

Given a subset of tasks $\tilde{G}\subseteq G$, one can write the induced GTSP
as

minimize $c(T^{k})$ _s.t._ $\displaystyle\forall g\in\tilde{G},\exists u\in
T^{k}:\gamma^{k}(u)=g$ (2a)

The best value obtained is refered as $GTSP(\tilde{G})$. A feasible solution
to the min-max MGTSP is an ordered set of tours $T=(T^{1},\ldots,T^{N_{A}})$,
such that the corresponding sets $(G^{1},\ldots,G^{N_{A}})$ form a partition
of $G$, _i.e._ $G^{k}\cap G^{l}=\emptyset,\forall k\neq l$ and
$\bigcup\limits_{k=1}^{N_{A}}G^{k}=G$. The cost associated to $T$ is given by

$c(T)=\max\limits_{k=1,\ldots,N_{A}}\\{c(T^{k})\\}$ (3)

The min-max MGTSP with _conflicts_ has to take also care of incompatible arcs
between the tours. A _conflict_ between two arcs $(s^{k}_{i},s^{k}_{j})$ and
$(s^{l}_{m},s^{l}_{n})$ exists if they are defined to be incompatible by the
modeled process. This is often introduced in order to model geometrical
collisions between robot paths or paths where the distance between the robots
is under the minimum allowed clearance threshold.

We can define a scheduled tour as a tour with time information, that is a re-
parametrization of the tour itself

$T^{k}_{scheduled}=\left(t[s^{k}_{1}]^{L},\ldots,t[s^{k}_{i}]^{A},t[s^{k}_{i}]^{L},\ldots,t[s^{k}_{1}]^{A}\right)$
(4)

where, $t[s^{k}_{i}]^{A}$ is the time agent arrives at $s^{k}_{i}$ and
$t[s^{k}_{i}]^{L}$ is the time the agent leaves the same point. Given
scheduled tours, two conflicting arcs $(s^{k}_{i},s^{k}_{i+1})$ in $T^{k}$ and
$(s^{l}_{m},s^{l}_{m+1})$ in $T^{l}$, with $k\neq l$, define an _active
conflict_ if
$\left[t[s^{k}_{i}]^{L},t[s^{k}_{i+1}]^{A}\right]\cap\left[t[s^{l}_{m}]^{L},t[s^{l}_{m+1}]^{A}\right]\neq\emptyset$,
_i.e._ if they overlap in time. At the same way, there is an active conflict
between a vertex $s^{k}_{i}$ and arc $(s^{l}_{m},s^{l}_{m+1})$ if
$\left[t[s^{k}_{i}]^{A},t[s^{k}_{i}]^{L}\right]\cap\left[t[s^{l}_{m}]^{L},t[s^{l}_{m+1}]^{A}\right]\neq\emptyset$,
and between vertices $s^{k}_{i}$ and $s^{l}_{m}$ if
$\left[t[s^{k}_{i}]^{A},t[s^{k}_{i}]^{L}\right]\cap\left[t[s^{l}_{m}]^{A},t[s^{l}_{m}]^{L}\right]\neq\emptyset$.
Note that when we identify the cost with time, then we have
$c(T^{k})=t[s^{k}_{1}]^{A}$.

An optimal solution to the min-max MGTSP with conflicts will consist in an
ordered set of scheduled tours, with minimum makespan, such that no active
conflict is present.

##  IV Optimal solution for min-max MGTSP without conflicts

Achieving optimality in absence of conflicts is very important in terms of
cycle time, because the solution obtained:

  * •

can result to be collision-free due to the geometrical properties of the
problems;

  * •

can be a very good initial guess for improvement-based heuristic algorithms;

  * •

can be used to dispatch the tasks among the agents and then resolving
conflicts by avoiding the possibility to redistribute tasks, with the aim to
decrease complexity;

  * •

can be velocity tuned to avoid conflicts and still be of very high quality;

  * •

from a practical perspective, can be modified to be collision-free by small
manual adjustments.

Therefore, in order to find these potential best solutions, we relax the
problem by assuming that there are no conflicts among the agents arcs. This
relaxation transforms the whole problem into a min-max MGTSP, _i.e_ without
conflicts. The problem instances faced in this work contain a small number of
tasks, under 40, and can well model some real world applications like
_geometry stations_. For larger problems, more sophisticated algorithms have
to be provided.

In this section, we discuss an optimal algorithm and show how this can be
easily implemented. One practical advantage is the ease of implementation, due
to the absence of bounds based on LP (Linear Programming) formulations, which
often require sophisticated implementations exploiting sparsity and dealing
with ill conditioned instances.

###  IV-A Branch and bound based algorithm

The general branch and bound (B&B) architecture in this work follows the lazy
pattern described in [42]. Each sub-problem consists of a GTSP: they are
solved to optimality, due to the limited size of the instances, by a dynamic
programming approach, which is a generalization of [43].

Along the search in the B&B tree, at each node, we maintain two sets, $G_{S}$
and $G_{N}$ that partition $G$. $G_{S}$ contains the tasks assigned to some of
the agents, whereas the remaining tasks, the ones not yet assigned to anyone,
are the elements of $G_{N}$. $G_{S}$ can be further partitioned into $N_{A}$
subsets indicated by $G^{k}_{S}$, with $k=1,\ldots,N_{A}$: each set
$G^{k}_{S}$ represents the tasks in $G_{S}$ assigned to agent $a_{k}$.

The _selection_ of the next node to be expanded in the B&B may be done by
depth-first, breadth-first, best-first strategies or by some other criterion.
Our experience has shown that depth-first performs well in terms of computing
times and has the enormous advantage of keeping the dimension of the nodes
queue limited. This avoids state space explosion.

The _branching_ step consists in assigning, if possible, an element $g_{i}\in
G_{N}$ to each of the agents, thus creating at most $N_{A}$ children nodes.
The group to be assigned can be randomly chosen in $G_{N}$, or decided based
on some heuristics that exploit information specific to the problem.

Algorithm 1 Creating children nodes in the branching step of B&B

$G_{N}\leftarrow G_{N}\setminus\\{g_{i}\\}$

for $k=1$ to $N_{A}$ do

if $D^{k}_{i}\neq\emptyset$ then

NEWNODE($G^{1}_{S},\ldots,G^{k}_{S}\cup\\{g_{i}\\},\ldots,G^{N_{A}}_{S},G_{N}$)

end if

end for

The bounding step computes a _lower bound_ for the current node. The first
part of the bound is solving $N_{A}$ GTSPs, one for each agent. The lower
bound $\underline{b}$ is then

$\displaystyle\underline{b}^{k}(G^{k}_{S})=$ $\displaystyle GTSP(G^{k}_{S})$
$\displaystyle\forall k=1,\ldots,N_{A}$ (5) $\displaystyle\underline{b}=$
$\displaystyle\max\limits_{k=1,\ldots,N_{A}}\\{\underline{b}^{k}(G^{k}_{S})\\}$

Note that if the triangle inequality holds among the vertices of the problem,
then it is possible to guarantee that the value obtained is actually a lower
bound on the optimal value, see Section 1.

When reaching a leaf in the tree, _i.e._ when
$G_{S}=\bigcup\limits_{k=1}^{N_{A}}G^{k}_{S}=G$, the bound obtained is equal
to an _upper bound_ $\overline{b}$ for the entire problem, _i.e._ the min-max
MGTSP, without conflicts. If the bound $\overline{b}$ computed is less than
the best bound $\overline{b}_{MIN}$ so far found, we set
$\overline{b}_{MIN}=\overline{b}$. This method will be referred as Approach 1.

Anyway, we can further improve the algorithm performance by using additional
knowledge about the process. Indeed, modeling operations like welding times,
on one hand, and modeling robot reachability constraints, on the other, can be
utilized to improve the lower bounds.

###  IV-B Tightening the bound

Tasks may consist in closing, welding and opening a gun, sealing or other
operations. The time it takes for the accomplishment of one such task is
modeled to be constant, namely $c_{G}$. This is due to the fact that, often,
the actual duration of such operations is considered to not heavily vary among
different robots and operations. The engineers responsible for the simulations
insert this value as input to the algorithms described in this article. This
characteristic may be exploited to improve the branching step and to tighten
the lower bounds. To do that, we use the following lemma:

###### Lemma 1

Given a best tour $T^{*}_{n}$ with length $c(T^{*}_{n})$ for an SGTSP $G_{n}$
with $n$ groups, and an SGTSP $G_{n+1}=G_{n}\cup\\{g_{n+1}\\}$, if the
triangle inequality holds, then

$c(T^{*}_{n+1})\geq\max\\{c(T^{*}_{n})+2\min(g_{n+1},G_{n})-\max(G_{n}),c(T^{*}_{n})\\}$
(6)

where $c(T^{*}_{n+1})$ is the length of an optimal tour in $G_{n+1}$,

$\displaystyle\min(g_{n+1},G_{n})$ $\displaystyle=$
$\displaystyle\min\limits_{\begin{subarray}{c}p_{i}\in G_{n}\\\ p_{k}\in
g_{n+1}\end{subarray}}c(p_{i},p_{k})$ $\displaystyle\max(G_{n})$
$\displaystyle=$ $\displaystyle\max\limits_{\begin{subarray}{c}p_{i}\neq
p_{j}\in G_{n}\end{subarray}}c(p_{i},p_{j})$

For a proof of the above Lemma, see Appendix A. The idea with the first part
is to avoid generating B&B nodes where vertices far away from small clusters
are added. The second part of the inequality, $c(T_{n+1}^{*})\geq
c(T_{n}^{*})$, is just a motivation for why the B&B is consistent.

Algorithm 2 Revised branching step exploiting geometrical information

$G_{N}\leftarrow G_{N}\setminus\\{g_{i}\\}$

for $k=1$ to $N_{A}$ do

$\delta=2\min(g_{i},G^{k}_{S})-\max(G^{k}_{S}$)

if $D^{k}_{i}\neq\emptyset$ and
$\max\\{\underline{b}^{k}(G^{k}_{S}),\underline{b}^{k}(G^{k}_{S})+\delta\\}$
$<\overline{b}_{MIN}$ then

NEWNODE($G^{1}_{S},\ldots,G^{k}_{S}\cup\\{g_{i}\\},\ldots,G^{N_{A}}_{S},G_{N}$)

end if

end for

Nonetheless, also the computation of the lower bound can be improved. In fact,
at each node in the B&B tree, the remaining tasks, _i.e._ the ones represented
as elements of $G_{N}$, have to be performed eventually by some agent. In the
classical formulation of the GVRP and of the MGTSP, all vertices can be
reached by all agents (vehicles or salesmen). In a typical multirobot station,
however, this is not always true. Tasks, for example, may be outside the
working space of some robots or collisions are not avoidable.

That leads to the following model, where the values obtained by the solution
of the GTSPs are increased with the minimal value it will take to cover all
remaining tasks.

minimize $c$ _s.t._ $\displaystyle c-c_{G}\sum\limits_{i:g_{i}\in
G_{N}}x_{ik}$ $\displaystyle\geq\underline{b}^{k}(G^{k}_{S})$
$\displaystyle\forall k=1,\ldots,N_{A}$ (7a)
$\displaystyle\sum\limits_{k=1}^{N_{A}}x_{ik}$ $\displaystyle=1$
$\displaystyle\forall i:g_{i}\in G_{N}$ (7b) $\displaystyle x_{ik}$
$\displaystyle=0$ $\displaystyle\forall i,k:D^{k}_{i}=\emptyset$ (7c)
$\displaystyle x_{ik}\in$ $\displaystyle\\{0,1\\}$ $\displaystyle\forall
i:g_{i}\in G_{N}$ $\displaystyle\forall k=1,\ldots,N_{A}$

The cycle time is indicated with $c$. The unknowns $x_{ik}$ take the value 1
if task $g_{i}$ is assigned to agent $a_{k}$, take the value 0 otherwise.
Constraints (7a) state that the cycle time should be greater than or equal to
the sum of the best tour among the groups in $G^{k}_{S}$, plus the sum of the
processing times for the task assigned to agent $a_{k}$. Equalities (7b)
establish that each task should be performed by exactly one agent. Constraints
(7c) model that task $g_{i}$ is unreachable for agent $a_{k}$.

Problem 7 can be solved by general purpose Mixed Integer Linear Programming
(MILP) packages. This method will be referred as Approach 2.

Anyway, in some cases, these packages do not return any optimal solution in a
reasonable time. Therefore, we relax problem (7) by eliminating (7c), meaning
that all agents can perform all tasks. The new problem becomes:

minimize $c$ _s.t._ $\displaystyle c-c_{G}n_{k}$
$\displaystyle\geq\underline{b}^{k}(G^{k}_{S})$ $\displaystyle\forall
k=1,\ldots,N_{A}$ (8a) $\displaystyle\sum\limits_{k=1}^{N_{A}}n_{k}$
$\displaystyle=|G_{N}|$ (8b) $\displaystyle n_{k}\in$
$\displaystyle\mathbb{N}$ $\displaystyle\forall k=1,\ldots,N_{A}$

In this formulation the unknowns $n_{k}$ are the number of tasks assigned to
agent $a_{k}$, which relates to the variables in (7) by
$n_{k}=\sum\limits_{i:g_{i}\in G_{N}}x_{ik}$. The problem can be optimally
solved in polynomial time by a greedy algorithm, _i.e._ by assigning one task
at the time to the agent with the shortest tour. For a closed form expression,
see Appendix A. This method will be referred as Approach 3.

###  IV-C Simulation Results

TABLE I: Comparison of B&B properties among different lower bounds used, on geometric instances from TSPLIB with four agents. |  | Approach 1 |  | Approach 2 |  | Approach 3 |
---|---|---|---|---|---|---|---
Problem name | Best solution | N. nodes | Time (s) | N. nodes | Time (s) | N. nodes | Time (s)
11eil51 | 72.804 | 60 | 0.000 | 47 | 0.015 | 46 | 0.093
14st70 | 190.036 | 425 | 0.015 | 360 | 0.016 | 345 | 0.468
16eil76 | 131.724 | 2912 | 0.047 | 2410 | 0.047 | 2325 | 6.334
16pr76 | 34001.008 | 360 | 0.031 | 263 | 0.031 | 258 | 0.515
20kroA100 | 7038.138 | 2925 | 0.094 | 2281 | 0.094 | 2137 | 6.115
20kroB100 | 6545.117 | 3533 | 0.156 | 2950 | 0.156 | 2471 | 4.977
20kroC100 | 5071.016 | 4499 | 0.124 | 3208 | 0.094 | 2679 | 12.792
20kroD100 | 6357.273 | 12919 | 0.250 | 9718 | 0.249 | 9022 | 18.720
20kroE100 | 6751.198 | 2263 | 0.078 | 1816 | 0.078 | 1766 | 5.273
20rat99 | 292.655 | 7866 | 0.187 | 6048 | 0.172 | 5476 | 13.759
21eil101 | 178.327 | 317570 | 4.524 | 266014 | 4.352 | 231098 | 476.115
21lin105 | 4622.421 | 7960 | 0.203 | 6499 | 0.219 | 6372 | 12.698
22pr107 | 15698.472 | 219 | 0.078 | 110 | 0.047 | 110 | 0.172
25pr124 | 22462.762 | 23355 | 0.530 | 16677 | 0.468 | 15022 | 24.898
26bier127 | 32507.758 | 6885 | 0.297 | 4144 | 0.265 | 3749 | 14.196
26ch130 | 1578.800 | 26701 | 0.561 | 10799 | 0.343 | 10615 | 32.448
28pr136 | 25694.316 | 109615 | 2.886 | 48242 | 1.622 | 47049 | 257.511
29pr144 | 22378.426 | 148788 | 3.557 | 103607 | 3.339 | 103024 | 95.129
30kroA150 | 7166.586 | 153316 | 4.181 | 79230 | 3.151 | 77698 | 609.933
30kroB150 | 7309.708 | 425300 | 11.809 | 241628 | 7.051 | 232518 | 1223.001
31pr152 | 28024.425 | 29786 | 6.177 | 19071 | 3.479 | 17303 | 55.318
32u159 | 13413.603 | 1607531 | 31.544 | 715318 | 17.301 | 706262 | 1796.913
39rat195 | 645.101 | 49179533 | 2184.513 | 4535693 | 335.761 | 4298503 | 44344.205
40d198 | 4779.423 | 1076643 | 1946.830 | 385042 | 824.200 | 371681 | 1599.696
40kroA200 | 8106.633 | 27993063 | 1012.104 | 3860702 | 185.283 | 3820629 | 24813.831

Here, we show some results from running this algorithm on modified SGTSP
instances from the TSPLIB, see [8], adapted to our problem. When it comes to
the MILP part to solve (7), we have used the COIN-OR package, [44].

Table I shows the results of the computations on geometrical instances of the
problem, based on the Euclidean distance. The table compares 3 different
bounds:

  1. 1.

approach 1 uses (5) to bound the function;

  2. 2.

approach 2 improves it by using Lemma 1, Algorithm 2, and the solution to
problem (7);

  3. 3.

approach 3 replaces the solution to problem (7) with the one to problem (8);

As expected, the number of expanded nodes decreases from approach 1, towards
approach 3 and approach 2, and so almost do the running times. Note that, for
the largest instance, computing times can almost be 10 times faster than with
the first approach.

Problem (7) has too large memory requirements when solved by COIN-OR, thus we
use it only if the number of tasks is less than 10 (using approach 3
otherwise). Even in that case, there are large computing times that do not
justify its use with respect to the greedy bound computed in Approach 3.

We have also run simulations where the processing time $c_{G}$ at each task
was decreased. We have noticed that, by decreasing $c_{G}$, the advantages of
approach 2 are not as evident, whereas the bound computed as in 2 gives some
little improvements.

Experiments by changing the initial solution have also proven its enormous
effect on the dimension of the explored state space.

These characteristics motivate, on the one side, the study for faster solution
to problem (7), and, on the other, the use of other neighborhoods and meta-
heuristics for the search, to improve initial solutions.

At the end of the B&B, an optimal solution with value $\overline{b}_{MIN}$ is
obtained, consisting of a tour for each agent
$T^{k}=(s^{k}_{1},\ldots,s^{k}_{N_{T_{k}}},s^{k}_{1})$, where the tasks have
been partitioned into $(G^{1},\ldots,G^{N_{A}})$, with
$G^{k}=\bigcup\limits_{i=1}^{N_{T^{k}}}s^{k}_{i}$.

##  V Synchronous routing and scheduling

Here, we fix the load for each agent to $(G^{1},\ldots,G^{N_{A}})$ and let
each $G^{k}$ be composed of exactly one element per group. This partition may
be the result of the B&B algorithm above, in which case each
$G^{k}=\\{s^{k}_{1},\ldots,s^{k}_{N_{T^{k}}}\\}$, or is due to a process
specific fixed assignment. On the other hand, reordering the groups within one
agent, and tuning the velocities are allowed. The goal is now to find a
sequence of vertices for each agent such that no collision occurs.

In order to cope with the high complexity of the problem, we introduce the
constraint that the agents move synchronously. This means that they will move
from one task to another, starting and finishing simultaneously. As special
case one or more agents can remain still, while others might move. This
assumption also finds concrete applications such as the one in [46], and in
the latest controllers family from some robot manufacturer, see ABB MultiMove
[49].

A natural way of modeling the problem is by using a synchronized state space
$S_{G}$, consisting of the Cartesian product of the sets $G^{k}$, _i.e._
$S_{G}=G^{1}\times\ldots\times G^{N_{A}}$.

This problem is not very common in the literature. The closest work found is
[46]. There, inspired by laser drilling applications for microelectronics
manufacturing system, the authors solve the coverage of planar points by
multiple robots, that is routing and scheduling with collision avoidance
constraints. Since the application has thousands of points, the problem is
divided into two sub-problems: a so called splitting problem and an ordering
problem, which are solved by heuristics.

When dropping the synchronous assumption, the problem arisen is the so called
laser sharing problem, see [45], where the assignment of jobs to robots is
fixed, and is indicated as ”RSP-J” in the same paper. The problem is solved by
a B&B that generates a large number of sequences.

The solution obtained in this way can always be improved by relaxing back to
the original problem, fixing the tours, and using a velocity tuning algorithm,
see Section V-C. We propose here solving this problem by creating a GTSP. The
original problem is re-modeled with the aim to apply a known solving algorithm
for the transformed problem.

Each group in the transformed model will be visited once, but the original
vertices might be reached more than once: this can facilitate the solution of
some problems.

###  V-A Two agents case

We now want to find a sequence of states in $S_{G}=G^{1}\times G^{2}$
$=\\{(s^{1}_{1},s^{2}_{1}),\ldots,(s^{1}_{m},s^{2}_{1}),\ldots,(s^{1}_{m},s^{2}_{n}),\ldots,(s^{1}_{N_{T_{1}}},s^{2}_{N_{T_{2}}})\\}$,
such that, when projecting a state onto $G^{1}$, respectively $G^{2}$, all
their elements are visited. In other words, the goal is to visit all tasks
modeled as element of $G^{1}$, resp. $G^{2}$.

$\displaystyle
G^{1}=\\{s_{1}^{1},s_{2}^{1},\ldots,s_{i}^{1},\ldots,s_{m}^{1}\ldots,s^{1}_{N_{T_{1}}}\\}$
$\displaystyle
G^{2}=\\{s_{1}^{2},s_{2}^{2},\ldots,s_{j}^{2},\ldots,s_{n}^{2},\ldots,s^{2}_{N_{T_{2}}}\\}$
(9)

Each state $s\in S_{G}$ can be uniquely identified by its two components
$s_{i}$ in $G^{1}$ and $s_{j}$ in $G^{2}$, and be referred as $s_{ij}$:

$S_{G}=\\{s_{11},\ldots,s_{i1},\ldots,s_{ij},\ldots,s_{N_{T_{1}}N_{T_{2}}}\\}$
(10)

The set $S_{G}$ is then provided with a set of arcs between its vertices, thus
forming a graph $(S_{G},S_{G}\times S_{G})$. Since the agents move
synchronously, the cost (time) for the arc $(s_{ij},s_{mn})$ is given by

$c(s_{ij},s_{mn})=\max\\{c(s^{1}_{i},s^{1}_{m}),c(s^{2}_{j},s^{2}_{n})\\}$
(11)

The objective is to find a sequence that minimizes cycle time while avoiding
conflicting states.

In order to obtain such a sequence, we present here a new transformation into
a well known problem, the GTSP, for which efficient algorithms are known. The
solution to the GTSP, then, can be always transformed back into a solution to
the original problem. To achieve that we generate a second layer $S_{G}^{2}$
that is a copy of the first one, $S_{G}$ (from now on referred as
$S_{G}^{1}$). The vertices are grouped together in the following way to create
a GTSP: vertices in $S_{G}^{1}$ are clustered in ”vertical” groups, one for
each column, whereas vertices in $S_{G}^{2}$ are clustered in ”horizontal”
groups, one for each row, see Fig. 2. More formally, the first layer consists
of $N_{T^{1}}$ groups $F_{i},\forall i=1,\ldots,N_{T^{1}}$

$F_{i}=\\{(s^{1}_{i},s^{2}_{j})^{1},\forall j=1,\ldots,N_{T_{2}}\\}$ (12)

The second layer will have $N_{T^{2}}$ groups $F_{j+N_{T^{1}}},\forall
j=1,\ldots,N_{T^{2}}$

$F_{j+N_{T^{1}}}=\\{(s^{1}_{i},s^{2}_{j})^{2},\forall i=1,\ldots,N_{T^{1}}\\}$
(13)

The idea behind this grouping is to create a problem aiming at fulfilling the
constraints of visiting each original vertex at least once. Then, to minimize
cycle time avoiding collisions, arcs are added between these vertices, and
weighted according to (11)

$c\left((s^{1}_{i},s^{2}_{m})^{l_{1}},(s^{1}_{j},s^{2}_{n})^{l_{2}}\right)=c(s^{l_{1}}_{ij},s^{l_{2}}_{mn})=\\\
=c(s_{ij},s_{mn}),\forall l_{1},l_{2}\in\\{1,2\\}$ (14)

If there is a conflict between these vertices, then the cost will be $\infty$.
Thus, a GTSP has been defined consisting of $2N_{T_{1}}N_{T_{2}}$ vertices and
$N_{T_{1}}+N_{T_{2}}$ groups.
Resuming: one can create a GTSP by defining

  * •

a set of vertices $S^{1}_{G}\cup S^{2}_{G}=$

$\\{s^{1}_{11},\ldots,s^{1}_{mn},\ldots,s^{1}_{N_{T_{1}}N_{T_{2}}},s^{2}_{11},\ldots,s^{2}_{mn},\ldots,s^{2}_{N_{T_{1}}N_{T_{2}}}\\}$
(15)

  * •

the distances among them:

$c(s^{1}_{ij},s^{1}_{mn})=c(s^{1}_{ij},s^{2}_{mn})=c(s^{2}_{ij},s^{1}_{mn})=c(s^{2}_{ij},s^{2}_{mn})=\\\
=\max\\{c(s^{1}_{i},s^{1}_{m}),c(s^{2}_{j},s^{2}_{n})\\}$ (16)

  * •

the organization of vertices into groups:

$\displaystyle F_{1}$ $\displaystyle=$ $\displaystyle\\{s^{1}_{1j},\forall
j=1,\ldots,N_{T_{2}}\\}$ $\displaystyle F_{2}$ $\displaystyle=$
$\displaystyle\\{s^{1}_{2j},\forall j=1,\ldots,N_{T_{2}}\\}$
$\displaystyle\ldots$ $\displaystyle F_{N_{T^{1}}}$ $\displaystyle=$
$\displaystyle\\{s^{1}_{N_{T^{1}}j},\forall j=1,\ldots,N_{T_{2}}\\}$
$\displaystyle F_{1+N_{T^{1}}}$ $\displaystyle=$
$\displaystyle\\{s^{2}_{i1},\forall i=1,\ldots,N_{T^{1}}\\}$
$\displaystyle\ldots$ $\displaystyle F_{N_{T^{2}}+N_{T^{1}}}$ $\displaystyle=$
$\displaystyle\\{s^{2}_{iN_{T^{2}}},\forall i=1,\ldots,N_{T^{1}}\\}$

$(s^{1}_{1},s^{2}_{3})^{1}$$(s^{1}_{1},s^{2}_{2})^{1}$$(s^{1}_{1},s^{2}_{1})^{1}$$(s^{1}_{2},s^{2}_{3})^{1}$$(s^{1}_{2},s^{2}_{2})^{1}$$(s^{1}_{2},s^{2}_{1})^{1}$$(s^{1}_{1},s^{2}_{3})^{2}$$(s^{1}_{2},s^{2}_{3})^{2}$$(s^{1}_{1},s^{2}_{2})^{2}$$(s^{1}_{2},s^{2}_{2})^{2}$$(s^{1}_{1},s^{2}_{1})^{2}$$(s^{1}_{2},s^{2}_{1})^{2}$
Figure 2: The synchronized state space with its clone, and the ”vertical” and
the ”horizontal” groups.

Note that each group $F_{i}$ is visited exactly once, but vertices in $G_{1}$
and $G_{2}$ may occur more than once. This makes the problem closer to a path
planning problem, where agents search for a collision-free way to move.
Indeed, by allowing that, one may obtain better solutions for instances
modeling very cluttered environments. This is the main advantage of this
algorithm against other meta-heuristic search algorithms, which explore a more
limited state space. Note that a predefined roadmap where agents can move on
is not present here as for AGV applications.

Consider, _e.g._ Fig. 3: there are two agents with circular shape at their
home position $H_{1}$, and $H_{2}$ respectively. After the B&B phase they are
assigned, respectively vertex $(A)$, and $(B,C)$, _i.e._ $G_{1}=\\{H_{1},A\\}$
and $G_{2}=\\{H_{2},B,C\\}$. By a standard approach, since they collide when
moving from $H_{1}$ to $A$, resp. from $B$ to $C$, the best cycle time will be
6.8. At the contrary, by transforming the problem into a GTSP, one obtains the
solution in Fig. 4, which has a total cost of 6: $(H_{1},H_{2})$, $(H_{1},B)$,
$(A,H_{2})$, $(H_{1},C)$, $(H_{1},H_{2})$. Note that the problem has been
transformed into a symmetric one, therefore also the reversed sequence is an
optimal one.

$H_{1}$$H_{2}$$A$$B$$C$$2$$1$$1$$1.8$ Figure 3: An example with two agents
where the transformation to a GTSP finds a better solution by re-routing and
scheduling the agents: geometrical view.
$(H_{1},C)^{1}$$(H_{1},B)^{1}$$(H_{1},H_{2})^{1}$$(A,C)^{1}$$(A,B)^{1}$$(A,H_{2})^{1}$$(H_{1},C)^{2}$$(A,C)^{2}$$(H_{1},B)^{2}$$(A,B)^{2}$$(H_{1},H_{2})^{2}$$(A,H_{2})^{2}$
Figure 4: An optimal solution to the problem represented in Fig. 3.

The advantage of this transformation is that a vertex is allowed to be visited
more than once, allowing to obtain shorter tours in some cases. However, some
problems remain still unsolvable, even if a feasible solution exists, see
Appendix B. Note that shorter tours may be obtained because the arcs in
$S_{G}\times S_{G}$ do not satisfy the triangle inequality. They are not the
result of an optimal path planning algorithm that exploits all the degrees of
freedom of both robots simultaneously.

###  V-B N agents case

The generalization to the case of $N_{A}$ agents is straightforward and can be
done by building $N_{A}$ layers. Each layer $S_{G}^{k}$ consists of
$|S_{G}|=\prod\limits_{k=1}^{N_{A}}N_{T^{k}}$ vertices, grouped in $N_{T^{k}}$
groups, in the following way. Group $F^{k}_{i}$ in the $k$-th layer will be:

$\displaystyle F^{k}_{i}=$
$\displaystyle\\{(s^{1}_{j_{1}},\ldots,s^{k}_{j_{k}},\ldots,s^{N_{A}}_{j_{N_{A}}})$
$\displaystyle:j_{k}=k,$ $\displaystyle j_{i}=1,\ldots,N_{T_{k}},$
$\displaystyle\forall i=1,\ldots,N_{A},i\neq k\\}$ $\displaystyle\forall
k=1,\ldots,N_{A}$

Thus, the GTSP built has a total of $N_{A}|S^{G}|$ vertices distributed within
$\sum\limits_{k=1}^{N_{A}}N_{T^{k}}$ groups. Note that the number of vertices
grows exponentially with the number of agents, making the approach applicable
especially for problems with a few number of agents, and that is the case for
car manufacturing stations where often two to four industrial robots are
present. Another issue worth to mention is that a pure centralized path
planning approach would require generating motions for a generalized robot
whose number of degrees of freedom (dofs) is the sum of the dofs of each
single robot. Considering that a typical station has often four robots and
that each industrial robot usually has six dofs, then it is likely to think
that iterative and decoupled approaches avoiding collisions are still
motivated. The algorithm described above, indeed, requires path planning for
each robot alone, not considering the others.

In Fig. 5 three tours are illustrated after running the B&B algorithm and in
Fig. 6 the tours after the synchronization.

The GTSP solver used is not an exact one, but it has good performance on the
SGTSP instances from TSPLIB, see [47]. It computes a starting tour and tries
to improve it, by some local search techniques: some of them are similar to
the ones described in [48]. Anyway, any other GTSP solver algorithm may be
used, among the many present in the literature, see [48]. Near optimal
solvers, especially if proved to be very powerful, may also lead to high
quality global solutions.

###  V-C Asynchronous smoothing

The results obtained through the synchronous routing may be further reduced
when dropping the ”synchronous” assumption. By doing that, the agents are
allowed to move asynchronously while keeping the constraint of never being
simultaneously on conflicting arcs and keeping the ordered sequence of
vertices that have been assigned. The problem can therefore be modeled with a
classical MILP formulation:

minimize $c$ s.t. $\displaystyle c\geq$ $\displaystyle c(T_{k})$
$\displaystyle\forall k=1,\ldots,N_{A}$ (19a) $\displaystyle c(T_{k})\geq$
$\displaystyle t[s^{k}_{N_{T_{k}}}]^{L}+c(s^{k}_{N_{T_{k}}},s^{k}_{1})$
$\displaystyle\forall k=1,\ldots,N_{A}$ (19b) $\displaystyle
t[s^{k}_{1}]^{L}\geq$ $\displaystyle 0$ $\displaystyle\forall
k=1,\ldots,N_{A}$ (19c) $\displaystyle t[s^{k}_{i+1}]^{A}\geq$ $\displaystyle
t[s^{k}_{i}]^{L}+c(s^{k}_{i},s^{k}_{i+1})$ $\displaystyle\forall
i=1,\ldots,N_{T_{k}}-1$ $\displaystyle t[s^{k}_{i+1}]^{L}\geq$ $\displaystyle
t[s^{k}_{i+1}]^{A}$ $\displaystyle\forall i=1,\ldots,N_{T_{k}}-1$
$\displaystyle\forall k=1,\ldots,N_{A}$ (19d) $\displaystyle
t[s^{k}_{i}]^{L}\geq$ $\displaystyle t[s^{l}_{j+1}]^{A}-Mb_{z}$ (19e)
$\displaystyle t[s^{l}_{j}]^{L}\geq$ $\displaystyle
t[s^{k}_{i+1}]^{A}-M(1-b_{z})$
$\displaystyle\forall(s^{k}_{i},s^{k}_{i+1}),(s^{l}_{j},s^{l}_{j+1})$ in
conflict $\displaystyle b_{z}\in$ $\displaystyle\\{0,1\\}$ (19f)

In this formulation there are priority constraints, see (19b), (19c), and
(19d), which simply state the order of points and their minimum timing
schedule. The mutual exclusion constraints (19e) imply that no active conflict
should be present in the optimal solution.

This problem can be solved by a MILP solver or by more specialized algorithms,
as the one in [50]. As investigated in [51], instances of this problem
involving only two agents can be efficiently solved in polynomial time by an
A* search algorithm.

Note that, for the example illustrated in Fig. 3, with solution in Fig. 4, by
letting the agents move asynchronously, the cycle time can be further reduced
from 6 to 4, whereas the value obtained in the classical way can not be
shortened at all, remaining 6.8. The test case illustrated in Fig. 5 has a
cycle time of 4273 when scheduling the tours according to formulation (19).
Note that, assuming the agents being circles, the area that each agent sweeps
is illustrated and collisions between paths can be identified, see Fig. 5 and
Fig. 6. After applying the synchronous routing algorithm, the tours whose
paths collide are changed such that the potential collision areas are
decreased. This result may or may not lead to better cycle times. In this
example, by maintaining these paths, dropping the synchronous hypothesis, and
by coordinating the paths, the cycle time is improved to 3963, see Fig. 6.

Figure 5: Solution to instance 20kroA100 from TSPLIB, after B&B: swept paths
are drawn. Figure 6: Solution to instance 20kroA100 from TSPLIB, after
synchronization: swept paths are drawn.

##  VI Robot Path Planning

So far, we have assumed that the robot paths durations have already been
generated. We have neglected collisions between the robots and the
environment, and computed the time $c(q_{S},q_{E})$ it takes for one robot to
move from a start configuration $q_{S}$ to an end configuration $q_{E}$ as

$c(q_{S},q_{E})=\max\limits_{i=1,\ldots,N_{J}}\left\\{\frac{|q_{S}(i)-q_{E}(i)|}{\omega_{MAX}(i)}\right\\}$
(20)

where, $q_{S}(i)$ and $q_{E}(i)$ are the $i$-th component for the start and
end configuration, respectively, and $\omega_{MAX}(i)$ is the maximum angular
velocity for the $i$-th robot joint. Note that often industrial robots have
six rotational joints.

In reality, however, straight paths in the robot joint space may collide with
the environment. This motivates the use of robot path planning techniques to
provide the robots with point-to-point collision-free paths. Several
approaches are suitable, see [4] and [5], mostly based on sampling methods:
two of these methods are Probabilistic Roadmap Method, see [52] and [53], and
Rapidly-Exploring Random Trees (RRT), see [54]. The performance of these
algorithms heavily rely on fast and accurate collision and distance
computations, which nowadays involve triangle models and point cloud, see [55,
56, 57]. Here, we have used the built-in robot path planning functionalities
in IPS, see [2].

Since collision handling may be very computationally expensive, we adopt here
a lazy strategy, [53]. By ”lazy” it is meant that computations known to be
demanding in terms of time and space complexity are postponed, run after the
faster parts of the algorithm. In our case, see Fig. 7, the computation of
collision-free paths, the detection of robot-robot collisions and the paths
coordination are the last steps of the iterative procedure. This strategy is
based on an optimistic view of the problem: if there were no collisions at
all, then solutions obtained after the branch and bound would be optimal. The
”update data” step consists in maintaining updated costs for the paths
computed, information about which paths are colliding with each other, and,
where possible, make some logical implications about the problem in order to
derive information about the state space not yet explored. The stop criteria
usually adopted in these cases are: maximum available time or maximum number
of iterations reached, or the gap between the solution after the ”Coordinate
paths” and the one after ”Branch and Bound” under a defined threshold.

Build initial solution start Branch and Bound Synchronize Plan robot paths
update data Detect robot-robot collisions Coordinate paths is stop criterion
fulfilled?stopnoyes Figure 7: Iterative lazy optimization.

##  VII Industrial test case

These algorithms have been interfaced and tested towards the simulation
software IPS, [2], in order to handle more realistic problem instances. Here,
we show the results obtained from the study of an industrial test case: robot
models and welding points are courtesy of a major Swedish automotive
manufacturer. This industrial case consists of a stud welding station with
four robots, their 4 home positions and 32 stud points. The first solution is
obtained neglecting collisions between robots, and afterwards coordinating the
paths obtained. The second solution, instead, allows re-sequencing the tasks
within each robot, by applying the algorithm in Section V and is smoothed by
coordinating the new paths obtained, see Section V-C. The cycle time is
improved from 8.14s to 6.66s, thus a circa 18% improvement. It is important to
note that, in this case, the cycle time is not improved when considering the
cost measure in (11), but only after the synchronous hypothesis is dropped,
see Section V-C. By looking at Fig. 9, it is possible to note how the cyan
path and the blue ones are shifted in such a way that smaller collision areas
are present, w.r.t. Fig. 8.

Figure 8: Four robots with their home positions and 32 stud weld tasks:
initial solutions. Figure 9: Four robots with their home positions and 32 stud
weld tasks: solutions after being synchronized.

Another test case adapted from a stud welding station has been successfully
solved. It consists of 3 robots sharing 20 tasks. The final paths are
illustrated in Fig. 10, where their respective TCP traces are drawn. The
robots will move in a collision free way among the assigned tasks avoiding
collisions with the environment and among each others, with the aim to
minimize the makespan. These solutions were obtained in less than ten
iterations.

Figure 10: Three robots with their home positions and 20 stud weld tasks:
solutions after being synchronized. Courtesy of Volvo Cars.

##  VIII Conclusion

In this work we have studied the problem of generating, routing and scheduling
robot motions in car manufacturing stations. We have proposed an iterative
decoupled approach, considering scenarios without robot collisions, and
scenarios for highly cluttered environments. An exact solution for the problem
without collisions is obtained by a B&B algorithm that does not require LP
based tools and is very easy to implement. The lower bounds used are a
generalization of previous approaches present in the literature.

The novel algorithm to resolve conflicts between robot paths has shown good
results when solving high cluttered environments problems, giving the
possibility to modify already computed motions to avoid collisions and even to
improve cycle time.

The possibility to interact with Computer Aided Manufacturing software and,
therefore, to generate robot programs has been demonstrated on industrial test
cases.

Good results have been shown for the optimization of stations with up to 40
tasks, 4 robots and tens of possibilities to perform each task, _e.g._
geometry stations. However, efficient heuristics are needed to handle other
types of stations with more than 100 tasks, _e.g._ assembly stations.

Several questions can be further investigated. Among them, a less decoupled
B&B where at each leaf the synchronous routing and the asynchronous smoothing
are solved.

## Acknowledgement

The authors would like to thank the entire personnel at the Geometry and
Motion Planning Group at the Fraunhofer-Chalmers Research Centre for
Industrial Mathematics, for useful discussions and support in the
implementation, in particular Daniel Segerdahl. The authors also thank the
anonymous referees and associate editor for their helpful suggestions.

## References

  * [1] P. Almström, C. Andersson, A. Muhammad, M. Winroth, “Achieving Sustainable Production through Increased Utilization of Production Resource”, _Proceedings of the 4th Swedish Production Symposium, SPS11_ , Lund, May 3-5, 2011.
  * [2] www.industrialpathsolutions.com
  * [3] P. Toth, D. Vigo, _The Vehicle Routing Problem_ , Monographs on Discrete Mathematics and its applications, SIAM, 2002.
  * [4] J.-C. Latombe, _Robot Motion Planning_ , Kluwer Academic Publishing, 1992.
  * [5] S. M. LaValle, _Planning Algorithms_ , Cambridge University Press, 2006.
  * [6] D. Applegate, W. Cook, S. Dash, and A. Rohe, “Solution of a min-max vehicle routing problem”, _INFORMS Journal on Computing_ , vol. 14, pp. 132-143, 2002.
  * [7] B. Na, “Heuristic approaches for the no-depot k-traveling salesmen problem with a minmax objective”, M.S. thesis, Texas A&M University, USA, 2003.
  * [8] http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/
  * [9] J. G. Carlsson, D. Ge, A. Subramaniam, Y. Ye, “Solving the min-max multi-depot vehicle routing problem”, _Lectures on global optimization_ , vol. 55, pp. 31-46, American Mathematical Soc., 2009.
  * [10] G. Ghiani and G. Improta, “An efficient transformation of the generalized vehicle routing problem”, _European Journal of Operational Research_ , no. 122 pages 11-17, 2000.
  * [11] G. Eek and C. Eriksson, “Effective methods for solving the balanced and synchronized multiple TSP using genetic algorithms”, M.S. Thesis, University of Gothenburg, Gothenburg, Sweden, 2009.
  * [12] D. Spensieri, F. Ekstedt, J. Torstensson, R. Bohlin, J. S. Carlson, “Throughput maximization by balancing, sequencing and coordinating motions of operations in multi-robot stations”, _Proceedings of the 2010 NordDesign conference_ , Gothenburg, Sweden, August 25-27, 2010.
  * [13] J. Segeborn, D. Segerdahl, J. S. Carlson, A. Carlsson, and R. Söderberg, “Load balancing of welds in multi station sheet metal assembly lines”, _Proceedings of the ASME 2010 International Mechanical Engineering Congress & Exposition_, Vancouver, British Columbia, Canada, November 12-18, 2010.
  * [14] J. Segeborn, D. Segerdahl, F. Ekstedt, J. S. Carlson, A. Carlsson, R. Söderberg, “A generalized method for weld load balancing in multi station sheet metal assembly lines”, _Proceedings of the ASME 2011 International Mechanical Engineering Congress & Exposition_, Denver, Colorado, USA, November 11-17, 2011.
  * [15] W. Welz, “Route planning for robot systems”, Diplomarbeit Technische Universitat Berlin, January 2010.
  * [16] J. Rambau and C. Schwarz, “On the benefits of using NP-hard problems in branch & bound”, _Operations Research Proceedings 2008_ , pp. 463-468, Springer, 2009.
  * [17] L. Qiu, W.-J. Hsu, S.-Y. Huang, H. Wang, “Scheduling and routing algorithms for AGVs: a survey”, _International Journal of Production Research_ , vol. 40, no. 3, pp. 745-760, 2002.
  * [18] T. Ganesharajah, N. G Hall, C. Sriskandarajah, “Design and operational issues in AGV-served manufacturing systems”, _Annals of Operations Research_ , vol. 76, pp. 109-154, 1998.
  * [19] N. Q. Wu, W. Q. Zeng, “Deadlock avoidance in an automated guidance vehicle system using a coloured Petri net model”, _International Journal of Production Research_ , vol. 40, no. 1, 223-238, 2002.
  * [20] N. Q. Wu, M. C. Zhou, “Modeling and deadlock control of automated guided vehicle systems”, _IEEE/ASME Transactions on Mechatronics_ , vol. 9, no. 1, 50-57, 2004.
  * [21] H. Hu, Z. Li, “Modeling and scheduling for manufacturing grid workflows using timed Petri nets”, _International Journal of Advanced Manufacturing Technology_ , vol. 42, no. 5-6, pp. 553–568, May 2009.
  * [22] H. Hu, M. Zhou, Z. Li, “Supervisor design to enforce production ratio and absence of deadlock in automated manufacturing systems”, _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ , vol. 41, no. 2, pp. 201-212, March 2011.
  * [23] H. Hu, M. Zhou, Z. Li, “Liveness and ratio-enforcing supervision of automated manufacturing systems using Petri nets”, _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ , vol. 42, no. 2, pp. 392-403, February 2012.
  * [24] H. Hu, M. Zhou, Z. Li, “Supervisor optimization for deadlock resolution in automated manufacturing systems with Petri nets”, _IEEE Transactions on Automation Science and Engineering_ , vol. 8, no. 4, pp. 794-804, October 2011.
  * [25] H. Hu, Y. Liu, “Supervisor synthesis and performance improvement in automated manufacturing systems using Petri nets”, _IEEE Transactions on Industrial Informatics_ , vol. 11, no. 2, pp. 450-558, April 2015.
  * [26] N. Wu, M. Zhou, Z. Li, “Resource-Oriented Petri Net for Deadlock Avoidance in Flexible Assembly Systems”, _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ , vol. 38, no. 1, pp. 56-69, January 2008.
  * [27] Z. Li, M. Zhou, N. Wu, “A Survey and Comparison of Petri Net-Based Deadlock Prevention Policies for Flexible Manufacturing Systems”, _IEEE Transactions on Systems, Man, and Cybernetics-Part C: Applications and Reviews_ , vol. 38, no. 2, pp. 173-188, March 2008.
  * [28] Z. Li, N. Wu, M. Zhou, “Deadlock control of automated manufacturing systems based on Petri nets – a literature review”, _IEEE Transactions on Systems, Man, and Cybernetics-Part C: Applications and Reviews_ , vol. 42, no. 4, pp. 437-462, July 2012.
  * [29] T. Nishi, Y. Tanaka, “Petri Net Decomposition Approach for Dispatching and Conflict-Free Routing of Bidirectional Automated Guided Vehicle Systems”, _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ , vol. 42, no. 5, pp. 1230-1243, September 2012.
  * [30] Z. Huang, Z. Wu, “Deadlock-free scheduling method using genetic algorithms and timed $S^{3}PR$ sets”, Proceedings of the _American Control Conference_ , Boston, June 30-July 2, 2004.
  * [31] M. P. Fanti, “Event-based controller to avoid deadlock and collisions in zone-control AGVS”, _International Journal of Production Research_ , vol. 40, no. 6, pp. 1453-1478, 2002.
  * [32] M. P. Fanti, M Zhou, “Deadlock Control Methods in Automated Manufacturing Systems”, _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ , vol. 34, no. 1, pp. 5-22, January 2004.
  * [33] N. Q. Wu, M. C. Zhou, “Shortest routing of bi-directional automated guided vehicles avoiding deadlock and blocking”, _IEEE/ASME Transactions on Mechatronics_ , vol. 12, no. 1, 63-72, 2007.
  * [34] T. Nishi, M. Ando, M. Konishi, “Distributed Route Planning for Multiple Mobile Robots Using an Augmented Lagrangian Decomposition and Coordination Technique”, _IEEE Transactions on Robotics_ , vol. 21, no. 6, 1191-1200, December 2005.
  * [35] T. Nishi, R. Maeno, “Petri Net Decomposition Approach to Optimization of Route Planning Problems for AGV Systems”, _IEEE Transactions on Automation Science and Engineering_ , vol. 7, no. 3, 523-537, July 2010.
  * [36] C. W. Kim, J. M. A. Tanchoco, “Conflict-free shortest bi-directional AGV routeing”, _International Journal of Production Research_ , vol. 29, 2377-2391, 1991.
  * [37] N. N. Krishnamurthy, R. Batta, M. H. Karwan, “Developing conflict-free routes for automated guided vehicles”, _Operations Research_ , vol. 41, 1077-1090, 1993.
  * [38] S. A. Reveliotis, “Conflict resolution in AGV systems”, _IIE Transactions_ , vol. 32, Issue 7, pp. 647-659. 2000.
  * [39] T. Nishi, Y. Hiranaka, I. E. Grossmann, “A bilevel decomposition algorithm for simultaneous production scheduling and conflict-free routing for automated guided vehicles”, _Computers and Operations Research_ , vol. 38, pp. 876-888, 2011.
  * [40] A. I. Correa, A. Lagevin, L-M. Rousseau, “Scheduling and routing of automated guided vehicles: a hybrid approach”, _Computers and Operations Research_ , vol. 34, pp. 1688-1707, 2007.
  * [41] G. El Khayat, A. Lagevin, D. Riopel, “Integrated production and material handling scheduling using mathematical programming and constraint programming”, _European Journal of Operational Research_ , vol. 175, pp. 1818-1832, 2006.
  * [42] J. Clausen, “Branch and bound algorithms - Principles and examples”, Tech. Report, Dept. of Computer Science, University of Copenhagen, Denmark, March 12, 1999.
  * [43] M. Held, R. M. Karp, “A Dynamic Programming Approach to Sequencing Problems”, _Journal of the Society of Industrial and Applied Mathematics_ , vol. 10, No. 1, pp. 196-210, March 1962.
  * [44] http://www.coin-or.org/
  * [45] J. Rambau and C. Schwarz, “How to avoid collisions in scheduling industrial robots?”, Preprint, University of Bayreuth, 2010.
  * [46] N. Chakraborty, S. Akella, and J. T. Wen, “Coverage of a planar point set with multiple robots subject to geometric constraints”, _IEEE Transactions on Automation Science and Engineering_ , vol. 7, no. 1, pp. 111-122, January 2010.
  * [47] F. Ekstedt and D. Spensieri, “A direct Lin-Kernighan heuristic for the Generalized Traveling Salesman Problem”, unpublished work, 2006.
  * [48] D. Karapetyan and G. Gutin, “Lin-Kernighan heuristic adaptations for the generalized traveling salesman problem”, _European Journal of Operational Research_ , 208, pp. 221-232, 2011.
  * [49] http://www.abb.com
  * [50] D. Spensieri, R. Bohlin, J. S. Carlson, “Coordination of robot paths for cycle time minimization”, _Proceedings of the IEEE International Conference on Automation Science and Engineering_ , Madison, USA, pp- 522-527, 2013.
  * [51] A. Kobetski, D. Spensieri, M. Fabian, “Scheduling algorithms for optimal robot cell coordination - a comparison”, _Proceedings of the IEEE International Conference on Automation Science and Engineering_ , pp. 381-386, Shangai, 8-10 October, 2006.
  * [52] L. E Kavraki, P. Svetska, J.-C. Latombe, and M. Overmars, “Probabilistic Roadmaps for Path Planning in high Dimensional Configuration Spaces”, in IEEE _Transactions on Robotics and Automation_ , vol. 12, no. 4, pp. 566-580, Aug. 1996.
  * [53] R. Bohlin, L. E. Kavraki, “Path Planning using Lazy PRM”, _Proceedings of the IEEE International Conference on Robotics and Automation_ , San Francisco, CA, pp. 521-528, 2000.
  * [54] S. M. LaValle, J. J. Kuffner, “Randomized Kinodynamic Planning”, _Proceedings of the IEEE International Conference on Robotics and Automation_ , 1999.
  * [55] S. Tafuri, E. Shellshear, R. Bohlin, J. S. Carlson, “Automatic collision free path planning in hybrid triangle and point models: a case study”, _Proceedings of the 2012 Winter Simulation Conference_ , Berlin, 2012.
  * [56] E. Shellshear, R. Ytterlid, “Fast Distance Queries for Triangles, Lines, and Points using SSE Instructions”, _Journal of Computer Graphics Techniques_ , vol. 3, no. 4, 2014.
  * [57] E. Shellshear, R. Berlin, J. S. Carlson, “Maximizing Smart Factory Systems by Incrementally Updating Point Clouds”, _IEEE Computer Graphics and Applications_ , vol. 35, no. 2, pp. 62-69, 2015.

## Appendix A

###### Lemma 1

Given a best tour $T^{*}_{n}$ with length $c(T^{*}_{n})$ for an SGTSP $G_{n}$
with $n$ groups, and an SGTSP $G_{n+1}=G_{n}\cup\\{g_{n+1}\\}$, if the
triangle inequality holds, then

$c(T^{*}_{n+1})\geq\max\\{c(T^{*}_{n})+2\min(g_{n+1},G_{n})-\max(G_{n}),c(T^{*}_{n})\\}$
(21)

where $c(T^{*}_{n+1})$ is the length for an optimal tour in $G_{n+1}$,
$\min(g_{n+1},G_{n})=\min\limits_{\begin{subarray}{c}p_{i}\in G_{n}\\\
p_{k}\in g_{n+1}\end{subarray}}c(p_{i},p_{k})$, and
$\max(G_{n})=\max\limits_{\begin{subarray}{c}p_{i}\neq p_{j}\in
G_{n}\end{subarray}}c(p_{i},p_{j})$.

###### Proof:

Suppose we are given a best tour
$T^{*}_{n+1}=(p^{*}_{1},\ldots,p^{*}_{n},p^{*}_{n+1},p^{*}_{1})$, with length
$c(T^{*}_{n+1})$, in $G_{n+1}$. A feasible tour in $G_{n}$ is
$T_{n}=(p^{*}_{1},\ldots,p^{*}_{n},p^{*}_{1})$, which has a total length of

$c(T_{n})=c(T^{*}_{n+1})-c(p^{*}_{n},p^{*}_{n+1})-c(p^{*}_{n+1},p^{*}_{1})+c(p^{*}_{n},p^{*}_{1})$
(22)

A best tour $T^{*}_{n}$, thus, satisfies:

$c(T^{*}_{n})\leq
c(T^{*}_{n+1})-c(p^{*}_{n},p^{*}_{n+1})-c(p^{*}_{n+1},p^{*}_{1})+c(p^{*}_{n},p^{*}_{1})$
(23)

We have, therefore, a lower bound for $T^{*}_{n+1}$. However, this bound uses
information about a best tour in $G_{n+1}$, which we do not want to compute.
In order to eliminate such dependency, one can use the fact that
$c(p^{*}_{n},p^{*}_{n+1})\geq\min(g_{n+1},G_{n})$,
$c(p^{*}_{n+1},p^{*}_{1})\geq\min(g_{n+1},G_{n})$, and
$c(p^{*}_{n},p^{*}_{1})\leq\max(G_{n})$. These yield to

$c(T^{*}_{n+1})\geq c(T^{*}_{n})+2\min(g_{n+1},G_{n})-\max(G_{n})$ (24)

The first part of the lemma has been proved without using the triangle
inequality. For the second part, by exploiting the triangle inequality among
$p^{*}_{n}$, $p^{*}_{n+1}$, and $p^{*}_{1}$, we have, from (22):

$c(T^{*}_{n+1})-c(T_{n})=c(p^{*}_{n},p^{*}_{n+1})+c(p^{*}_{n+1},p^{*}_{1})-c(p^{*}_{n},p^{*}_{1})\geq
0$ (25)

It follows directly:

$c(T^{*}_{n+1})\geq c(T_{n})\geq c(T^{*}_{n})$ (26)

∎

###### Lemma 2

The optimum $d^{*}$ for the following problem

minimize $d$ s.t. $\displaystyle d-\sum\limits_{i:g_{i}\in G_{N}}x_{ik}$
$\displaystyle\geq d_{k}$ $\displaystyle\forall k=1,\ldots,N_{A}$ (27a)
$\displaystyle\sum\limits_{k=1}^{N_{A}}x_{ik}$ $\displaystyle=1$
$\displaystyle\forall i:g_{i}\in G_{N}$ (27b) $\displaystyle x_{ik}\in$
$\displaystyle\\{0,1\\}$ $\displaystyle\forall i:g_{i}\in G_{N}$
$\displaystyle\forall k=1,\ldots,N_{A}$

is given by

$d^{*}=\begin{cases}d_{\overline{k}}&\text{if $|G_{N}|\leq\Delta N$}\\\
\min\limits_{j=1,\ldots,N_{A}}\left\\{d_{j}+\left\lceil{\frac{|G_{N}|+\sum\limits_{k=1}^{N_{A}}\lceil
d_{k}-d_{j}\rceil}{N_{A}}}\right\rceil\right\\}&\text{otherwise}\end{cases}$
(28)

where $d_{\overline{k}}=\max\limits_{k=1,\ldots,N_{A}}\\{d_{k}\\}$ and $\Delta
N=\sum\limits_{i=1}^{N_{A}}\lfloor d_{\overline{k}}-d_{i}\rfloor$.

## Appendix B

Some problems remain still unsolvable, even if a feasible solution exists, see
Fig. 11. A feasible solution is, _e.g._ , $(H_{1},H_{2})$, $(A,H_{2})$,
$(A,B)$, $(A,H_{2})$, $(H_{1},H_{2})$, $(H_{1},C)$, $(H_{1},H_{2})$, whereas
by the transformation to a GTSP it is not possible to find any.

$H_{1}$$A$$B$$H_{2}$$C$ Figure 11: An example with two agents where there
exists a feasible solution, but by the transformation to a GTSP it is not
possible to find any.


# Complete Description of Measures Corresponding to Abelian Varieties over
Finite Fields

Nikolai S. Nadirashvili and Michael A. Tsfasman  CNRS, France  Institute for
Information Transmission Problems, Moscow, Russia<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Date: 2021)

###### Abstract.

We study probability measures corresponding to families of abelian varieties
over a finite field. These measures play an important role in the
Tsfasman–Vlăduţ theory of asymptotic zeta-functions defining completely the
limit zeta-function of the family. J.-P. Serre, using results of R.M. Robinson
on conjugate algebraic integers, described the possible set of measures than
can correspond to families of abelian varieties over a finite field. The
problem whether all such measures actually occur was left open. Moreover,
Serre supposed that not all such measures correspond to abelian varieties (for
example, the Lebesgue measure on a segment). Here we settle Serre’s problem
proving that Serre conditions are sufficient, and thus describe completely the
set of measures corresponding to abelian varieties.

###### Key words and phrases:

Abelian varieties over finite fields, Weil numbers, asymptotic zeta-function

###### 1991 Mathematics Subject Classification:

11G10, 11G20

Supported in part by ANR project FLAIR (ANR-17-CE40-0012)

##  1\. Preliminaries

###  1.1. Abelian varieties over a finite field

Let us recall the notation and background from [Se] and [Ts 1]. We fix once
and for all a finite field $\mathbb{F}_{q}$. Given an algebraic variety $A$ of
dimension $g$ over it, consider the multiset $\Omega=\\{\omega_{j}\\}$ of its
inverse Frobenius roots. These $2g$ roots are

  * —

pairwise complex-conjugate (real roots having even multiplicity);

  * —

of absolute value $\sqrt{q}$ ;

  * —

roots of a monic polynomial with integer coefficients (i.e., they are
algebraic integers and if $\omega\in\Omega$ then all its
$\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-conjugates are also in
$\Omega$).

Call the multiset $\\{\omega_{j}\\}$ with these properties an _integer Weil
system_ or just a _Weil system_. The Weil system determines the zeta-function
of $A$.

Instead of considering Weil systems $\Omega=\\{\omega_{i}\\}$ on the circle we
can sum pairs of complex conjugate ones and consider
$X=\\{x_{i}\\},i=1,...,g$, where $x_{i}=\omega_{i}+\bar{\omega}_{i}$.

To the Weil system $\Omega$ there corresponds the probability measure

$\mu_{\Omega}={\frac{1}{g}}\sum_{x_{i}\in X}\delta_{x_{i}},$

where $\delta_{x_{i}}$ is the Dirac skyscraper measure at $x_{i}$. We work in
the weak topology space of measures dual to continuous functions on the
segment $I=[-2\sqrt{q},2\sqrt{q}]$, writing the pairing as $\langle
f,\mu\rangle=\int fd\mu$.

We are interested in the corresponding asymptotic problem.

For a family of Weil systems $F=\\{\Omega_{j}\\}$ of growing cardinality
$g_{i}$ we consider the limit in the weak topology

$\mu_{F}=\lim_{j}\mu_{\Omega_{i}}.$

Of course, it need not exist. When it exists, we say that the family is
_asymptotically exact_. Such is for example any tower of embedded abelian
varieties. Each family contains an asymptotically exact subfamily, since the
space of positive measures of given mass on the circle being compact. This
notion is compatible with that of asymptotic exactness for curves and abelian
varieties used in [Ts 2, Ts/Vl 1, Ts/Vl 2].

We face the following natural questions, cf. [Ts 1]:

  1. Q1.

_Describe the set of all limit measures corresponding to asymptotically exact
families of abelian varieties over $\mathbb{F}_{q}$._

  2. Q2.

_Describe the set of all limit measures corresponding to asymptotically exact
families of curves over $\mathbb{F}_{q}$._

The second question is treated in [Ts 2], [Ts/Vl 1], [Ts 1], and — together
with its analog for number fields — in [Ts/Vl 2]. The answer a priori looks
rather strange. In particular, all limit measures have continuous density. But
this is specific for curves. Moreover, for curves — not even to speak about
number fields — the complete answer is still unknown (see [Ts 1]).

For abelian varieties the situation is quite different. For example, the limit
measure can be atomic.

###  1.2. Algebraic integers

By the Honda–Tate theorem [Ta] for any $\Omega$ there exist an abelian variety
over $\mathbb{F}_{q}$ whose set of inverse Frobenius roots is a multiple of
$\Omega$. (The word multiple here is important, $\Omega$ itself is not always
realisable.)

The set $\Omega=\\{\omega_{j}\\}$ of cardinality $2g$ uniquely determines and
is uniquely determined by the set
$\bar{\Omega}=\\{\omega_{j}+\bar{\omega}_{j}\\}$ of cardinality $g$, where for
each pair of conjugate elements of $\Omega$ we take their sum. $\bar{\Omega}$
is a set of $g$ totally real algebraic numbers on the segment
$I=[-2\sqrt{q},2\sqrt{q}]$. Thus the description of possible sets $\Omega$
becomes the question in algebraic number theory:

  * —

describe algebraic intergers such that all their
$\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-conjugates are in
$[-2\sqrt{q},2\sqrt{q}]$.

For each $\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-stable multiset
$X=\\{x_{i}\\},i=1,...,g$ of totally real algebraic numbers lying on the given
segment $I$ we define the probability measure (i.e., measure of total mass 1),
normalizing the sum of atomic measures at these numbers,
$\mu_{X}={\frac{1}{g}}\sum_{i=1}^{g}\delta_{x_{i}}$ . Consider a sequence
${\bf X}=\\{X_{j}\\}$ of such sets with growing $g_{j}$ and the limit measure
$\mu_{\bf X}=\lim{\frac{1}{g_{j}}}\mu_{X_{j}}$. Of course, such measure need
not exist. If the limit exists, we say that the family $\bf X$ is
_asymptotically exact_. (For us a _family_ means just a _sequence with
growing_ $g_{j}$.)

Any family contains an asymptotically exact subfamily. Any tower (family with
$X_{j}\supset X_{j-1}$) is asymptotically exact.

We come to the following problem in algebraic number theory, equivalent to the
above question Q1.

  1. Q3.

Describe the set of all limit probability measures corresponding to
asymptotically exact families of
$\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-stable multisets of totally
real algebraic integers lying on the given segment $I$.

If we take some $\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-stable
multiset $X$ and set $X_{n}=nX,n\to\infty$, then the limit is just $\mu_{X}$
which is an atomic measure. (This corresponds to the family of powers of a
given abelian variety.)

Any atomic $\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-stable measure
supported in totally real algebraic integers on $I$ can be obtained in this
way.

Recall that every positive measure $\mu$ on a segment is a sum of $\mu_{a}$
and $\mu_{d}$, where $\mu_{a}$ is atomic and $\mu_{d}$ is atom-free (or
diffuse), i.e., $\mu_{d}(x)=0$ for any point $x$.

So our question is reduced to the following

  1. Q4.

Describe the set of all limit _diffuse_ probability measures corresponding to
asymptotically exact families of
$\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-stable multisets of totally
real algebraic integers lying on the given segment $I$.

###  1.3. Capacity and equilibrium measure

Following [Se] we introduce some notation. Consider a compact
$K\subset\mathbb{C}$. For a positive measure $\mu$ on $K$ define

$I(\mu)=\int\int_{K\times K}\ln|x-y|\mu(x)\mu(y).$

Note that $I(\mu)\in\mathbb{R}\cup\\{-\infty\\}$. Set $v(K)=\sup_{\mu}I(\mu)$,
where $\mu$ runs over all positive probability measures supported in $K$. Then
the _capacity_ of $K$ is defined as ${\rm{cap}}(K)=e^{v(K)}$. Its logarithm
$v(K)$ is called the _logarithmic capacity_ of $K$.

For a non-compact $Y\subset\mathbb{C}$ the capacity is defined as supremum of
capacities of its compact subsets.

For a compact $K$ with non-zero capacity there exists a unique positive
probability measure $\mu$, such that $I(\mu)=\ln\rm{cap}(K)$. This measure is
diffuse. Its support can be smaller than $K$, but for $K\subset\mathbb{R}$
their difference has capacity $0$. This measure $\mu=\mu_{K}$ is called the
_equivilibrium measure_ of $K$.

###  1.4. Algebraic integers with all conjugates in a given compact

For a set $V\subset\mathbb{C}$ let ${\rm{Irr}}_{V}$ be the set of irreducible
monic polynomials in $\mathbb{Z}[x]$ of degree at least $1$ such that all
their roots lie in $V$. For such a polynomial $P(x)$ of degree $g$ let
$\mu_{P}$ be the corresponding probability measure supported in its roots,
$\mu_{P}={\frac{1}{g}}\sum_{i=1}^{g}\delta_{x_{i}}$. Now let
$K\subset\mathbb{C}$ be compact. There are two quite different cases (Fekete,
Szegë, Robinson, cf. [Se]), depending on the capacity of $K$.

  * —

If $\rm{cap}(K)<1$, i.e., $v(K)<0$, then $\rm{Irr}_{K}$ is finite.

  * —

If $K\subset\mathbb{R}$ is a union of finite number of segments and
$\rm{cap}(K)>1$, i.e., $v(K)>0$, then $\rm{Irr}_{K}$ is infinite.

  * —

If $K\subset\mathbb{C}$ is
$\operatorname{Gal}(\bar{\mathbb{Q}}/\mathbb{Q})$-stable and $\rm{cap}(K)\geq
1$, i.e., $v(K)\geq 0$, then for any open $U\supset K$ the set $\rm{Irr}_{U}$
is infinite.

###  1.5. Space of measures corresponding to integer polynomials

Consider probability measures $\mu_{P}$, corresponding to irreducible
polinomials with all roots in $K$. The set of such measures is enumerable;
introduce any order on it and let ${\bf M}_{n}$ be the closure in the weak
topology of the convex envelope of all such measures except for the first
$n-1$ ones. Then

${\bf M}_{1}\supset{\bf M}_{2}\supset{\bf M}_{2}\supset\dots$

Let ${\bf M}_{\infty}=\bigcap_{n}{\bf M}_{n}$. We know [Se] :

  * —

The set ${\bf M}_{\infty}$ is non-empty, convex, and compact.

  * —

Any measure $\mu\in{\bf M}_{\infty}$ is diffuse.

  * —

The capacity $\rm{cap}(\rm{Supp}(\mu))$ of its support is greater than or
equal to $1$.

  * —

If $\rm{cap}(\rm{Supp}(\mu))=1$ then $\mu=\mu_{K}$ is the equilibrium measure
of $K$.

###  1.6. Potential

Potential $p_{\mu}:\mathbb{C}\to\mathbb{R}\cup\\{-\infty\\}$ of the measure
$\mu$ is defined as

$p_{\mu}(z)=\int_{\mathbb{C}}\ln|w-z|\mu(w).$

Potential $p_{K}$ of the equilibrium measure $\mu_{K}$ of $K$ is everywhere in
$\mathbb{C}$ greater than or equal to the logarithmic capacity $v(K)$.
Moreover, it equals $v(K)$ on $K$ and is harmonic outside $K$.

##  2\. Criterion

###  2.1. Main result

We are going to prove the following

######  Theorem 2.1.

Let $I$ be a segment $[a,b]\subset\mathbb{R}$ and let $K$ run over all compact
subsets of $I$. For a positive probability meausre $\mu$ whose support lies in
$I$ the following conditions are equivalent:

  1. (1)

$\mu$ lies in the weak topology closure of equilibrium measures $\mu_{K}$ for
$K$ of capacity $\rm{cap}(K)\geq 1$, i.e., with $v(K)\geq 0$;

  2. (2)

$\mu$ lies in the weak topology closure of equilibrium measures $\mu_{K}$ for
$K$ of capacity $\rm{cap}(K)=1$, i.e., with $v(K)=0$;

  3. (3)

$p_{\mu}(z)\geq 0$ for all $z\in\mathbb{C}$.

###  2.2. Application to ${\bf{M}}_{\infty}$

Recall the theorem of Serre ([Se], thm 1.6.3), following from the above
results of Robinson (see 1.4).

######  Theorem 2.2.

Let $K\subset\mathbb{R}$ be a finite union of segments with $\rm{cap}(K)>1$,
$E\subset K$ be its compact subset with $\rm{cap}(E)\geq 1$, and $\mu_{E}$ be
its equilibrium measure. Then $\mu_{E}\in{\bf{M}}_{\infty}$.

From this theorem and Theorem 2.1 we immediately derive the following
criterion.

######  Theorem 2.3.

Let $K\subset\mathbb{R}$ be a finite union of segments, then a positive
probability measure $\mu$ on $K$ lies in ${\bf{M}}_{\infty}$ if and only if
its potential is non-negative, i.e., $p_{\mu}(z)\geq 0$ for all
$z\in\mathbb{C}$.

######  Corollary 2.4.

A positive difuse probability measure on the segment
$I_{q}=[-2\sqrt{q},2\sqrt{q}]$ is the limit measure for an asymptotically
exact family of abelian varieties if and only if $p_{\mu}(z)\geq 0$ for all
$z\in\mathbb{C}$.

######  Remark 2.5.

In particular, we get an affirmative answer to the question (1.6.8) of Serre
[Se]: the normalized Lebesgue measure of mass $1$ of $I_{q}$ lies in ${\bf
M}_{\infty}$ for any $q\geq 2$. Indeed, for the Lebesgue measure supported on
the segment $[-a,a]$ the potential is

$p_{L}(y)=(y+a)(\ln{(y+a)}-1)+(a-y)(\ln{(a-y)}-1).$

It reaches its minimum at $y=0$ and $p_{L}(0)=2a(\ln a-1)$, which is positive
for $a\geq e$, and $2\sqrt{2}>e$.

##  3\. Generalities on potential

We set ${\mathbb{C}}={\mathbb{R}}^{2}$, $z=x+iy$. By
$D_{r}\subset{\mathbb{R}}^{2},\,r>0,$ we denote the disk $|z|<r$, and by
$D^{x}_{r}$ the disk $|z-x|<r$; by $S_{a}\subset{\mathbb{R}}^{2},\,a>0$, the
slab $-a<y<a$; by $R(a,b,c)\subset{\mathbb{R}}^{2},\,c>0,$ the rectangle
$a<x<b,\,-c<y<c$; by $I\subset{\mathbb{R}}^{2}$ the segment $[(-Q,0),(Q,0)]$.

Let $E{\subset\mathbb{R}}^{2}$ be compact. Denote $u_{E}(z)=p_{E}(z)-v(E)$,
where $p_{E}(z)$ is the logarithmic potential of the equilibrium measure of
the set $E$. We call function $u_{E}(z)$ the normalised potential of the set
$E$. Define harmonic functions $u_{E}^{R}(z),\,R>Q$ in $D_{R}\setminus E$
satisfying the following boundary conditions: $u_{E}^{R}(z)=0$ on the regular
points of $\partial E$ and $u_{E}^{R}(z)=1$ on $\partial D_{R}$. Then $(\ln
R)u_{E}^{R}(z)$ tends to $p_{E}(z)=u_{E}(z)+v(E)$ as $R\to\infty$.

The basic properties of the logarithmic capacity and of the solution of the
exterior Dirichlet problem are well known see, e.g. [Do], cf. (Part 1, and
especially XIII.18)

(1) Let $E\subset I$ and $(-a,a)\subset E$. Then $u_{E}^{R}(z)<C(\ln|z|-\ln
a)$ for $|z|>2a$, where $C>0$ is an absolute constant.

Let $\mu$ be a probability measure on $I$, and

$p_{\mu}(z)=\int_{\mathbb{C}}\ln|w-z|\mu(w)$

its potential. The following properties of $p_{\mu}$ are easy to check:

(2) $p_{\mu}(z)$ is subharmonic. It is harmonic on
${\mathbb{R}}^{2}\setminus\rm{Supp}(\mu)$.

(3) $p_{\mu}(z)$ is an even function of $y$.

(4) $p_{\mu}(z)$ is a monotonically increasing function of $y$ for $y>0$ and a
monotonically decreasing function of $y$ for $y<0$.

(5) $p_{\mu}(z)=\ln|z|+o(1)$ as $|z|\to\infty$, and $\mu=2\pi\Delta
p_{\mu}(z)$.

(6) Assume that $p_{\mu}(z)\geq 0$. Since $p_{\mu}(z)>1$ on
${\mathbb{R}}^{2}\setminus D_{r}$, if $r>Q+e$, then (by the Harnack inequality
for harmonic functions) for any $a>0$ there exists $c(a)>0$ such that
$p_{\mu}(z)>c(a)$ in ${\mathbb{R}}^{2}\setminus S_{a}$.

(7) Let $\\{\mu_{n}\\}$ be a sequence of probability measures on $I$. Then,
$\mu_{n}\to\mu$ in the weak topology if and only if $p_{\mu_{n}}(z)\to
p_{\mu}(z)$ pointwise.

(8) Let $E\subset{\mathbb{R}}^{2}$ be a compact nonempty set with a smooth
boundary. Let $f$ be a continous function on $\partial E$. Denote $M^{-}=\inf
f,\,M^{+}=\sup f$. Choose a constant $m$ such that $M^{-}\leq m\leq M^{+}$.
For any large enough $R$ we can solve the Dirichlet problem:

$\Delta u_{R}=0\;{\rm in}\;D_{R}\setminus E,$ $u_{R}=f\;{\rm on}\;\partial
E,\;\;u_{R}=m\;{\rm on}\;\partial D_{R}.$

Choose a convergent sequence $u_{r_{n}},\;r_{n}\to\infty$ in
${\mathbb{R}}^{2}\setminus E$. Then the function

$u=\lim u_{r_{n}}$

is a solution of the external Dirichlet problem

$\Delta u=0\;in\;{\mathbb{R}}^{2}\setminus E,$ $u=f\;{\rm on}\;\partial E.$

By the maximum principle each of the functions $u_{r_{n}}$ satisfies the
inequalities

$M^{-}\leq u_{r_{n}}\leq M^{+},$

and hence the function $u$ satisfies the same inequalities

$M^{-}\leq u\leq M^{+}.$

We show now that the last Dirichlet problem has a unique bounded solution.
Assume by contradiction that for a given function $f$ there are two different
bounded solutions $u_{1}$ and $u_{2}$. Then $u=u_{1}-u_{2}$ is a bounded
harmonic function vanishing on $\partial E$.

Assume now that $D_{r_{0}}\subset E\subset D_{r_{1}}$. Let $v_{R}$ be a
solution of the Dirichlet problem

$\Delta v_{R}=0\;{\rm in}\;D_{R}\setminus E,$ $v_{R}=0\;{\rm on}\;\partial
E,\;\;v_{R}=\ln(R/r_{0})\;{\rm on}\;\partial D_{R}$

Then $\ln|x|/r_{1}<v(x)<\ln|x|/r_{0}$ and we can choose a sequence
$r_{n}\to\infty$ such that $v_{r_{n}}\to u$ in ${\mathbb{R}}^{2}\setminus E$.
Then $v$ is a positive harmonic function in ${\mathbb{R}}^{2}\setminus E$,
vanishing on $\partial E$ and tending to infinity at infinity. Thus for any
positive $\varepsilon$ for sufficiently large $R$ in $D_{R}\setminus E$ by the
maximum principle we have $u<\varepsilon v$. Since the last inquality holds
for any positive $\varepsilon$ we get $u=0$.

(9) As an immediate consequence of (8) we have the following maximum principle
for harmonic functions in the unbounded domains in dimension $2$. Let
$E\subset{\mathbb{R}}^{2}$ be compact, and let $u$ be a bounded harmonic
function in ${\mathbb{R}}^{2}\setminus E$ continuous in its closure. Then the
following maximum principle for the exterior domain holds:

$\inf_{\partial E}u<u<\sup_{\partial E}u.$

(10) Let $s=[-\varepsilon,\varepsilon]$ be a segment on the axis $x$, and $u$
be a harmonic function in $D_{1}\setminus s$ such that $u=0$ on $s$ and $u=1$
on $\partial D_{1}$. Then $u<C\ln K/\ln(1/\varepsilon)$ in $D_{K\varepsilon}$
for any $1<K<1/\varepsilon$, here $C>0$ is an absolute constant.

(11) Let $x_{1},\dots,x_{n}\in[0,1],\;l_{1},\dots,l_{n}<\varepsilon$. Let $X$
be the union of segments

$X=\cup_{i=1}^{n}[x_{i}-l_{i},x_{i}+l_{i}]$

and

$U_{K}=\cup_{i=1}^{n}D_{Kl_{i}}^{x_{i}},$

where $K>1$. Assume $v(X)<0$. Since $p_{X}(z)<\ln(1+R)$ on $D_{R},R>2$ then
there is an absolute constant $C>0$ such that $p_{X}(z)<v(X)-Cv(X)\ln K$ on
$U_{K}$. Now, from (10) we get

$v(U_{K})<v(X)-Cv(X)\ln K.$

##  4\. Proof of Theorem 2.1

Any measure $\mu$ can be appoximated by a sequence $\mu_{n}$ of atomic
measures. Our goal is to prove that if for all $z$ the potential
$p_{\mu}(z)\geq 0$, then for every large enough $n$ we can approximate these
$\mu_{n}$ by $\mu^{\prime}_{n}$, where $\mu^{\prime}_{n}$ is the equilibrium
measure of a set $S_{n}$ which is a union of small segments containing
${\rm{Supp}}(\mu_{n})$.

(12) Let $\mu_{n}$ be a sequence of probability measures, each supported in a
finite set of points of $I$, $I=[a,b]$, such that $\mu_{n}$ tends to $\mu$
weakly. Then $p_{\mu_{n}}(z)\to p_{\mu}(z)$ in ${\mathbb{R}}^{2}\setminus
S_{\delta}$ uniformly for any $\delta>0$.

(13) Assume that $p_{\mu}(z)\geq 0$. Denote
$G_{n}:=\\{z\in{\mathbb{R}}^{2},p_{\mu_{n}}(z)<0\\}$.

From (12) it follows that for any $\varepsilon>0$ there exists some
$N=N(\varepsilon)$ such that $G_{n}\subset S_{\varepsilon}$ for any $n>N$.

(14) Let $i_{a},0\leq a\leq 0$ be the segment $[-a,a],y=0$, $v(z)$ be a
harmonic in $D_{1}\setminus i_{a}$ function, $v(z)=0$ on $i_{a}$ and $v(z)=1$
on $\partial D_{1}$. Then $v(z)<C\ln(|z|/a)/\ln(1/|a|)$ in $D_{1}\setminus
D_{2a}$ where $C>0$ is an absolute constant.

(15) The set $I_{n}=G_{n}\cap I$ is a finite union of segments, which we
denote $[(a_{i}^{n},0),(b_{i}^{n},0)]$. Let $p_{I_{n}}(z)$ be the potential of
the equilibrium measure of the set $I_{n}$. In the following sections we prove
that $u_{I_{n}}(z)$ tends pointwise to $p_{\mu}(z)$ as $n$ goes to infinity.
Obviously it is sufficient to prove that

$|p_{\mu_{n}}(z)-u_{I_{n}}(z)|\to 0\;.$

Let $H_{n}\subset D_{R}$ be a sequence of open sets such that $I_{n}\subset
H_{n}$. Note that by (9) the result will follow if we prove that there is a
sequence $\varepsilon_{n}>0,\,\varepsilon_{n}\to 0$ such that

$|p_{\mu_{n}}(z)-u_{I_{n}}(z)|<\varepsilon_{n}\;\;{\rm on}\;\;\partial
H_{n}\;.$

Note that since $u_{I_{n}}(z)\geq 0$, we have $u_{I_{n}}(z)-p_{\mu_{n}}(z)\geq
0$ on $I_{n}$ and hence everywhere. If this is proved, then we have
$u_{I_{n}}(z)-p_{\mu_{n}}(z)\to 0$ on an open set, hence from the last
inequality, from (9), and from the Harnack inequality
$u_{I_{n}}(z)-p_{\mu_{n}}(z)\to 0$ pointwise everywhere.

We still have to choose the sets $H_{n}$ and to prove the above inequality on
them.

(16) Fix $k$. Let $G_{k}\subset S_{\varepsilon}$. Then

$G_{k}\subset\cup_{i}R(a_{i}^{k},b_{i}^{k},\varepsilon)\;.$

Set $l_{i}=b_{i}^{k}-a_{i}^{k}$.

For $l_{i}<\varepsilon$ we want to prove the inequality

$p_{\mu_{k}}(z)>-C_{1}e^{-C_{2}A}\;\;{\rm on}\;\;\partial
R(a_{i}^{k},b_{i}^{k},Al_{i}),$

where $C_{1},C_{2}$ are absolute constants, $A>1$.

Indeed, denote $G^{i}=G_{k}\cap R(a_{i}^{k},b_{i}^{k},\varepsilon)$. Then
$G^{i}$ is a compact subdomain of the rectangle
$R(a_{i}^{k},b_{i}^{k},\varepsilon)$ and $p_{\mu_{k}}(z)=0$ on $\partial
G^{i}$. Since the equilibrium measure $\mu_{k}$ of the segment
$[(a_{i}^{k},0),(b_{i}^{k},0)]$ is less than or equal to $1$, it follows that

$p_{\mu_{k}}(z)>-C_{0}\;\;{\rm on}\;\;G^{i}\setminus
R(a_{i}^{k},b_{i}^{k},l_{i}),$

where $C_{0}>0$ is an absolute constant. By the maximum principle for harmonic
functions
$p_{\mu_{k}}(z)>-2C\cos(2\pi(x-(b_{i}^{k}+a_{i}^{k})/2)/l_{i})e^{-y/l_{i}}$,
where $z=(x,y)$ and we get the desirable inequality.

(17) We have proved a lower estimate for the function $p_{\mu_{k}}$ on the set
$\partial G_{k}$. In the following sections we shall get an upper estimate for
the function $u_{I_{n}}$ on this set and, as a consequence, the convergence
demanded in (15).

(18) First we prove

$\lim\inf v(I_{n})>-\infty$

as $n\to\infty$.

Assume the opposite, i.e.,

$\lim\inf v(I_{n})=-\infty.$

Denote $l_{i}^{k}=b_{i}^{k}-a_{i}^{k}$,
$G^{\prime}_{n}:=\\{z\in{\mathbb{R}}^{2},p_{\mu_{n}}(z)<-C_{0}\\}$. From (16)
it follows that

$G^{\prime}_{k}\subset\cup_{i}R(a_{i}^{k},b_{i}^{k},l_{i}^{k}).$

.

From our assumption and (11) it follows

$\lim v(\cup_{i}R(a_{i}^{k},b_{i}^{k},l_{i}^{k}))=-\infty$

and hence

$\lim v(G^{\prime}_{n})=-\infty.$

Since $v(G_{n})=0$ and $v(G^{\prime}_{n})=v(G_{n})+C_{0}=C_{0}$, the result
follows.

(19) Let $U(z)$ be the normalized potential of the equilibrium measure of the
segment $e=[-1,1]$ on the real axis. Then

$U(z)<C+\ln|z|$

near infinity and

$U(z)<C\sqrt{{\rm dist}(z,e)}$

near the segment $e$. The last inequality is well known and it follows for
example from the explicit expression for the equilibrium measure of the
segment $[a,b]$ (cf. [Se], A2.4)

$\mu_{[a,b]}=\frac{1}{\pi}\frac{dx}{\sqrt{(b-x)(x-a)}}\;.$

We will use the last inequalities for potential $U$ to get upper bounds for
the potentials $u_{I_{k}}(z)$.

Fix $k$ and let $l_{i}=b_{i}^{k}-a_{i}^{k}$. Denote by $c_{i}$ the centre of
the segment $[a_{i},b_{i}]$, $c_{i}=((a_{i}+b_{i})/2,0)$. Denote by $H$ an
annulus centered at $c_{i}$ with inner radius $l_{i}$ and exterior radius
$2R$. Then as it follows from (18), $u_{I_{k}}(z)<C$ on the exterior boundary.
Thus

$u_{I_{k}}(z)<CU((z-c_{i})/l_{i})/\ln(1/l_{i})\;.$

The above inequalities for the function $U$ imply the following estimates. Fix
$\varepsilon>0$.

(a) For $l_{i}>\varepsilon$ the normalized potential
$u_{I_{k}}(z)<C_{\varepsilon}\sqrt{\varepsilon/l_{i}}/\ln{(1/l_{i})}$ on
$R(a_{i}^{k},b_{i}^{k},\varepsilon)$, where $C>0$ is an absolute constant.

(b) If $\varepsilon/A<l_{i}\leq\varepsilon,\,A>1$ then
$u_{I_{k}}(z)<C_{A}\ln(A)/\ln({1/l_{i}})$ on
$R(a_{i}^{k},b_{i}^{k},\varepsilon)$, where $C_{A}>0$ is a constant depending
on $A$.

(20) In this section we define the sets $H_{k}$ to verify the inequality of
(15). By (13) there is a sequence $c_{k}\to 0$ such that $G_{k}\subset
S_{c_{k}}$. For any $\varepsilon>0$ later on we shall define a sufficiently
large constant $C>0$ and a sufficiently small constant $\delta>0$. Denote
$l_{i}^{k}=b_{i}^{k}-a_{i}^{k}$. For $l_{i}\leq\delta$ set

$R_{i}^{k}=R(a_{i}^{k},b_{i}^{k},Cl_{i}^{k}).$

For $l_{i}>\delta$ set

$R_{i}^{k}=R(a_{i}^{k},b_{i}^{k},d_{i}^{k}),$

where $d_{i}^{k}=\min\\{l_{i}^{k}/C,c_{k}\\}$. Now let

$H_{k}=\cup_{i}R_{i}^{k}\;.$

Then for a sufficiently large constant $C>0$ and a sufficiently small constant
$\delta>0$ by (16) and (19) we get the inequality

$|p_{\mu_{n}}(z)-u_{I_{n}}(z)|<\varepsilon\;\;{\rm on}\;\;\partial H_{k}\;.$

(21) We show that $\Delta p_{\mu_{n}}(z)-\Delta u_{I_{n}}(z)\to 0$ as
$n\to\infty$ in the sense of distributions. Since $\Delta p_{\mu_{n}}(z)$,
$\Delta u_{I_{n}}(z)$ are probability measures and since the space of
probability measures on a segment is weakly compact, we can choose a
subsequence such that the following limits exist : $\Delta
p_{\mu_{n}}(z)\to\mu_{1},\;\Delta u_{I_{n}}(z)\to\mu_{2}$ weakly as
$n\to\infty$, where $\mu_{1},\mu_{2}$ are probability measures on $I$. Since
$p_{\mu_{n}}(z)$ is a potential of a probability measure and since weak
convergence of measures implies convergence of their potentials, we may assume
convergence $p_{\mu_{n}}(z)\to h(z),\,z\in R^{2}\setminus I$, where $h$ is a
harmonic function in the complement of $I$. Hence by (20) $u_{I_{n}}(z)\to
h(z)$. Thus $\mu_{1}=\mu_{2}$. Therefore $\Delta p_{\mu_{n}}(z)-\Delta
u_{I_{n}}(z)\to 0$ as $n\to\infty$ in the sense of distributions.

Denote by $\nu_{n}$ the equilibrium measure of the set $I_{n}$. Then
$\nu_{n}=\Delta u_{I_{n}}$. Thus $\nu_{n}-\Delta p_{\mu_{n}}(z)\to 0$ weakly.
Therefore, it follows that $p_{I_{n}}(z)\to p_{\mu}(z)$ weakly, and since
$p_{I_{n}}(z)=u_{I_{n}}(z)+v(I_{n})$ we get $v(I_{n})\to 0$ and
$|p_{\mu_{n}}(z)-p_{I_{n}}(z)|\to 0$. We can choose a sequence $\delta_{n}\to
0$ such that if $I^{\prime}_{n}$ is the image of $I_{n}$ after delation of the
axis $x$ to $(1+\delta_{n})x$ then ${v}(I^{\prime}_{n})=0$. Since
$\delta_{n}\to 0$, we have $p_{I^{\prime}_{n}}(z)\to p_{\mu}(z)$ weakly.

The theorem is proved.

## References

  * [Do] J.L. Doob Classical Potential Theory and Its Probabilistic Counterpart, Springer, 1983
  * [Se] J.-P. Serre, “Distribution asymptotique des valeurs propres des endomorphismes de Frobenius”, Sém. Bourbaki 70:1146 (mars 2018).
  * [Ta] J. Tate, “Classes d’isogénie des variétés abéliennes sur un corps fini (d’après Honda),” Sém. Bourbaki, 358, 1968-69.
  * [Ts 1] M. A. Tsfasman, “Serre’s Theorem and Measures Corresponding to Abelian Varieties over Finite Fields”, Moscow Math. J., 2019, v.19, n.4, pp. 789-806
  * [Ts 2] M. A. Tsfasman, “Some remarks on the asymptotic number of points,” Springer Lect. Notes in Math., v. 1518, pp. 178–192, 1992.
  * [Ts/Vl 1] M. A. Tsfasman and S. G. Vlǎduţ, “Asymptotic properties of zeta-functions.”, J. Math. Sci. (New York), 1997, v.84:5, pp. 1445-1467
  * [Ts/Vl 2] M. A. Tsfasman and S. G. Vlǎduţ, “Infinite global fields and the generalized Brauer–Siegel theorem”, Moscow Math. J., 2002, v.2, n.2, pp. 329-402
  * [We] A. Weil, Variétés Abéliennes et Courbes Algébriques, Hermann, Paris, 1948\.