# Meta-operators for Enabling Parallel Planning
Using Deep Reinforcement Learning

Ángel Aso-Mollar1, Eva Onaindia1

###### Abstract

There is a growing interest in the application of Reinforcement Learning (RL)
techniques to AI planning with the aim to come up with general policies.
Typically, the mapping of the transition model of AI planning to the state
transition system of a Markov Decision Process is established by assuming a
one-to-one correspondence of the respective action spaces. In this paper, we
introduce the concept of meta-operator as the result of simultaneously
applying multiple planning operators, and we show that including meta-
operators in the RL action space enables new planning perspectives to be
addressed using RL, such as parallel planning. Our research aims to analyze
the performance and complexity of including meta-operators in the RL process,
concretely in domains where satisfactory outcomes have not been previously
achieved using usual generalized planning models. The main objective of this
article is thus to pave the way towards a redefinition of the RL action space
in a manner that is more closely aligned with the planning perspective.

## Introduction

Integrating AI planning and Machine Learning algorithms is a hot research
topic that has reported significant advances in many different scenarios, such
as training general policies with supervised or unsupervised Learning,
(Ståhlberg, Bonet, and Geffner 2022a), (Ståhlberg, Bonet, and Geffner 2022b),
Generalized Planning using Reinforcement Learning (RL), (Gehring et al. 2022),
(Rivlin, Hazan, and Karpas 2020), (Francès et al. 2019), learning Action
Models using Reinforcement Learning (Ng and Petrick 2019), and many others.

Specifically, strands of this discipline have integrated AI Planning with RL
usually to train, with synthesized experience, heuristics that guide the
search process of a planning agent, through rectification of domain-
independent heuristics (Gehring et al. 2022) or through direct policy training
(Rivlin, Hazan, and Karpas 2020), for example. The training process is usually
based on giving the agent a reward of 1 when it reaches the goal of the
planning problem and 0 otherwise. This introduces the so-called sparse reward
problem wherein we observe that the agent rarely receives rewards from the
environment, thus hindering its ability to learn. It is also known that some
works attempt to overcome these problems using, for instance, Reward Machines
(Icarte et al. 2022).

In addition, it is common for papers that use RL with AI planning to make use
of the very action definition in planning by applying it directly into
Reinforcement Learning: to establish the mapping between the transition model
of AI planning and the state transition system of a Markov Decision Process,
it is assumed that there is a one-to-one correspondence between their
respective action spaces. This is what we refer to as the 1-1 mapping between
planning and RL actions. In this paper, we propose a new perspective based on
the concept of meta-operator, with which we decouple this 1-1 mapping often
seen in the literature, therefore allowing to apply several actions at the
same time (in parallel) in a concrete time step. Parallel plans increase the
efficiency of planning because increased parallelism leads to a reduction of
the plan length or the number of time points. In this paper, we also
investigate if general parallel policies can be obtained with a generalized
planning framework and if such policies improve general sequential policies.

By defining the concept of meta-operator as the application of several
planning actions at the same time, we can force the RL training process to
simulate parallelism. In addition, we can also use meta-operators to further
guide the training process, alleviating the sparse reward problem by
establishing a small reward when the meta-operator contains more than one
atomic action.

We incorporate meta-operators in the context of Generalized Planning and to
this aim our approach is based on a modified version of the architecture
proposed in (Rivlin, Hazan, and Karpas 2020), where Graph Neural Networks are
used to generate a compact representation of the planning states. We will
report results showing that the inclusion of meta-actions allows almost a 100%
coverage in domains where it is difficult to reach generalization with
traditional learning methods such as logistics or depots, as stated in
(Ståhlberg, Bonet, and Geffner 2022a).

We tested our models in problem instances from the International Planning
Competitions111https://www.icaps-conference.org/competitions/ (IPC) and in
randomly generated ones, and we observed that the coverage using meta-
operators improves with respect to not using them. It is also observed that
the length of the plan, as expected, is also reduced when we include meta-
operators.

This paper is structured as follows: several basic concepts with which we rely
our work on are going to be defined, such as Planning, RL and their
integration, including Generalized Planning. Next, the core element of this
article will be defined: the meta-operator, and then will be integrated into
the RL structure. Finally, several different experiments that justify the use
of meta-operators in RL will be conducted, followed by a discussion of whether
to use this mechanism and the benefits against other models that do not
include this approach.

## Background

### Planning

Classical planning is the problem of finding a sequence of actions that when
applied to an initial state lead to the achievement of a goal or set of goals.
The two main components of a planning task are the domain and the problem.

A planning domain is defined as $D=\langle\mathcal{F},\mathcal{A}\rangle$,
where $\mathcal{F}$ is a set of fluents that describe lifted properties of the
domain objects and their relations, and $\mathcal{A}$ is a set of action
schemas in which every $a\in\mathcal{A}$ is defined by a triplet
$\langle\textup{Pre}(a),\textup{Add}(a),\textup{Del}(a)\rangle$,
$\textup{Pre}(a),\textup{Add}(a),\textup{Del}(a)\subseteq\mathcal{F}$.

A planning problem instance linked to a domain
$D=\langle\mathcal{F},\mathcal{A}\rangle$ is defined as a tuple $P=\langle
F,O,I,G\rangle$ where $F$ is the set of facts that result from grounding the
fluents $\mathcal{F}$ to the objects of the problem $P$, and $O$, called
operators, is the result of grounding the action schemas $\mathcal{A}$. A
state $s\subseteq F$ is a set of facts, and $I,G\subseteq F$ are,
respectively, the initial and the goal state of the problem.

Every operator $o\in O$ has its corresponding grounded triplet defined by
$\langle Pre(o),Add(o),Del(o)\rangle$. $Pre(o)\subseteq F$ are the
preconditions that must be true in a state for the operator to be applicable;
that is, an operator $o$ is applicable in a state $s\subseteq F$ if
$Pre(o)\subseteq s$. $Add(o)\subseteq F$ are the effects of the operator that
assert a positive literal after the application of $o$; and $Del(o)\subseteq
F$ is the set of negative effects, i.e., the set of literals which become
false after the operator is applied. Applying $o\in O$ in $s$ leads to a new
state $s^{\prime}=(s\setminus Del(o))\cup Add(o))$.

A solution to a planning problem $P=\langle F,O,I,G\rangle$ is a sequence of
operators or plan $\rho=\langle o_{1},o_{2},\ldots,o_{k}\rangle,\ o_{i}\in O\
\forall i$, such that the result of applying the operator sequence $\langle
o_{1},o_{2},\ldots,o_{k}\rangle$ to the initial state $I$ leads to a state
$s\subseteq F$ that satisfies $G\subseteq s$, that is to say, that the state
$s$ contains every goal in $G$.

### Reinforcement Learning

Reinforcement Learning (RL) is a computational approach to learning from
environmental interaction (Sutton and Barto 1998). The main objective of RL is
to learn a policy, or behaviour, that maximizes the positive interaction that
an agent gets through time.

RL scenarios are often modeled as finite Markov Decision Processes (MDP)
(Puterman 1990). An MDP is a control process that stochastically models
decision-making scenarios; an agent continuously interacts with an environment
by executing actions that change its internal state, and these actions are
rewarded consequently. The main objective of the agent is to improve its
performance so that it can progressively maximize the received reward.

Formally, an MDP is defined as $M=\langle S,A,R,Pr\rangle$, where $S$ is a set
of states, $A$ is a set of actions, $R:S\times A\times S\rightarrow\mathbb{R}$
is a reward function that values how good or bad it is to take an action
$a_{t}$ at a certain state $s_{t}$ that leads to another state $s_{t+1}$,
$R(s_{t},a_{t},s_{t+1})=r_{t}$, and $Pr$ is a transition probability function
such as $Pr(s^{\prime},s,a)=p(s^{\prime}|s,a)$, which is usually unknown.

At each time step $t$, an agent takes an action in the current state $s_{t}$
among a set of applicable actions $A_{t}\subseteq A$ in $s_{t}$ following a
probability distribution called policy, represented by $\pi:S\rightarrow
p(\cdot|\ S=s_{t})$. The objective of RL is to reward decisions that reach the
target as quickly as possible, by prioritizing immediate rewards.

The RL optimization problem is formally defined as

$\max_{\pi}\mathbb{E}_{\pi}G_{t}$ (1)

where $G_{t}=\sum_{k=0}^{\infty}\gamma^{k}r_{t+k+1}$ is the sum of the
discounted reward of a trajectory sampled from the policy $\pi$ from the
current state $t$ to the infinite. Thus, formula (1) measures the reward
received by all possible trajectories multiplied by the probability of taking
the decisions of the trajectories under the policy $\pi$, where
$\gamma\in[0,1]$ is the discount factor of future rewards.

### Integrating RL on top of planning

The integration of Planning and RL is a promising approach for solving control
problems. Planning is effective to reason over a long horizon but assumes
access to a local schema; in contrast, RL is effective at learning policies
and relative values of states but fails to plan over long horizons (Eysenbach,
Salakhutdinov, and Levine 2019). Combining planning and RL can thus be
regarded as learning a control policy for the planner in lieu of using a
planning heuristic function.

Combining Planning and RL enables addressing a variety of tasks like learning
generalist policies to solve a large number of problems (Celorrio, Aguas, and
Jonsson 2019), learning action models based on interaction (Moerland et al.
2023), training intelligent heuristics that guide the agent in the control
task (Gehring et al. 2022), or even training the agent to solve concrete
planning problems. Generally, the RL module is integrated inside the planning
agent to help it make a decision, although other approaches use it the
opposite way (Lee et al. 2022). In this work, we will follow the first
approach.

As the formulation of both RL and Planning use the concept of state and
action, we must establish the equivalence between RL actions and planning
operators, so that the planning and RL modules can communicate with each
other.

It is usually assumed in the literature that for a planning problem $P=\langle
F,O,I,G\rangle$ the equivalence is established by simply considering the
planning operators $O$ as the RL actions $A$ ($A=O$). States in RL ($S$) are
abstract concepts that simply represent situations that occur in the
environment, so we can inherit the concept of planning state and apply it
directly to RL by making every $s\in S$ as $s\subseteq F$; i.e., we can assume
RL states as planning states. Hence, the notion of an action being applicable
to a state $s$ in RL follows directly from the above considerations and uses
the same principle of applicability of an operator in the planning state $s$.

Since this one-to-one correspondence between actions and operators yields huge
action spaces, some researchers adopt certain considerations to reduce the
number of possible actions. For example, in PDDLGym (Silver and Chitnis 2020),
authors consider only a partial grounding of the action schema in order to
reduce the number of possible actions. Particularly, they claim that certain
operators contain parameters that correspond to free objects (in terms of
controlling the agent), that is, objects whose properties can potentially
change, leaving ungrounded the non-free agents. For example, in a transport
domain in which we want to transport packages between cities using trucks, the
packages and trucks would be free objects because their location changes, but
the cities would not. In this scenario, an action schema like (move ?t - truck
?c1 - city ?c2 - city) will only bound the parameter truck to a concrete
object as the two parameters city refer to static objects that only denote the
structure of the movement network.

We can think of this consideration as a transformation $g$ of the planning
operator space, $g(O)$, which is then used as the RL action space, $A=g(O)$.
Note that in PDDLGym this specification is hard-coded, while we propose an
automated method, as we will discuss later.

We already mentioned that the notion of a planning state is directly
translated into an RL state. We need to take into account that sometimes, and
depending on the RL method, the planning states must be translated into some
encoding that the RL method supports. For example, if we are using a NN to
represent the policy, we need a vectorial representation of the state. This
representation can be obtained, for example, by using Graph Neural Networks
(Zhou et al. 2020) or Neural Logic Machines (Dong et al. 2019), in conjunction
with intermediate structures such as graphs.

Sometimes we will say that we apply a learned policy $\pi$ to a planning
problem $P$, meaning that the planning operator corresponding to the RL action
given by the policy is successively applied from the initial state of the
problem to the goal, and states will change accordingly through the planning
operators application. This way, as mentioned above, the planning states are
mapped to the RL states and (hegemonically) planning operators are mapped to
RL actions.

### Generalized Planning

Generalized Planning is a discipline in which we aim to find general policies
that explain a set of problems, usually of different size. Specifically, given
a set of problems $\\{\langle
F_{i},O_{i},I_{i},G_{i}\rangle\\}_{i=1}^{N},N>0$, rather than searching for a
solution plan $\rho_{i}=\langle o_{1},o_{2},\ldots,o_{k}\rangle,\ o_{j}\in
O_{i}$ for each individual problem $P_{i}$ by applying a policy $\pi_{i}$
previously trained for $P_{i}$, we aim to find a general policy $\pi$ that
when applied to every problem of the set returns a solution plan to all of
them.

## Defining meta-operators

As we discussed in the previous section, a bijective correspondence is usually
considered between the set of planning operators and the set of RL actions
when solving a planning problem $P=\langle F,O,I,G\rangle$ of a domain $D$.
That is, $A=O$. We have also seen that we can transform this operator space
into a different one $g(O)$. In this work, unlike PDDLGym that performs a
reduction of the operator space, we will define a transformation function to
enrich it.

For this purpose, we define the notion of a meta-operator for a problem
$P=\langle F,O,I,G\rangle$ as a new operator that, given $L>1$ and
$o_{1},...,o_{L}\in O$, which we will refer to as atomic operators, is
represented as

$\bigoplus_{i=1}^{L}o_{i}$

with the following considerations:

  * •

$Pre\left(\bigoplus_{i=1}^{L}o_{i}\right)=\bigcup_{i=1}^{L}Pre(o_{i})$

  * •

$Add\left(\bigoplus_{i=1}^{L}o_{i}\right)=\bigcup_{i=1}^{L}Add(o_{i})$

  * •

$Del\left(\bigoplus_{i=1}^{L}o_{i}\right)=\bigcup_{i=1}^{L}Del(o_{i})$

  * •

No pair of actions $o_{i}$, $o_{j}$ from all the atomic actions
$o_{1},...,o_{L}$ that form the meta-operator conflict with each other. Two
atomic operators $o_{i}$ and $o_{j}$ conflict with each other if one of these
two conditions hold:

    * –

$\exists\ p\in Pre(o_{i})$ such that $p\in Del(o_{j})$.

    * –

$\exists\ p\in Add(o_{i})$ such that $p\in Del(o_{j})$.

Broadly speaking, a meta-operator is nothing but a synthetic operator
resulting from the union of atomic ones that can be executed in any order. The
resulting operator will therefore inherit the union of its sets $Add$, $Pre$
and $Del$, always bearing in mind that all the atomic operators involved can
be executed at the same time, that is, that they do not conflict with each
other.

This means that they are inconsistent if they compromise the consistency of
the resulting state when applying the operator, i.e., if the resulting state
changes if the sequential order of application of the atomic operators
changes. These two considerations are equivalent to the notions of
inconsistent effects and interference in the calculation of a mutex realtion
between two actions in GraphPlan (Blum and Furst 1997).

We define the set of meta-operators of degree $L\in\mathbb{N}$ in a problem
$P$ of a domain $D$ as the union of every possible meta-operation:

$\mathcal{O}^{L}=\bigcup_{o_{i}\in O}\left[\bigoplus_{i=1}^{L}o_{i}\right]$

where $o_{i}$ are actions from $O$ that do not interfere with each other when
defining a single meta-operator and $\mathcal{O}^{1}=O$.

## Including meta-operators in RL

Meta-operators are then added to the RL action space, enriching and enabling
the application of parallel planning operators within this sequential RL
action space. We can therefore define a new transforming function $g_{L}$ as
the union of the base set with meta-operators of degree $L$:

$g_{L}(O)=\bigcup_{i=1}^{L}\mathcal{O}^{i}$

and train our RL algorithms using this consideration, i.e., $A=g_{L}(O)$.

This integration of meta-operators is calculated online, as in Algorithm 1, at
each time step $t$ and current state $s_{t}$, for all applicable RL actions
$A_{t}=\\{o\in O:o\textup{ is applicable at }s_{t}\\}$. We define that a meta-
operator is applicable at state $s_{t}$ if every atomic operator is also
applicable at $s_{t}$ and if operators do not interfere with each other, which
is followed by definition.

Algorithm 1 Calculate all applicable RL actions

Input:
A set of applicable planning operators $O$
Degree of meta-operators $L$
Output:
A set of applicable RL actions $A$

1: $A=O$

2: $N=\emptyset$

3: # Generate all conflicts

4: for $a\in O$ do

5: for $p\in Pre(a)$ do

6: for $b\in O$ do

7: if $p\in Del(b)$ then

8: $N=N\cup\\{(a,b),(b,a)\\}$ # Interference

9: for $p\in Add(a)$ do

10: for $b\in O$ do

11: if $p\in Del(b)$ then

12: $N=N\cup\\{(a,b),(b,a)\\}$ # Inconsistent effects

13: # Combine pairs, triplets, until $L$, of operators

14: # that do not contain any pair of operators in $N$.

15: for $i\in\\{2,...,L\\}$ do

16: $A=A\cup\textup{MakeMetaOperators}(O,L,N)$

17: return A

We used the Generalized Planning RL training scheme proposed in (Rivlin,
Hazan, and Karpas 2020) to observe the effects of the inclusion of a meta-
operator in domains that often struggle to generalize in the literature, such
as logistics or depots (Ståhlberg, Bonet, and Geffner 2022a). That
architecture also uses GNNs for state representation and gives a reward of 1
to the agent if it reaches the goal, and 0 in any other case, as usual.

This last decision is highly criticized by the planning community because it
introduces the so-called sparse rewards problem, that is, the agent receives
information from the environment at very specific moments, thus hindering the
learning process. The inclusion of meta-operators opens up the possibility of
defining a certain reward for actions that include more than one atomic
operator, briefly alleviating the aforementioned problem.

In that sense, we want to test how this inclusion of different specific
amounts of reward in meta-operators affects to the learning process, so we
will do an analysis using different amounts of reward and we will see how well
they train and how parallel generated plans are compared with each other.

We also want to test whether the inclusion of meta-operators actually improves
the coverage for problems of domains we analyzed, compared to the coverage of
a sequential model, trained with the same domain but without using meta-
operators ($A=O$).

## Experiments

In this section, we present a series of experiments that support the inclusion
of meta-operators in Generalized Planning using RL. In particular, we are
interested in two things: (1) analyzing the impact of an extra reward when a
meta-operator is applied in the learning process, and (2) checking whether the
inclusion of meta-operators improves the results in terms of coverage (number
of solved problems).

Specifically, we conduct two experiments. Experiment 1 is designed to measure
the degree of parallelism of the solution plans using different rewarding in
meta-operators. Experiment 2 evaluates the performance of our model against
two different defined datasets.

### Domains

We will use two domains that are widely used in the IPC and also known for
their complexity, logistics and depots, and a third domain which is an
extension of the well-known blocksworld domain.

#### Multi-blocksworld.

This domain is an extension of the blocksworld domain that features a set of
blocks on an infinite table arranged in towers, with the objective of getting
a different block configuration by moving the blocks with robot arms. Blocks
can be put on top of another block or on the table, and they can be grabbed
from the table or from another block. We have defined two robot arms.

#### Logistics.

This domain features packages located at certain points which must be
transported to other locations by land or air. Ground transportation uses
trucks and can only happen between two locations that are within the same
city, while air transportation is between airports, which are special
locations of the cities. The destination of a package is either a location
within the same city or in a different city. In general, ground transportation
is required to take a package to the city’s airport (if the package is not at
the airport). The package is then carried by air between cities, and finally
using ground transportation the package is delivered to the final destination
if its destination is not the arrival airport.

#### Depots.

This domain consists of trucks that are used for transporting crates between
distributors, and hoists to handle the crates in pallets. Hoists are only
available at certain locations and are static. Crates can be stacked/unstacked
onto a fixed set of pallets. Hoists do not store crates in any particular
order. This domain slightly resembles the multi-blocksword domain as there is
a stacking operation, though crates do not need to be piled up in a specific
order, and to the logistics domain as to the existence of agents that
transport crates from one point to another.

### Data generation

RL algorithms need a large number of instances in order to converge. That is
why, for the training process, it was necessary to use automatic generators of
planning problems. For logistics and depots domains, we used generators of the
AI Planning Github (Seipp, Torralba, and Hoffmann 2022), while for the multi-
blocksworld domain we created a new generator based on the generator for the
blocksworld domain (Seipp, Torralba, and Hoffmann 2022).

Domain | Train size | Total objects train | Validation/Test size | Total objects test
---|---|---|---|---
Multi-blocksworld | 5-6 blocks | 5-6 | 10-11 blocks | 10-100
Logistics |  | 2-4 airplanes
---
2-4 cities
2-4 trucks
2-4 locations per city
1-3 packages
9-10 |  | 3-4 airplanes
---
6-7 cities
3-4 trucks
6-7 locations per city
6-7 packages
24-29
Depots |  | 1-2 depots
---
2-3 distributors
2-3 trucks
3-5 pallets
2-4 hoists
3-5 crates
13-22 |  | 5-6 depots
---
5-6 distributors
5-6 trucks
5-6 pallets
5-6 hoists
5-6 crates
30-36
Table 1: Sizes used for the problem generation, in terms of general and
specific objects.

Table 1 illustrates the size distribution of the problems used in this work;
it shows the number of objects of each type involved in all three domains. We
generated a dataset of random problems out of the distributions shown in Table
1, which we will refer to as Dataset 1 (D1).

#### Dataset 1 (D1)

It consists of problems uniformly sampled from the test distribution of Table
1 and generated with the aforementioned generators. We generated 460 problems
for the multi-blocksworld domain, 792 problems for the logistics domain and
640 for the depots domain, as a result of creating ten problems for each
configuration in the test distribution.

Additionally, we created a second collection of samples that we will refer to
as Dataset 2 (D2) from a renowed planning competition.

#### Dataset 2 (D2)

It consists of problems that were used in the IPC (specifically, in the
IPC-2000 and IPC-2002). We used the 35 first instances for the blocksworld
domain in IPC-2000, with a slight modification to introduce the two robot
arms; the 30 first instances of logistics from IPC-2000; and the 22 instances
of depots from IPC-2002. This set of instances was chosen in order to compare
our results with those obtained in the work (Ståhlberg, Bonet, and Geffner
2022b).

### Setup

We opted for using a $L=2$ degree meta-operator to come up with a feasible
extension of the action space. As the inclusion of meta-operators increases
the action space, we need to find a balance between size and performance.
Using two-degree meta-operators is sufficient to fulfill the two objectives
mentioned at the beginning of this section, namely analyzing the impact of
rewarding meta-operators and evaluating the coverage of the models. We will
also evaluate how much does the action space rise with our approach compared
to a sequential model.

The RL training was conducted on a machine with a Nvidia GeForce RTX 3090 GPU,
a 12th Gen Intel(R) Core(TM) i9-12900KF CPU and Ubuntu 22.04 LTS operating
system, and the same hyperparameter configuration than (Rivlin, Hazan, and
Karpas 2020). A similar training process as the one proposed in (Rivlin,
Hazan, and Karpas 2020) was followed here, all policies are trained for 900
iterations, each with 100 episodes and up to 20 gradient update steps, using
Proximal Policy Optimization RL training algorithm, with a discount factor of
0.99.

### Experiment 1: Rewarding of meta-operators

In this experiment, we aim to observe how the reward granted to the
application of a meta-operator in the RL training influences the learning
process and the quality of the plans. We are interested in measuring the
effect of meta-operators in terms of the plan length or the number of time
steps of the plan. To that end, we define the parallelism rate of a solution
plan of a problem as:

$\frac{\\#\ parallel\ operators}{\\#\ total\ plan\ timesteps}$

where $\\#\ parallel\ operators$ is the number of meta-operators that appear
in the plan by applying the learned policy to the problem, and $\\#\ total\
plan\ timesteps$ is the total number of timesteps of the plan. This is a
measure of how frequently parallel operators appear in the decisions made by
the planner agent.

We trained a series of models giving different reward values to meta-
operators. This experiment can be thought of as a way of tuning the meta-
operators reward, which can therefore be regarded as a hyperparameter. Since
we primarily aim to find the most appropriate reward for the use of meta-
operators in this experiment, we decided to focus only on the training
distribution.

We trained five models from the train distribution of Table 1 with reward
values to meta-operators of $0.0,0.1,0.01,0.001$ and $0.0001$, respectively.
Subsequently, the five models were run on a fixed sample, generating 10
problems for each element from the train distribution, and results were
analyzed in order to obtain the average parallelism rate of all plans.

During the experiment execution, rewards and the number of parallel actions at
each time step are monitored so as to balance out the reward coming from
parallel actions and coming from achieving a solution plan. In other words, we
want to avoid situations in which parallel actions are just added for the sake
of reward, which may deviate plans towards a large number of parallel actions
sacrificing reaching the objective.

Reward | Multi-blocks | Logistics | Depots
---|---|---|---
0.0 | 0.550 | 0.243 | 0.530
0.1 | - | - | -
0.01 | 0.701 | 0.851 | 0.857
0.001 | 0.559 | 0.582 | 0.381
0.0001 | 0.557 | 0.768 | 0.294
Table 2: Average parallelism rate for all models trained with specified reward
for the application of meta-operators.

The results of this experiment are shown in Table 2: all models are able to
fulfill the aforementioned objective (100% of coverage in training) except for
the model that gives a reward of 0.1 to meta-operators (no results are shown
because the model did not converge). Intuition tells us that there are certain
values that reward too much parallelism, even above reaching the problem goal
itself, resulting in potentially infinite plans that execute parallel actions
in a loop (until the maximum episode time is reached). This means that a
reward value of 0.1 for meta-operators outputs policies that yield more than
ten parallel actions per plan, which exceeds the value given to reaching the
goal, thus making the algorithm converge to a situation in which no goal is
reached but lots of meta-operators of degree $L>2$ appear in the plan.
Ultimately, RL is about optimizing a reward function, and if adding meta-
operators produces more reward, this will be the path taken by the model.

According to Table 2, we observe that the model that gives the best results in
terms of parallelism rate for all domains within the plans generated is the
0.01 reward model. This indicates that, in order to obtain potentially better
results, a balance must be established between the amount of reward given to
parallelism and the amount of reward given to the goal.

For example, a somewhat more conservative proposal, which we know for sure
would not exceed the goal reward, is to establish a meta-operator reward of
$\frac{GOAL\\_REWARD}{MAX\\_ITERATIONS}$, where $GOAL\\_REWARD$ is the reward
given to reaching the goal (generally 1) and $MAX\\_ITERATIONS$ is the maximum
number of applications of the policy before stating that the goal cannot be
achieved. Generally, this approach is excessively limiting and does not
encourage parallelism. This is evident from Table 2, which shows that greater
rewards lead to improved parallelism.

In fact, the amount of reward provided in meta-operators is also dependent on
the average length of the plans we want to test. That means, perhaps if the
problems we want to test have a larger average plan length than the ones we
trained, it would be wiser to test with a model that has been trained with a
slightly lower reward in order to not “overflow” the reward and fall into the
undesired scenario of policies that produce parallel actions with no goal
termination. This would be a problem that would probably occur in Generalized
Planning, for example, where we train models with a problem size smaller than
the problem size on which we will test the results.

All in all, it has been found that the amount of reward given to meta-
operators is significant in terms of quality and convergence of plans.

### Experiment 2: Performance in Generalized Planning

In this experiment, we compare the original sequential model with versions of
the parallel model obtained with different reward values. We note that the aim
is to test the performance of the models when dealing with new inputs of a
larger size. The trained policy for each domain is then analyzed as in the
Generalized Planning literature by testing the problems in datasets D1 and D2.
Particularly, for each model, we measure the coverage and the average length
of the generated plans for the problems in D1 and D2.

Domain - D1 | OR | R = 0.0 | R = 0.001 | R = 0.01
---|---|---|---|---
Multi-blocks (460) | 268 | 439 | 408 | 406
Logistics (792) | 131 | 701 | 717 | 317
Depots (640) | 287 | 572 | 640 | 552
Multi-blocks | 79.99 | 76.43 | 65.92 | 92.84
Logistics | 200.50 | 112.61 | 109.05 | 110.25
Depots | 338.49 | 106.68 | 124.46 | 121.80
Domain - D2 | OR | R = 0.0 | R = 0.001 | R = 0.01
Multi-blocks (35) | 2 | 35 | 34 | 34
Logistics (30) | 11 | 26 | 28 | 27
Depots (22) | 20 | 17 | 20 | 20
Multi-blocks | 172.00 | 43.17 | 44.59 | 54.12
Logistics | 136.64 | 115.92 | 120.64 | 121.00
Depots | 127.45 | 83.59 | 110.20 | 114.25

Table 3: Results for datasets D1 and D2. The top part of the tables shows
coverage out of the total number of instances shown between parenthesis. The
bottom part of the tables indicates the average plan length. OR is the
original sequential model in (Rivlin, Hazan, and Karpas 2020); R=0.0 is the
model with no meta-operator reward; R=0.1 is the model with reward of $0.01$,
and R=0.001 is the model with reward of $0.001$.

Table Experiment 2: Performance in Generalized Planning shows the results
obtained with the sequential (OR) model (Rivlin, Hazan, and Karpas 2020), the
parallel model trained with no reward (R=0.0), with a reward of 0.01 on meta-
operators (R=0.01) and with a reward of 0.001 on meta-operators (R=0.001) for
the International Planning Competition (D2) and randomly generated (D1)
datasets. With these experiments we aim to illustrate how the results vary
from one model to another depending on the reward, as stated in the previous
section.

The table is divided in two halves, one for each set of problems. In the top
part of each half we show results for coverage of the analyzed models with
respect to the problems of each set, while in the bottom part of each half the
average length of plans for each problem of the set is analyzed. In the top
part of each half we show within parentheses the total number of problems that
compound that set, and then in each column the number of those for which the
model under analysis has managed to reach the target. The bottom part of each
half corresponds to the number of time steps with which the models managed to
solve each set of problems. We present the average number of actions taking
into account only the solved problems.

Results of Table Experiment 2: Performance in Generalized Planning show that
coverage from the models that use meta-operators improves with respect to the
coverage of the sequential model. Moreover, it is observed that for the
logistics domain, which has been found in the literature to be hard to address
for RL approaches (Ståhlberg, Bonet, and Geffner 2022b), a very good result is
obtained with respect to that of the sequential model.

Results are correlated to the enrichment of the action space. Table 4 shows
the size of the action space visited during training of the sequential model
and the parallel ones. There is only one result in Parallel because the action
space is the same for all three models, only the amount of reward given to the
meta-operator application varies. We observe that in the multi-blocks domain
where there are only two agents (two robot-arms), the increase in the number
of actions is much less significant than in logistics or depots domain, in
which many more agents (trucks, planes, etc.) coexist at the same time.

Plans also tend to be shorter in terms of average actions when applying meta-
operators. In Table Experiment 2: Performance in Generalized Planning, results
for R=0.0 tend to give short plan lengths, which makes sense because as a
sparse reward problem plans tend to reach the goal faster thanks to the
discount factor and reward only coming from reaching the goal. This sometimes
even result in better coverage, such as in multi-blocks, while in R=0.001 and
R=0.01 the reward also comes from the meta-operators and not only from
reaching the goal.

| Multi-blocks | Logistics | Depots
---|---|---|---
Sequential | 100 | 108 | 228
Parallel | 1140 | 3960 | 8519
Table 4: Action space or number of RL actions (planning operators and, when
makes sense, meta-operators) visited during training.

The use of a model with lower reward R=0.001 has reported significantly better
coverage results than the model with R=0.01. In the Generalized Planning task
there is a variance in the size of problems tested, resulting also in a
variance in the length of its correspondent plans. As the model R=0.01 gives a
high reward to parallelism, if the size of the plans is too large, parallelism
is being rewarded too much. For example, 11 meta-operators would already mask
the objective’s reward, which is 1, i.e., $11\cdot 0.01=1.1>1$.

## Discussion

Although there is not yet a standard neural network architecture that suits
perfectly the planning frame, the approaches that leverage deep learning
techniques that automatically extract structure from high-level data are
mostly based on graph convolution to learn state embeddings like Graph Neural
Networks (Rivlin, Hazan, and Karpas 2020; Ståhlberg, Bonet, and Geffner 2022b)
or first-order logic neural models like Neural Logic Machines (Gehring et al.
2022).

Overall, the aforementioned approaches show to be competitive with baseline
sequential models that use planning heuristics or state-of-the-art
implementations of classical planners in terms of coverage (number of solved
problems). More interestingly, some works report the inability of NN-based
heuristics to outperform classical planning engines in transport-like domains
wherein a tight coupling between the different objects of a planning task
exists (Rivlin, Hazan, and Karpas 2020; Ståhlberg, Bonet, and Geffner 2022b;
Gehring et al. 2022).

Frameworks that combine AI planning and RL establish a one-to-one
correspondence between operators in planning and actions in RL. It is
observed, as stated before, that low performance usually arises in domains
with a high number of agents that need to collaborate with each other to reach
the goal, also called tightly-coupled domains (Torreño, Onaindia, and Sapena
2014), for example logistics or depots.

Specifically, speaking about logistics, in (Ståhlberg, Bonet, and Geffner
2022b) authors prove that it is not possible to achieve satisfactory results
using a GNN architecture. The authors identify transitive relations without
which the locality of the objects is lost. For example, when a package is
inside a truck and the truck is in a city, the information about the city
where the package is found fails to propagate in a transitive manner, i.e., it
fails to infer that the package is indeed in the city. Results within this
perspective are improved by hard-coding extra predicates (derived atoms) and
including them to guarantee the transitive relations. In the mentioned
example, a new atom is added to denote that a package is in a city.

Interestingly, we found that in tightly-coupled domains such as logistics,
depots or multi-blocksworld test results from our perspective are very
satisfactory when including meta-operators. In fact, if we compare the same
IPC instances of the logistics domain from (Ståhlberg, Bonet, and Geffner
2022b) with our approach we observe that, while they have a coverage of 17
(28) without using the hard-coded atoms, we obtain 26 (28). It should be
emphasized that our process is automatic, i.e., it is not necessary to
manually enrich any problem or domain to reach the goal.

It is interesting to think about why adding meta-operators actually improves
the training process when it comes to using RL for planning. It is not only
the fact that it allows several atomic operators to coexist at the same time,
decentralizing as already mentioned in other sections the purely sequential
nature of RL, but also the fact that when adding the extra operators the
algorithm converges faster and yields better results.

We hypothesize that the inclusion of meta-operators opens a way of (virtual)
communication between agents. As stated in (Rivlin, Hazan, and Karpas 2020),
problems in tightly-coupled domains arise when certain resources need to be
used by several agents at the same time. Since a policy is nothing more than
the “brain” of the planning agent, when several entities can be executed in
parallel (e.g. moving different vehicles) or change their state (location,
mode, etc.) at the same time, the sequential nature of the RL may interfere
with this notion of collaboration, thus obstructing the convergence of RL
methods for this type of problems.

By introducing parallel actions we allow several entities to change their
state at the same time, so that they can “prepare” for what is next to come. A
clear example is in the logistics domain: if we have a package being
transported in a truck, which must then be picked up by another entity that
transports it elsewhere, in a sequential nature the truck would leave the
package at one location, and then the other entity would have to pick it up
and take it away. However, how can one distinguish, given that the only known
information is that the package is at one location, that the package is to be
picked up by the next entity and not by the same entity? By introducing
parallelism, the second entity could go to the exchange point before the
package actually reaches that exchange point. We believe that our approach
provides some sense of “entity” to the independent objects of a planning
problem and improves agent collaboration.

On the other hand, we also wanted to hypothesize about why meta-operators also
favor exploration in the state space of planning. If we think about state
graphs as in (Bonet and Geffner 2022), which are graphs that outline the whole
structure of actions and states that the agent can take at each moment, we can
consider adding meta-operators to the RL action space as introducing extra
edges to these state graphs.

In the end, training with RL is done on a trial-and-error basis, choosing
actions randomly at first and seeing if it brings some benefit. By including
meta-operators we allow our planning agent to do in one sampling step what
previously had to be done in more (as long as the actions can coexist in
parallel, obviously). This helps the convergence and thus the generalization
of our methods.

To illustrate this last affirmation, we want to give another example. Let’s
assume an initial training iteration with a policy that is initialized with a
uniform probability $p$ to each RL action. If the agent decides to apply a
meta-operator of degree $L$, we arrive to a state $s$ with a probability of
$p$, but in a sequential environment we would have to arrive to $s$ (following
the same path) with $L$ different transitions, each one corresponding to the
atomic operators that conform the meta-operator, with a $p^{L}$ probability.
This summarizes why we think it is interesting to apply meta-operators.

## Conclusion and Future work

In conclusion, we defined the notion of meta-operator as a way to enrich the
RL structure that has hegemonically been used to attack the planning problem.
Furthermore, we saw that it improves results from other sequentially trained
models, both with randomly generated problems and with problems from the
International Planning Competition. Finally, we wanted to highlight the
discussion that may arise from this paper, with special emphasis on the fact
that RL structures can be enriched to better match the true nature of
planning.

Finally, we are aware that the inclusion of meta-operators makes the action
space larger, so a balance between parallelism and size must also be found so
that training times are not excessively long. For future work we would like to
explore different ways of attacking RL problems with much larger, even
continuous, action spaces in order to computationally be able to use meta-
operators of any degree. We would also like to explore other different ways of
representing planning states and how this may influence what was analyzed and
to try to explain how does tuning the meta-operators reward affect in several
different domains.

## Acknowledgments

This work is supported by the Spanish AEI PID2021-127647NB-C22 project and
Ángel Aso-Mollar is partially supported by the FPU21/04273.

## References

  * Blum and Furst (1997) Blum, A.; and Furst, M. L. 1997.  Fast Planning Through Planning Graph Analysis.  _Artif. Intell._ , 90(1-2): 281–300.
  * Bonet and Geffner (2022) Bonet, B.; and Geffner, H. 2022.  Language-Based Causal Representation Learning.
  * Celorrio, Aguas, and Jonsson (2019) Celorrio, S. J.; Aguas, J. S.; and Jonsson, A. 2019.  A review of generalized planning.  _Knowl. Eng. Rev._ , 34: e5.
  * Dong et al. (2019) Dong, H.; Mao, J.; Lin, T.; Wang, C.; Li, L.; and Zhou, D. 2019.  Neural Logic Machines.  _CoRR_ , abs/1904.11694.
  * Eysenbach, Salakhutdinov, and Levine (2019) Eysenbach, B.; Salakhutdinov, R.; and Levine, S. 2019.  Search on the Replay Buffer: Bridging Planning and Reinforcement Learning.  In _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , 15220–15231.
  * Francès et al. (2019) Francès, G.; Corrêa, A. B.; Geissmann, C.; and Pommerening, F. 2019\.  Generalized Potential Heuristics for Classical Planning.  In Kraus, S., ed., _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019_ , 5554–5561. ijcai.org.
  * Gehring et al. (2022) Gehring, C.; Asai, M.; Chitnis, R.; Silver, T.; Kaelbling, L. P.; Sohrabi, S.; and Katz, M. 2022.  Reinforcement Learning for Classical Planning: Viewing Heuristics as Dense Reward Generators.  In Kumar, A.; Thiébaux, S.; Varakantham, P.; and Yeoh, W., eds., _Proceedings of the Thirty-Second International Conference on Automated Planning and Scheduling, ICAPS 2022, Singapore (virtual), June 13-24, 2022_ , 588–596. AAAI Press.
  * Icarte et al. (2022) Icarte, R. T.; Klassen, T. Q.; Valenzano, R. A.; and McIlraith, S. A. 2022.  Reward Machines: Exploiting Reward Function Structure in Reinforcement Learning.  _J. Artif. Intell. Res._ , 73: 173–208.
  * Lee et al. (2022) Lee, J.; Katz, M.; Agravante, D. J.; Liu, M.; Tasse, G. N.; Klinger, T.; and Sohrabi, S. 2022.  Hierarchical Reinforcement Learning with AI Planning Models.
  * Moerland et al. (2023) Moerland, T. M.; Broekens, J.; Plaat, A.; and Jonker, C. M. 2023.  Model-based Reinforcement Learning: A Survey.  _Found. Trends Mach. Learn._ , 16(1): 1–118.
  * Ng and Petrick (2019) Ng, J. H. A.; and Petrick, R. P. A. 2019.  Incremental Learning of Planning Actions in Model-Based Reinforcement Learning.  In Kraus, S., ed., _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019_ , 3195–3201. ijcai.org.
  * Puterman (1990) Puterman, M. L. 1990.  Chapter 8 Markov decision processes.  In _Stochastic Models_ , volume 2 of _Handbooks in Operations Research and Management Science_ , 331–434. Elsevier.
  * Rivlin, Hazan, and Karpas (2020) Rivlin, O.; Hazan, T.; and Karpas, E. 2020.  Generalized Planning With Deep Reinforcement Learning.  _CoRR_ , abs/2005.02305.
  * Seipp, Torralba, and Hoffmann (2022) Seipp, J.; Torralba, Á.; and Hoffmann, J. 2022.  PDDL Generators.  https://doi.org/10.5281/zenodo.6382173.
  * Silver and Chitnis (2020) Silver, T.; and Chitnis, R. 2020.  PDDLGym: Gym Environments from PDDL Problems.  _CoRR_ , abs/2002.06432.
  * Ståhlberg, Bonet, and Geffner (2022a) Ståhlberg, S.; Bonet, B.; and Geffner, H. 2022a.  Learning General Optimal Policies with Graph Neural Networks: Expressive Power, Transparency, and Limits.  In Kumar, A.; Thiébaux, S.; Varakantham, P.; and Yeoh, W., eds., _Proceedings of the Thirty-Second International Conference on Automated Planning and Scheduling, ICAPS 2022, Singapore (virtual), June 13-24, 2022_ , 629–637. AAAI Press.
  * Ståhlberg, Bonet, and Geffner (2022b) Ståhlberg, S.; Bonet, B.; and Geffner, H. 2022b.  Learning Generalized Policies without Supervision Using GNNs.  In Kern-Isberner, G.; Lakemeyer, G.; and Meyer, T., eds., _Proceedings of the 19th International Conference on Principles of Knowledge Representation and Reasoning, KR 2022, Haifa, Israel. July 31 - August 5, 2022_.
  * Sutton and Barto (1998) Sutton, R. S.; and Barto, A. G. 1998.  _Reinforcement learning - an introduction_.  Adaptive computation and machine learning. MIT Press.  ISBN 978-0-262-19398-6.
  * Torreño, Onaindia, and Sapena (2014) Torreño, A.; Onaindia, E.; and Sapena, Ó. 2014.  FMAP: Distributed cooperative multi-agent planning.  _Applied Intelligence_ , 41(2): 606–626.
  * Zhou et al. (2020) Zhou, J.; Cui, G.; Hu, S.; Zhang, Z.; Yang, C.; Liu, Z.; Wang, L.; Li, C.; and Sun, M. 2020.  Graph neural networks: A review of methods and applications.  _AI Open_ , 1: 57–81.


# The spectral dimension of human mobility

Lei Dong1,5†, Kevin O’Keeffe1†, Paolo Santi1,2∗, Mohammad Vazifeh1,
Samuel Anklesaria1, Markus Schläpfer3,4, Geoffrey West4, Carlo Ratti1

1Senseable City Lab, Massachusetts Institute of Technology, Cambridge, MA
02139, USA
2Istituto di Informatica e Telematica del CNR, Pisa 56124, Italy
3Future Cities Lab, ETH Zurich, Zurich 8092, Switzerland
4Santa Fe Institute, Santa Fe, NM 87501, USA
5Institute of Remote Sensing and Geographical Information Systems,
School of Earth and Space Sciences, Peking University, Beijing 100871, China
†These authors contributed equally to this work.
∗To whom correspondence should be addressed<EMAIL_ADDRESS>(P.S.).

> Human mobility patterns are surprisingly structured (?, ?, ?, ?, ?, ?). In
> spite of many hard to model factors, such as climate, culture, and
> socioeconomic opportunities, aggregate migration rates obey a universal,
> parameter-free, ‘radiation’ model (?). Recent work (?) has further shown
> that the detailed spectral decomposition of these flows – defined as the
> number of individuals that visit a given location with frequency $f$ from a
> distance $r$ away – also obeys simple rules, namely, scaling as a universal
> inverse square law in the combination, $rf$. However, this surprising
> regularity, derived on general grounds, has not been explained through
> microscopic mechanisms of individual behavior. Here we confirm this by
> analyzing large-scale cell -phone datasets from three distinct regions and
> show that a direct consequence of this scaling law is that the average
> ‘travel energy’ spent by visitors to a given location is constant across
> space, a finding reminiscent of the well-known travel budget hypothesis of
> human movement (?). The attractivity of different locations, which we define
> by the total number of visits to that location, also admits non-trivial,
> spatially-clustered structure. The observed pattern is consistent with the
> well-known central place theory in urban geography (?), as well as with the
> notion of Weber optimality in spatial economy (?), hinting to a collective
> human capacity of optimizing recurrent movements. We close by proposing a
> simple, microscopic human mobility model which simultaneously captures all
> our empirical findings. Our results have relevance for transportation, urban
> planning, geography, and other disciplines in which a deeper understanding
> of aggregate human mobility is key.

Individuals make regular visits to different places at a wide range of
distance and visiting frequencies. This frequency depends on the type of
activity performed at a destination locationn (eateries, shopping malls, work
places etc) at a certain distance from an individual’s origin place (often an
individual’s home location) (?, ?). In a recent study, we have shown that the
number of visiting individuals follows an inverse square law of the production
of frequency and distance (?). More precisely, we can group the visitors to a
given location $c$ by frequency of visitation $f$ during a reference period
$T$, and consider the spectral flow rates $N_{c,f}(r)$: the total number of
visitors who visit location $c$ from distance $r$ for $f$ times in $T$. The
total number of individuals to $c$ is then $N_{total,c}=\sum_{f}\int
N_{c,f}(r)rdr$. Here, we have computed $N_{c,f}(r)$ using datasets from three
different regions: Greater Boston Area (the United States), Dakar region
(Senegal), and Abidjan (Ivory Coast), see Supplementary Material (SM) and
Table S1 for details.

Following the approach in (?) and defining a high-resolution grid with cells
of size $1km\times 1km$, we construct the user’s movement in two main steps
(see SM for details): $1)$ identify the home cell for each user, which we
define as the grid cell where the user spent the most time at night (see
Fig.1A-C); $2)$ for each (user, cell) pair, compute the number of monthly
visits $f$, and travel distance $r$, from the home cell to the given cell by
the given user, where a cell is considered visited if the user resides there
for a minimum time of $\tau_{\text{min}}=2~{}\text{hours}$. $r$ is defined as
the geographical distance between the center of the user’s home cell and the
center of the visited cell. The desired $N_{c,f}(r)$ are then easily
calculated from the data.

Fig. 1D-F show different frequency groups – hereafter called $f$-groups – have
different flow rates: for fixed travel distance $r$, $N_{c,f}(r)$ declines
with $f$; the frequent visitors to a cell are outnumbered by the infrequent
visitors. Strikingly, under the simplest transformation $r\rightarrow
rf\,(n=1)$, the data collapse to a single, universal curve (?), so that the
visitation density from distance $r$ to a cell $c$, $\rho_{c,f}(r)$, can then
be approximated as $\rho_{c,f}(r)=N_{c,f}(r)/(2\pi
r)={\mu_{c}/2\pi(rf)^{-2}}$, where $\mu_{c}$ is a cell dependent
‘attractivity’ measuring how popular a given cell is (Fig. 1G-I). This tells
us that, in contrast to net migration rates (?, ?) – which the gravity and
radiation models endeavor to explain –, the main parameter governing spectral
flow rates is not the distance $r$ but rather the product $rf$. Since it
measures the total distance traveled by an individual during a given reference
period, we interpret $E:=rf$ as a travel energy (or alternatively, a travel
budget). Our finding, then, is that the common structure between the spectral
flow rates is the travel energy. Or put another way, though their radial
distributions are different, the energy distributions of each frequency group
($f$-group) at a given cell are identical. Hence,
$\rho_{c,f}(r)\propto\mu_{c}/(rf)^{\eta}=\mu_{c}/E^{\eta}$, where $\eta\approx
2$.

A surprising consequence of this finding is that the average travel energy per
visitor to a cell, $\langle E\rangle=E_{total}/N_{total}$, where $E_{total}$
is the total energy spent by visitors to a cell and $N_{total}$ is the total
number of visitors, is spatially invariant – a kind of conservation law of
human mobility:

$\langle E\rangle=\frac{\sum_{f}\int_{r_{min}}^{r_{max}}rfN_{c,f}(r)2\pi
rdr}{\sum_{f}\int_{r_{min}}^{r_{max}}N_{c,f}(r)2\pi
rdr}=\frac{\sum_{f}\int_{r_{min}}^{r_{max}}(rf)^{-1}2\pi
rdr}{\sum_{f}\int_{r_{min}}^{r_{max}}(rf)^{-2}2\pi rdr},$ (1)

where $r_{min},r_{max}$ are the minimum and maximum distances traveled by
walkers in our datasets. We see the only cell dependent quantity, $\mu_{c}$,
cancels out. Fig. 2 shows the conservation law is confirmed by our datasets.
The spatial invariance of $\langle E\rangle$ is surprising because one might
think that more attractive locations in a city would, on average, receive more
travel energy from their visitors. In fact, more attractive places differ only
in the number of visitors they receive, not the travel energy per visitor.

The spatial homogeneity of $\langle E\rangle$ led us to investigate the
spatial distribution of the cell attraction parameters $\mu_{c}$. Recall these
encode how popular, in terms of number of visitors, a given cell $c$ is. Fig.
3A shows $\mu_{c}$ for the Boston dataset have a clustered, spatial structure
where the sizes of the clusters form a hierarchy. The emergence of clusters is
expected: they form from the agglomeration effect of cities, – that is, from
the tendency of services and facilities to locate around city centers or sub-
centers – a finding consistent with the literature on urban structures (?, ?,
?, ?, ?), as well as previous empirical studies of urban mobility (?, ?). The
emergence of the hierarchy of cluster sizes is likely a result of another well
known law of Zipf’s (?). To test this, we investigated if the cluster sizes
are power law distributed. We used the City Clustering Algorithm (CCA) (?) to
compute the clusters from data, which works as follows (see SM for details).
First, the values of all cells with $\mu_{c}$ less than a threshold
$\mu_{c}^{*}$ are set to zero. The values of all remaining cells are set to
$1$. Second, the cells with value $1$ that are contiguous in space are merged
recursively, until ‘islands’ of $1$’s surrounded by $0$’s are formed, giving
the desired set of clusters. Thus, given a threshold $\mu_{c}^{*}$, a set of
clusters is generated. We chose the threshold $\mu_{c}^{*}$, by plotting the
ratio of the area of the largest cluster to the sum of the areas of all the
clusters formed in the Boston data for different $\mu_{c}^{*}$ (Fig. 3C). As
seen, there is a critical value of $\mu_{c}^{*}\approx 10^{2}$ where the area
ratio is minimized; this marks the onset of the emergence of a giant cluster
and serves as a natural choice of $\mu_{c}^{*}$. Fig. 3D shows the
distribution of cluster sizes at this $\mu_{c}^{*}$ do indeed follow Zipf’s
law (?), a law fundamental in city science (?). We show a spatial plot of the
clusters selected at $\mu_{c}^{*}$ in Fig. 3B.

We now take stock of our findings: (i) the universal energy distribution and
its associated conservation law, and (ii) the clustered spatial pattern of
attractivity parameters $\mu_{c}$ whose size distribution match Zipf’s law.
Current models of human mobility cannot simultaneously account for both these
observations. The popular exploration and preferential return model (EPR) (?),
which we will discuss shortly, accounts for (i) but not (ii) (Fig. S4). Here,
building on the EPR model, we develop a model that can produce both (i) and
(ii).

The EPR model is a random walk-like model. At each step with a certain
probability, the walker chooses to explore a previously unvisited location via
a Lévy jump (?), namely, with a radial jump $\Delta r\sim(\Delta
r)^{-1-\alpha}$ and uniformly chosen angle $\theta\sim(2\pi)^{-1}$. If the
walker does not choose to explore she returns to a previously visited location
with a certain probability (A detailed description of the EPR model is given
in SM).

Notice the EPR model describes the motion of a single, independent walker: in
a population of walkers following the EPR model, the individuals do not
interact. In reality, however, individuals’ motions do interact (?): the
motions are correlated through common attraction points and activity hubs.
That is, people do not choose destinations that are entirely independent of
other peoples’ destinations; they tend to visit ‘popular’ places – places
visited frequently by other people. Thus, in ignoring this coupling between
walkers’ motion, the EPR model is unable to reproduce observation (ii): the
clustered distribution of attractivity parameters $\mu_{c}$. As shown in Fig.
S4, the EPR model’s $\mu_{c}$ are uniform across space, in stark contrast to
real data (Fig. 3A).

To account for clustered $\mu_{c}$, we introduce the notion of preferential
exploration, resulting in a modification of the EPR model that we call
preferential exploration and preferential return (PEPR). Preferential
exploration is achieved by coupling the walkers’ motion. When exploring a new
location, a walker is preferentially attracted to popular places, i.e., places
visitors have spent large amounts of energy getting to. The radial jump
distances $\Delta r$ are still sampled from $P(\Delta r)\sim(\Delta
r)^{-1-\alpha}$ but the angle $\theta$ the walker chooses to jump in is no
longer drawn uniformly at random. Instead, angles which correspond to regions
of high visitation are selected preferentially. Let, as before, $E_{total}$ be
the aggregate energy spent getting to the cell by all visitors to that cell.
Further, let the diffused aggregate energy $E_{c}(\theta;R)$ of cell $c$ be
the sum of the aggregate energy of all cells within distance $R$ of $c$
between angles $\theta$ and $\theta+d\theta$. Then walkers following the PEPR
model sample $\theta$ from $P(\theta;R,\nu)\sim E_{c}(\theta;R)^{\nu}$. We
show a schematic of the PEPR model in Fig. 4A.

Figs. 4BC show the PEPR model reproduces finding (i), the spectral flow rates
and their scaling collapse, and more importantly finding (ii), realistic
hierarchical visitation patterns: a qualitatively similar spatial pattern of
clusters (Fig. 4D) and a quantitatively accurate cluster size distribution
(Fig. 4F). Regarding the spatial patterns, we say “qualitatively similar”
since the exact layout of the model clusters is different to that of real
data. For example, in the real data there is a large cluster located on the
coast (corresponding to Boston city) surrounded by multiple smaller clusters,
which is different to the simulation data (Fig. 4D). Reproducing the clustered
spatial patterns at this level of accuracy is however beyond the scope of the
PEPR model since it ignores many complexities which likely influence the
development of human towns/cities such as natural resources, rivers,
topography, etc (see SM). Furthermore, the PEPR model was run on a square
lattice, whereas Boston has an irregular geometry.

Our results support the well-known Central Place Theory (?) of urban science
which to date is (at large-scale) empirically unsupported. The theory asserts
that ‘urban centers’ form an orderly hierarchy arranged in space, where larger
centers, which provide more ‘high-level’ services (e.g., shopping centers,
museums, theaters), are surrounded by smaller centers, which provide ‘local-
level’ services (e.g., groceries, primary schools, clinics). The rationale
behind the theory is that such an arrangement minimizes the total distance
traveled by the population, and is in that sense optimal. Our work
corroborates both aspects of Central Place Theory: the clustered spatial
pattern of $\mu_{c}$ we observed (Fig. 3) is consistent with the hierarchical
structure, and in SM we show the conservation law $\langle E\rangle=const$
across space accords with the minimum-distance optimality. In addition, we
show that the average distance traveled by individuals for a given visiting
frequency $\langle r\rangle_{f}$ obeys the relation $\langle
r\rangle_{f}=K/f$, which also serves as a validation of the Central Place
Theory (Fig. S5).

Central Place Theory is rooted on an individual-level least-effort principle
(?), and an emerging self-organized optimality (?). To strengthen the evidence
for this intriguing possibility, we computed the Fermat-Toricelli Weber (?)
metric of our dataset. This is a metric used in spatial economy to quantify
optimality from the perspective of the activity centers in a city (buildings,
shops etc). Each cell $c$ is assigned an index
$\Delta{\mathcal{D}}_{\text{total}}/\mathcal{D}_{total}\in[0,1]$, where
$\mathcal{D}_{total}$ is the total distance traveled by the reference
population that visits $c$, and $\Delta{\mathcal{D}_{total}}$ is the
improvement in overall distance traveled by the reference population gained by
relocating the destination cell to another position on the grid. If the
location of a cell is already optimal for the reference population,
$\mathcal{D}_{total}$ cannot be reduced by relocating that cell and therefore
the index is $0$. If the location of the cell is suboptimal, the index is
close to 1. Remarkably, Fig. S6 shows most cells in our Boston dataset are
close to their “Weber optimal” locations, having
$\Delta{\mathcal{D}}_{\text{total}}/\mathcal{D}_{total}\approx 0$. We give a
full account of FTW theory and our computations in SM.

This study provides evidences of self-organized optimality of a human
collective behavior, namely, day-to-day mobility. In contrast, many results in
game theory show that collective behavior is non-rational and far from the
socially desired outcome (?, ?). This non-rationality is thought to be due to
cognitive limitations, that is, from the inability of the human mind to
completely understand the complex system in which the human operates (?). The
results of this study stand as a clear counter example to this. They
demonstrate that collectively, humans are able to overcome their cognitive
bounds and achieve optimal group-level behavior – an important and hopeful
finding for the human mind.

## References

  * 1. G. K. Zipf, American Sociological Review 11, 677 (1946).
  * 2. S. Erlander, N. F. Stewart, The gravity model in transportation analysis: theory and extensions, vol. 3 (Vsp, 1990).
  * 3. D. Brockmann, L. Hufnagel, T. Geisel, Nature 439, 462 (2006).
  * 4. M. C. Gonzalez, C. A. Hidalgo, A.-L. Barabási, Nature 453, 779 (2008).
  * 5. C. Song, T. Koren, P. Wang, A.-L. Barabási, Nature Physics 6, 818 (2010).
  * 6. F. Simini, M. C. González, A. Maritan, A.-L. Barabási, Nature 484, 96 (2012).
  * 7. M. Schläpfer, M. Szell, Salat, C. Ratti, G. West, arXiv preprint arXiv:2002.06070  (2020).
  * 8. P. L. Mokhtarian, C. Chen, Transportation Research Part A: Policy and Practice 38, 643 (2004).
  * 9. W. Christaller, Die zentralen Orte in Süddeutschland (Jena: Gustav Fischer, 1933).
  * 10. M. Fujita, P. R. Krugman, A. J. Venables, The Spatial Economy: Cities, Regions, and International Trade (MIT Press, 2001).
  * 11. G. K. Zipf, Human Behavior and the Principle of Least Effort (Addison-Wesley, 1949).
  * 12. A. Anas, R. Arnott, K. A. Small, Journal of Economic Literature 36, 1426 (1998).
  * 13. M. Batty, Science 319, 769 (2008).
  * 14. V. Henderson, J.-F. Thisse, Handbook of Regional and Urban Economics: Cities and Geography, vol. 4 (Elsevier, 2004).
  * 15. A. Bertaud, Order Without Design: How Markets Shape Cities (MIT Press, 2018).
  * 16. T. Louail, et al., Scientific Reports 4, 5276 (2014).
  * 17. C. Zhong, et al., Urban Studies 54, 437 (2017).
  * 18. H. D. Rozenfeld, et al., Proceedings of the National Academy of Sciences 105, 18702 (2008).
  * 19. G. B. West, Scale: the Universal Laws of Growth, Innovation, Sustainability, and the Pace of Life in Organisms, Cities, Economies, and Companies (Penguin, 2017).
  * 20. V. Zaburdaev, S. Denisov, J. Klafter, Reviews of Modern Physics 87, 483 (2015).
  * 21. A. Strandburg-Peshkin, D. R. Farine, I. D. Couzin, M. C. Crofoot, Science 348, 1358 (2015).
  * 22. M. Batty, The New Science of Cities (MIT Press, 2013).
  * 23. A. Weber, C. J. Friedrich, et al. (1929).
  * 24. D. M. Kreps, Game Theory and Economic Modelling (Oxford University Press, 1990).
  * 25. R. B. Myerson, Game Theory (Harvard University Press, 2013).
  * 26. R. Brubaker, The Limits of Rationality (Routledge, 2013).
  * 27. Y.-A. de Montjoye, Z. Smoreda, R. Trinquart, C. Ziemlicki, V. D. Blondel, arXiv preprint arXiv:1407.4885 (2014).
  * 28. V. D. Blondel, et al., arXiv preprint arXiv:1210.0137 (2012).
  * 29. H. D. Rozenfeld, D. Rybski, X. Gabaix, H. A. Makse, American Economic Review 101, 2205 (2011).
  * 30. W. Cao, L. Dong, L. Wu, Y. Liu, arXiv preprint arXiv:1910.12593 (2019).
  * 31. H. Barbosa, et al., Physics Reports 734, 1 (2018).

### Acknowledgments

We thank W.P. Cao for assistance to perform the CCA analysis, and all the
members of the MIT Senseable City Lab Consortium for supporting this research.
L.D. acknowledges funding from the National Natural Science Foundation of
China (No. 41801299).

### Author contributions

C.R., P.S., G.W., L.D., and K.O. designed the research. L.D., K.O., and P.S.
performed the research. L.D., K.O., M.V., S.A., and M.S. analyzed data. L.D.,
K.O., P.S., and M.V. constructed the model and wrote the paper. All authors
reviewed the paper.

### Competing interests

The authors declare no competing interest.

### Data and code availability

The data and code to replicate this research can be requested from the
authors.

Figure 1: Universality in the distance-frequency patterns of human movements.
The home locations for Greater Boston Area (A), Dakar (B), and Abidjan (C). We
show how $N_{c,f}(r)$ is calculated in ((A)). For a given cell, we count
visitors from origin distance within $[r,r+\Delta r]$, see SM for details.
(D-F) The number of visitors $N_{c,f}(r)$ making visits from distance $r$
averaged over a group of cells. Different values of visiting-frequency $f$ are
shown in different colors. (G-I) Re-scaling of the same data with visiting-
frequency, $f$. This confirms the prediction and analysis of ref. (?) which
showed that the visit density for a center, $\rho_{c,f}(r)$, can be well-
approximated by a single function
$\rho_{c,f}(r)=\mu_{c}/(rf)^{-{\eta}},\eta\simeq 2$, implying that the single
parameter, $rf$, is sufficient to express the interplay between distance and
the visiting-frequency, uncovered in ref. (?). Here, data from Abidjan has
been added to further confirm this result ($R^{2}s>0.97$ and standard errors
of $\eta$s are shown in parentheses).

, whi

Figure 2: Constant travel energy per visitor. (A-C) The average energy
$\langle E\rangle$ spent by an individual to visit a cell manifests uniformity
across space consistent with the notion of travel budget per visitor as
discussed in the paper. Note that the southern part of Dakar is an important
port for Senegal, thus a lot of non-local visitors travel to this place,
making the travel distance higher than the remaining places (but still within
the same order of magnitudes). (D) The scatter plots of number of visitors and
travel distance per visitor. The $R^{2}$s of linear regression between number
of visitors and distance per visitor are very small (Greater Boston Area,
$R^{2}=0.0167$, n = 14,273, $p$-value $<$ 0.005; Dakar, $R^{2}<0.001$, n =
173, $p$-value = 0.996; Abidjan, $R^{2}=0.005$, n = 183, $p$-value = 0.355).
Figure 3: Hierarchical structure of attractiveness, $\mu_{c}$. (A)
Geographical pattern of $\mu_{c}$ in Greater Boston Area. We derive cell
specific $\mu_{c}$ by fitting Eq. (1) with the ordinary least squares
regression. We set different thresholds $\mu_{c}^{*}$ for $\mu_{c}$ and then
use the City Clustering Algorithm proposed in (?) to derive the continuous
clusters with $\mu_{c}$ over the threshold (B). We calculate the area ratio of
the area of the largest cluster to the sum of the areas of all clusters (C),
derive the coefficient of the rank-size distributions at the critical value of
$\mu_{c}^{*}\approx 10^{2}$ (vertical dashed line in (C)), and present the
detected clusters in (B) with different colors. When $\mu_{c}$ is very small,
the whole Greater Boston Area would be connected to a single cluster,
resulting in the area ratio $\approx 1$. When $\mu_{c}$ is very large, only
one cluster (Boston downtown) would exist, also resulting in the area ratio
$\approx 1$. (D) Statistical summary of the rank-size regression at the
critical value of $\mu_{c}^{*}$: slope = -1.05 (0.012), $R^{2}=0.975$,
indicating a well-fitted Zipf’s law. Figure 4: PEPR Model and simulation
results. (A) Schematic of the PEPR model. (B-F) Simulation results on a
lattice (with parameters $\alpha=0.55$, $\rho=0.6$, $\gamma=0.21$, $R=10$,
$\nu=4$, and the number of agents = $1\times 10^{5}$). (B) Relations of the
number of visitors $N_{c,f}$ and $r$ with different $f$. (C) Similar to Fig.
1G-I we rescale (B) with visiting-frequency $f$, and all data points collapse
onto a straight line ($\eta\simeq 2$, $R^{2}=0.992$). (D) Attractiveness
$\mu_{c}$ generated by our model shows some significant spatial clusters. (E)
The energy landscape based on the simulation results, which support the
constant energy hypothesis in Eq. (2). (F) We repeat 50 simulations with the
same parameters, and calculate the coefficient of the rank-size distribution
at the critical value of $\mu_{c}$ ($\mu_{c}^{*}=10$). The mean value of the
coefficient is -1.14 (the red dashed line) and the 95% confidential interval
is [-1.39, -0.894] (the black dashed lines), showing a well-fitted Zipf’s law,
which is also similar to the empirical finding (Fig. 3D).

## Supplementary Materials

  * •

Materials and Methods

  * •

Tables S1

  * •

Figures S1-S9

## Materials and Methods

### Boston data

Individuals’ movements in Greater Boston Area are inferred from mobile phone
Call Detailed Records (CDR) data collected over a span of 4 months. The
dataset is provided by a company, and has been used in our previous studies
(?). The raw data contains about 2 million anonymized users.

### Dakar data

The Dakar dataset is based on anonymized Call Detailed Records (CDR) provided
by the Data for Development (D4D) Challenge. The detailed information of this
dataset is provided in (?). Here, we use the SET2, which includes individual
trajectories for 300,000 sampled users in Senegal, and after the preprocess,
we have 173,000 users and 173 cells in Dakar region during two weeks of
January, 2013. We also use the datasets of March, June, and August of 2013 to
verify the robustness of the observed universality of distance-frequency, see
Fig. S3.

### Abidjan data

The Abidjan dataset is also based on anonymized CDR provided by the D4D
Challenge. The structure of the data and the data preprocessing method are
detailed in (?). Here, we use the SET2 of the original dataset. It contains
individual trajectories for 50,000 random sampled users in the Ivory Coast,
and after the preprocess, we have 18,000 users and 183 $1km\times 1km$ cells
in Abidjan during two weeks of December, 2011.

### Data preprocessing

CDR are generated only for voice calls, text messages or data exchanges and
therefore have limited resolution in time. The geographic location of the cell
towers and their density determines the accuracy of location measurements
through triangularization techniques. Therefore, the trajectories extracted
from CDRs constitute a discrete approximation of the moving population
$M(x;y;t)$. There are several steps in preprocessing of the data before it can
be suitable for use in our analysis.

The main steps are: i) Partitioning of the study area. The area under study is
partitioned into a rectangular grid. ii) For each grid cell of size $1km\times
1km$, we identify the individuals that have visited the location with a given
frequency $f$, for instance $f=5$ distinct days in a month for Boston (or bi-
weekly for Dakar and Abidjan), while staying there for a minimum time
$\tau_{min}=2{\text{h}}$. Performing a robustness analysis shows that the
result of our study is not sensitive to small changes in $\tau_{min}$. iii)
For each person, we determine the home location as the grid cell which has
been visited during most nights, i.e. between $7$pm and $7$am of local time.
By summing over all days in a given time window (one month for Boston, and two
weeks for Dakar and Abidjan), one can find the home cell with high level of
confidence for the majority of subjects. The resident population $P_{i}$ of a
given cell $i$ is then defined as the the total number of assigned persons to
that cell. The number of visitors for each cell is defined as the total number
of distinct, non-resident individuals visiting that cell. The number of visits
for each cell is the total number of times that cell has been visited during
the time window of interest. In Fig. S1, we present the visitation
distributions for Boston, Dakar and Abidjan, respectively.

The duration of stay criterium on defining cell visits yields a list of cells
visited by that subject during a day. By aggregating those visits over the
course of a month (or two weeks) for each subject, we obtain a visiting-
frequency vector of dimension $N_{\text{c}ells}$ which is equal to the number
of cells on the geographical grid. The $i$-th component of this vector
represents the number of times the $i$-th cell has been visited by that
subject. We then construct the overall visit matrix for each month
$\mathcal{M}$. The $ij-$component of this matrix is the number of times $j$-th
cell have been visited by the $i$-th subject. Although this matrix is huge in
dimensions, its sparseness allows fast computation to derive various aggregate
mobility related measures.

Here, the distance between cells is calculated by the haversine formula, which
derives the great-circle distance between two points on a sphere. To count the
number of visitors that cell $c$ received from origin distances $[r,r+\Delta
r]$, we take $\Delta r=2km$ for Boston, and $\Delta r=1km$ for Dakar and
Abidjan as the latter two regions are much smaller compared with Boston area.
Meanwhile, to reduce the noise of the ‘tail’ part of the aggregated visit, we
take log-bins for distances over $20km$ in Boston dataset and over $10km$ in
Dakar and Abidjan datasets (Fig. 1 D-F).

### Quantifying spatial structure

We use City Clustering Algorithm (CCA) to derive spatial clusters of cell
attractiveness. CCA, proposed in (?, ?), defines a ‘city’ as a maximal,
spatial continuous area with granular population data. The algorithm takes
three steps: First, set a population threshold $P_{*}$ and binarize the study
area into 0, 1 values – cells with population over $P_{*}$ are set to be 1,
otherwise to be 0. Second, the algorithm picks a populated cell (value = 1)
randomly and adds the nearest populated cells recursively until all the
nearest neighbors are unpopulated cells (value = 0). Third, repeat the picking
and merging process until all populated cells belong to one specific cluster.
This method is intuitive and can divide the US metro area into different
clusters as shown in (?).

In fact, CCA is not limited to use population as the input layer. However, no
matter what kind of input layers used to perform CCA, the common problem is
finding the proper threshold $P_{*}$ to binarize the urban area. A recent
study proposes to employ percolation theory to solve the parameter selecting
problem of CCA (?). The paper has demonstrated that tuning the threshold
$P_{*}$, a giant cluster would emerge as $P_{*}$ reaches a certain point in
datasets of population, nighttime light, and road networks, which is in line
with the two-dimensional percolation process (?). We also find similar
behavior when tuning the threshold of attractiveness $\mu_{c}$ in our case,
which is likely to reflect the self-organization nature of urban systems. By
setting $\mu_{c}$ at the critical value and performing the CCA, we have a
giant cluster and a large number of smaller ones (Fig. 3B). To test the Zipf’s
law, we then run the ordinary least squares (OLS) regression between the
cluster size and its corresponding rank among all detected clusters:

$\log_{10}Size_{i}=\beta_{0}+\beta_{1}\log_{10}Rank_{i}+\epsilon_{i}~{},$ (2)

where $Rank_{i}$ is the size rank of cluster $i$. We derive the parameter of
interest $\beta_{1}$, and report the regression results in the main text (Fig.
3D). Zipf’s law is considered to be a rank-size distribution of
$\beta_{1}=-1$.

### Model and simulation

#### EPR model

The EPR model is a random walk-like model. At each step, the walker decides
whether to explore a new, previously unvisited location with probability
$P_{\text{new}}=\rho S^{-\gamma}$ where $S$ is the number of locations she has
visited so far and $\rho$, $\gamma$ are model parameters. If the walker
decides to explore, she jumps a distance $\Delta r$ sampled from $P(\Delta
r)\propto|\Delta r|^{-1-\alpha}$ – $\alpha$ is another model parameter – at an
angle $\theta$ chosen uniformly at random $P(\theta)=(2\pi)^{-1}$ (i.e., does
a Lévy flight (?), to make the jump sizes consistent with empirical data (?)).
But if the walker does not choose to explore (which occurs with probability
probability $1-P_{\text{new}}$) she returns to a previously visited location
with probability proportional to the number of previous visits in each
location.

#### PEPR model and simulation

We simulated $1\times 10^{5}$ agents moving according to the model’s rules on
a $300\times 300$ grid of cells. Home locations for these agents were assigned
uniformly at random across the grid. Analysis was performed only on the
$100\times 100$ center region to eliminate boundary effects. The model
parameters for data shown in Fig. 4BC were found by strobing over a grid in
parameter space and selecting the parameters which led to the desired scaling
collapse (exponent 2). $\alpha$, $\rho$, and $\gamma$ used here are also
consistent with empirical findings (?). Similarly, those from Fig. 4D-F were
those which led to a cluster size distribution which followed Zipf’s law.

### Weber equilibrium

We here analyze the role places of particular importance play in the formation
of self-organized patterns of urban settlements by analyzing the efficiency of
these cells from the point of view of spatial-economic theory. According to
this theory, the urban population distribution is driven by centrifugal and
centripetal forces which exist due to the economic competition in minimizing
the transportation costs to important resources. We can investigate the
efficiency of attractor cells quantitatively, by studying how close the total
transportation cost of incoming visits to each cell is to the optimum
transportation cost defined as the minimum possible value for the total
transportation distance from visitor’s home-location. This problem can be
formulated as the Fermat-Torricelli-Weber (FTW) problem on a square grid.

We define a bi-directed visit OD-flow spatial network in which the nodes
correspond to geographical cells and the directed edges are weighted according
to the number of visits exchanged between pairs of nodes. The visits flow
matrix is an asymmetric square matrix which contains all the information about
the visiting patterns and is defined as

$\mathcal{V}_{ij}=\text{ total visits
from}\;\mathcal{C}_{i}\;\text{to}\;\mathcal{C}_{j}$ (3)

Note that in general $\mathcal{V}_{ij}\neq\mathcal{V}_{ji}$.

In the Weber problem in location theory, the optimal point is a point which
minimizes the total distance from $n$ points on a plane – Fig. S6. One can
consider this problem on a grid where the optimum location can be chosen from
a finite number of points corresponding to the centre of cells on a
geographical grid, and the optimum is a cell which minimizes the overall
transportation distances from where the visits originate from.

We define the Weber matrix as follows:

$\mathcal{W}_{ij}=T[\mathcal{C}_{i}\rightarrow\mathcal{C}_{j}]$ (4)

where $T[\mathcal{C}_{i}\rightarrow\mathcal{C}_{j}]$ is the total distance
travelled by visitors of $i$-th cell if this cell where placed in the location
of $j$-th cell. Using the distance matrix $\mathcal{D}$ and the visits flow
matrix $\mathcal{V}$ we can compute the Weber matrix

$\mathcal{W}_{ij}=\sum_{k}\mathcal{D}_{jk}\mathcal{V}_{ki}=[{\mathcal{D}\cdot\mathcal{V}}]_{ji}$
(5)

Each row of the Weber matrix contains all the possible values of the objective
function defined according to the Weber problem on a grid for the
corresponding cell. The question is how close the value of the actual total
transportation distance for each cell, which correspond to the value on the
diagonal axis of the matrix, is to the minimum possible value for each row,
corresponding to the FTW cell. One way to measure this closeness is to see how
much improvement can be gained for each cell if we move each cell to its FTW
location. We define the fractional improvement as the ratio of the total
energy improvement gained for each cell to the actual energy which is given as

$\frac{\Delta\mathcal{D}_{i}}{\mathcal{D}_{i}}=\frac{\mathcal{W}_{ii}-\text{min}({\mathcal{W}_{i*}})}{\mathcal{W}_{ii}}$
(6)

The above quantity is always between zero and one. The value zero corresponds
to the extremum case where no improvement can be gained, meaning that the
cell’s location coincides with the optimal transportation location. The value
would be equal to one when the transportation distance can be reduced to zero
by moving the cell. Since the number of visits a cell gets does not change as
we relocate the cell, the fractional distance per visit improvement, i.e.,
$\frac{\Delta\mathcal{D}_{iv}}{\mathcal{D}_{iv}}$ is equal to the fractional
total distance improvement,

$\frac{\Delta\mathcal{D}^{\text{per visit}}_{i}}{\mathcal{D}^{\text{per
visit}}_{i}}=\frac{\Delta\mathcal{D}_{i}}{\mathcal{D}_{i}}$ (7)

The average distance per visit can be quite large yet the highly important
cells are very close to their FTW cell in the majority of the cases. In Fig.
S6 we plot the total received visits by each cell versus the fractional
improvement which can be gained by relocating them to their FWT point. As seen
from the figure, the majority of the highly visited cells have low fractional
improvement. The exceptions to this pattern are the few cells in the yellow
box in Fig. S7. To see why these cells were exceptional, we checked the
location of the cells on the map and found that in the majority of cases, they
correspond to tourism attraction points near beaches, lakes, etc, which
explains why they are anomalous – these locations having an intrinsic reason
to be located where there are, as opposed to being there so as to optimize
their FTW score.

### Topography

In the main text we showed that the spatial pattern of the visitation rates of
the PEPR model (Fig. 4D) were different to those of the real Boston data (Fig.
3A). This mismatch was not surprising since the PEPR model only models how the
interactions between individuals influences their movements and is blind to
the terrain on which the individuals move. Presumably, different types of
terrain would attract / repel people with different strengths. For example,
areas with lots of natural resources would naturally attract settlement, as
would rivers and coasts be attractive since they influence trade. Put
generally, human movement would be influenced by topography, an effect which
the PEPR model does not strive to capture.

A thorough study of the role of topography in human movement is beyond the
scope of the present work. We here however take a first step in this direction
by running the PEPR model on a non-trivial geometry to see if it leads to more
realistic spatial visitation patterns. In the main text, the PEPR model was
run on a square lattice, a crude approximation of the irregular geometry of
Boston (Fig. 3A). In Fig. S8 we show the visitation pattern of the PEPR model
when run on a lattice with a simple perturbation: a rectangular chunk removed
from one side. As seen, the spatial visitation pattern is not qualitatively
altered, demonstrating that other topographical features are needed to recover
the hierarchical pattern observed in real data (Fig. 4D main text).

Region | Country | # of CDR users | CDR date | Population | Area
---|---|---|---|---|---
|  | after preprocess |  | $million$ | $km^{2}$
Greater Boston Area | US | 340,000 | 2009 | 4.73 | 11,700
Dakar Region | Senegal | 173,000 | 2013 | 2.96 | 547
Abidjan | Ivory Coast | 18,000 | 2011-12 | 4.70 | 422

Table S1 $|$ Statistical summary of three regions. Population of Abidjan was
derived from 2014 census of Ivory Coast. Population and area data of Dakar
region were derived from 2013 census of Senegal. Population and area data of
Greater Boston Area were derived from Wikipedia.

Fig. S1 $|$ Geographical distributions of the total visitation. A, Boston; B,
Dakar; and C, Abidjan.

Fig. S2 $|$ The cumulative distribution of the number of visitors visiting
from distance within $r$ radius for various visiting frequencies $f$. A,
Boston; B, Dakar; and C, Abidjan. The curves for different frequencies
collapse approximately into a single curve after rescaling the distance with
frequency, $r\rightarrow rf$.

Fig. S3 $|$ Robustness check of the universal $rf$ in different time period.
Dakar datasets of A, March, 2013; B, June, 2013; and C, August, 2013. The
variation density for a wide range of $rf$ can be well-approximated by a
single function $\rho_{c,f}(r)=\mu_{c}\cdot(rf)^{-{\eta}},\eta\simeq 2$
($R^{2}>0.98$).

Fig. S4 $|$ Simulation results of the EPR model. Visitations (A) and
attractivity parameters $\mu_{c}$ (B) generated by the EPR model are uniform
across space, which is in contrast to real data (Fig. 3 and Fig. S1). C, Dd,
Similar to Fig. 4B, C, EPR model can reproduce Eq.(1).

Fig. S5 $|$ $\langle r\rangle=K/f$ and the Central Place Theory. A, Schematic
figure of the Central Place Theory. The spatial arrangement of three Tiers of
centers in a two dimentional space. This hierarchical arrangement of central
places results in the most efficient transport network. C-D, The distance
travelled per visit to perform activities with visiting-frequency $f$ averaged
over all the individuals: Boston (B), Dakar (C), and Abidjan (D). $\langle
r\rangle=K/f$ curve fits very well with the empirical observation and supports
the notion of universal travel-budget.

Fig. S6 $|$ Fermat-Torricelli-Weber (FTW) efficiency of collective human
movements.. A) The schematic figure shows how the FTW efficiency is computed.
The total distance travelled by visitors of a specific cell (red dot) can be
minimized by moving the destination cell on the grid. The efficiency is
$\Delta{\mathcal{D}_{total}}/{\mathcal{D}_{total}}$, which is the ratio
between $\Delta{\mathcal{D}_{total}}$, i.e., the improvement gained in
reducing the aggregate distance travelled by moving the cell from its actual
location to the optimum FTW point, and the actual aggregate distance travelled
by visitors to that cell, ${\mathcal{D}_{total}}$. B) Each density plot
represents the number of cells with a particular number of visits and FTW
efficiency for Greater Boston Area based on CDR for the month of August 2009.
The FTW efficiency is computed for each cell based on visits made by visitors
who live at distances larger than $r_{c}$. These plots compare how density
changes by increasing $r_{c}$ from 0 to $10~{}km$. For $r_{c}=0$ the density
is particularly high where the FTW efficiency is very high. As the number of
visits is increased, the distribution becomes narrower and the FTW efficiency
increases. This pattern still survives but becomes weaker as $r_{c}$ is
increased as described in Supplementary Material.

Fig. S7 $|$ Transportation optimality of centers and tourism outliers for
Greater Boston Area based on CDR. A, The fractional distance improvement
versus cell’s number of visits. The higher the number of visitors, the higher
is the chance that the cell is transportation efficient. As shown in B, C, the
outliers to this pattern (the cells corresponding to points in the yellow box
in (a)) are located near shores, lakes, etc., and are well-known touristic
locations.

Fig. S8 $|$ Visitation pattern of PEPR model on non-square lattice.