# A Novel Framework for Decentralized Dynamic Resource Allocation Using
Voronoi Tessellations

Bhagyashri Telsang and Seddik Djouadi This paper was supported in part by the
National Science Foundation under grant NSF-CMMI-2024111.B. Telsang and S.
Djouadi are with the Department of Electrical Engineering and Computer
Science, University of Tennessee Knoxville, USA<EMAIL_ADDRESS>
###### Abstract

In this work, we approach the problem of resource allocation in a team of
agents through the framework of Centroidal Voronoi Tessellations. CVTs provide
a natural way to embed a desired global trend in the team through probability
distributions, and in one-dimensional spaces, CVTs offer an inherent line
structure allowing for a simple communication graph and scalability. We first
consider the amount of resource to be allocated to be a constant and provide
an analytical solution to such static resource allocation problem by embedding
the allocation constraint within the distribution through a system of
nonlinear equations. Using the solution of such a constrained CVT minimization
problem as an initialization step, we propose a decentralized dynamic resource
allocation solution that employs a one-step update when the desired
distribution is Gaussian. We introduce a “civility model” for negotiations
between the agents to allow for flexibility in local preferences and
maintaining robustness against local disturbances. We demonstrate the
effectiveness of the proposed method by considering the application of demand-
response in smart grids through the problem of power allocation in a group of
building thermal loads.

##  I INTRODUCTION

Often times we see a conflicting, paradoxical problem around us in the world.
There is too much and yet there is not enough; obesity and hunger coexisting,
overpopulation and population scarcity coexisting, floods and drought
coexisting few hundred miles apart, vacant house and homeless people outside
of them. In each of the scenario, there is a resource that is abundant in one
sector, be it location or a group of people, but scarce in another. It makes
one wonder if we can alleviate the problem by allocating the resources in a
“right” manner.

Taking roots in the field of Economics through [1] in $1970^{\prime}s$, the
resource allocation problem has broadened to the field of engineering in more
recent decades. Mathematically, the resource allocation problem can be framed
as, [2]:

$\displaystyle\min_{z_{i}\in\mathbb{R}^{n}}\frac{1}{N}\sum_{i\in
I_{N}}f_{i}(z_{i})$ $\displaystyle\text{such that,}\hskip 14.22636pt\sum_{i\in
I_{N}}z_{i}=r$ (1)

In the resource allocation problem, an $r$ amount of resource is to be
allocated among $N$ agents while minimizing the sum of their individual cost
functions $\\{f_{i}\\}_{i\in I_{N}}$. Typically, in engineering problems, the
agents are local controllers tasked to maintain local interests while equipped
with capabilities to communicate with other agents.

Simultaneously seeming trivial and complex, the nature of (1) can be broken
down into the following aspects: the information structure in the group of
agents, the separability of the objective function, and the global constraint.
Due to the separability of the objective function, each agent can minimize the
(global) cost function without any dependance on other agents. However,
because of the global constraint imposed on the team, the team information
structure becomes a significant aspect.

Like most of the work on the resource allocation problem, the authors in [2]
assume the individual cost functions to be convex. In the case where the cost
functions are differentiable, they propose a gradient descent consensus
algorithm. And when the cost functions are not necessarily differentiable,
they present a sub-gradient based algorithm. While they let the team
information structure be dynamic, they impose reasonable mild conditions on
the team information structure like connectedness, and start at an initial
feasible condition.

While [2] proposes the gradient descent algorithm where the agents trade
resources in proportion to the gradient difference for their individual cost
functions, [3] takes up the allocation problem (1) to focus on choosing the
proportional weights (to the resource trading) to obtain sufficient conditions
for the convergence of the algorithm, and to further improve the rate of
convergence. [4] considers the dual of the resource allocation problem and
derives two methods using the alternating direction method of multipliers
(ADMM) algorithm. Also considering the dual problem, [5] includes
uncertainties in the individual cost functions and solves the problem using
sub-gradient methods on the distributed Lagrangian.

Mixing economics in the team, [6] considers a stochastic system in which
agents allocate shared system resources in response to customer requests that
arrive stochastically over time and introduces the notion of a transfer
contract to specify compensation among agents when resources are traded. Each
agent has a model of how resources are shared by others and it makes
allocation decisions by maximizing its utility function subject to such model.
However, it is worth noting that in most of the work on decentralized resource
allocation, the amount of resource to be allocated is fixed over the
iterations; the agents begin at a feasible solution and move along the
resource allocation constraint through the feasible solutions to only minimize
the cost function in (1).

In this work, we approach the resource allocation problem through Centroidal
Voronoi Tessellations (CVTs). Even though they date centuries, Voronoi
tessellations (VTs) have been found to be immensely helpful in various
applications ranging from health to computer graphics to natural sciences. The
first documented application of Voronoi tessellations appeared in [7] on the
1854 cholera epidemic in London in which it is demonstrated that proximity to
a particular well was strongly correlated to deaths due to the disease [8]. In
more recent decades, VTs have almost become a common basis tool for path
planning algorithms by multi-robot systems in the field of coverage control
[9] to such an extent that the VT-based coverage control has been generalized
using optimal transport-based control [10]. An adaptive coverage controller is
proposed in [11] where the leader in the leader-follower strategy therein
distributes the followers within its obstacle-free sensing range, and the
optimized distribution is obtained through a CVT. In their study on optimality
of multi-robot coverage control, the authors in [12] draw a relationship
between CVT configurations and the sufficient condition for optimality through
the spatial derivative of the density.

Despite the wide range of applications of CVTs, to the best of our knowledge,
they have not been employed in the resource allocation problem. Our main
motivations for employing them in the context of resource allocation problem
are minimal communication requirement, robustness, flexibility, scalability
and generalizability offered by the CVT framework. To allocate one-dimensional
resources, a line graph for the communication network in the team information
structure is sufficient to obtain the global optima. We will delve deeper into
the advantages and our motivation for using CVT in the resource allocation
problem in Section II.

To demonstrate our solution to the resource allocation problem using the CVT
framework, we consider the application of demand-response in smart grids. We
are in the era of a volatile energy market with a looming energy crisis. While
the underlying occurrences like geopolitical transitions that cause such
crises are beyond local control, the effects are certainly received through
the spectrum. In such cases, what one can do locally, on smaller scales is to
better employ the resources available at hand. The field of demand-response in
smart grids aims to maintain robust operation during such times.

As many parts of the world are gradually moving towards competitive
transactive energy markets as a means to generate and procure electricity
alongside many of the support services required to operate a power system,
many countries are pushing the reform of the electricity power sector very
positively. For example, Chile pioneered in the 1980s the deregulation of the
electric power industry. In today’s U.S. retail electricity market, fourteen
states have already adequate retail competition with Texas, Illinois, and Ohio
respectively having 100$\%$, 60$\%$, and 50$\%$ of their residential customers
receiving service from electricity suppliers [13]. But, even today, most of
the customers have very limited “direct” participation in supporting the grid.

Through the developments in the transactive energy market, there have been
some interesting and innovative proposals. [14] proposes a data-driven method
to forecast the electricity demand for decentralized energy management. On the
consumption end, although mostly work in progress, peer-to-peer (P2P)
electricity trading is gaining momentum with time. Analogous to internet
servers and clients, P2P electricity trading is the platform where the end
consumer becomes a prosumer (functioning as both energy producer and consumer)
and exchanges the remaining electricity with other consumers in the power grid
[15]. A detailed review of existing P2P trading projects is carried out in
[16]. Another major proposal in this direction is load aggregation. As defined
in [17] “An aggregator is a new type of energy service provider which can
increase or moderate the electricity consumption of a group of consumers
according to the total electricity demand on the grid. An aggregator can also
operate on behalf of a group of consumers producing their own electricity by
selling the excess electricity they produce.” A detailed review of the value
of aggregators in the electricity market can be found in [18].

While such a centralized framework is beneficial in certain applications the
cost and the risk of the associated communication overhead can be too high.
Such a center-heavy approach also makes the framework vulnerable to attacks
due to a single point of failure. The advantages of having communication
capabilities between agents reduces such associated risks. One then needs to
develop a framework to model the flow of information among the agents and
design control laws at global and local levels such that the global and local
objectives are achieved. We consider a certain power, generated or negotiated
between the aggregator and the utilities, as the resource to be allocated
among a team of agents that are building loads, HVACs to be specific. Such
thermal loads, due to their latency, inherently allow for further flexibility
in the CVT framework.

The paper organization is as follows. We begin with review of some definitions
and preliminaries, along with our motivation for employing CVTs to solve this
problem in Section II. Like most work on resource allocation problems, we
consider the static allocation problem where the amount of resource to be
allocated is fixed, and solve it using a system of non-linear equations in
Section III. We then move to varying the allocation amount and solve the
corresponding dynamic resource allocation problem in Section IV. In Section V
we demonstrate the developed decentralized dynamic resource allocation method
on a demand-response problem of power allocation in a group of building loads.
Finally, we draw conclusions in Section VI and present some lines of future
work.

##  II Preliminaries

In this Section, we will first review some definitions and background of CVTs
in Section II-A, followed by a brief review on iterative and analytical
methods to compute CVTs in Section II-B. Then in Section II-C, we define some
notations on communication and resource graph employed in this paper.

###  II-A Centroidal Voronoi Tessellations

Consider a region $\Omega\in\mathbb{R}^{n}$ with density $\rho(.)$. Denote a
team of $N\in\mathbb{N}$ agents indexed by the set $I_{N}=\\{1,2,\ldots,N\\}$.

  * 1

Tessellation: $\\{V_{i}\\}_{i\in I_{N}}$ is a tessellation of $\Omega$ if
$V_{i}\cap V_{j}=\emptyset$ for $i\neq j$, and $\cup_{i\in
I}{V}_{i}={\Omega}$.

  * 2

Voronoi region and generators: The Voronoi region ${V}_{z_{i}}$ of the Voronoi
generator $z_{i}$ is ${V}_{z_{i}}=\\{x\in\Omega:||x-z_{i}||<||x-z_{j}||,\
i\neq j\ \text{and}\ i,j\in I_{N}\\}$.

  * 3

Voronoi tessellation: The set of Voronoi regions
$\textbf{V}_{\textbf{z}}=\\{V_{z_{i}}\\}_{i\in I_{N}}$ of $\\{z_{i}\\}_{i\in
I_{N}}$ is called a Voronoi tessellation
$\\{\textbf{z},\textbf{V}_{\textbf{z}}\\}$.

The mass centroid of a region $V_{i}\subset\Omega$ under the probability
density function $\rho(.)$ is defined as:

$z_{V_{i},\rho}^{c}=\frac{\int_{V_{i}}x\rho(x)dx}{\int_{V_{i}}\rho(x)dx}$ (2)

A Voronoi tessellation in which the generators are the mass centroids of their
respective Voronoi regions is called a Centroidal Voronoi Tessellation (CVT),
[19]. The CVT obtained for 3 generators in the region $\Omega=[0,15]$ under
Uniform and Normal distributions – $\mathcal{U}(0,15)$ and
$\mathcal{N}(7.5,1)$ – are shown in Fig. 1. The generators under Uniform ad
Normal distribution over $\Omega$ are marked in star and square symbols
respectively.

Figure 1: Centroidal Voronoi Tessellations of $[0,15]$ under Uniform and
Normal distributions, denoted in star and square symbols respectively.

Consider the functional $\mathcal{F}$ with any $N$ points $\\{z_{i}\\}_{i\in
I_{N}}\in\Omega$ and any tessellation $\\{V_{i}\\}_{i\in I_{N}}$ of $\Omega$
as its input arguments:

$\mathcal{F}((z_{i},V_{i}),i\in I_{N})=\sum_{i\in I_{N}}\int_{x\in
V_{i}}||x-z_{i}||^{2}\rho(x)dx$ (3)

Proposition $3.1$ in [19] states that a necessary condition for the function
$\mathcal{F}$ to be minimized is that $\\{V_{i}\\}_{i\in I_{N}}$ are the
Voronoi regions corresponding to $\\{z_{i}\\}_{i\in I_{N}}$, and
simultaneously, $\\{z_{i}\\}_{i\in I_{N}}$ are the centroids of their
respective Voronoi regions. In other words, the minimizer of $\mathcal{F}$ is
a Centroidal Voronoi Tessellation.

Additionally, if the tessellation in (3) is fixed to be the Voronoi
tessellation of $\\{z_{i}\\}_{i\in I_{N}}$, then the following functional
$\mathcal{K}$ has the same minimizer as $\mathcal{F}$, [19].

$\mathcal{K}((z_{i}),i\in I_{N})=\sum_{i\in I_{N}}\int_{x\in
V_{z_{i}}}||x-z_{i}||^{2}\rho(x)dx$ (4)

This functional $\mathcal{K}$ is also referred to as the energy of the
tessellation or the quantization energy.

Including the resource allocation constraint from (1) in the functional
$\mathcal{K}$, we obtain the following constrained CVT minimization problem:

$\displaystyle\min_{z_{i}}\sum_{i\in I_{N}}\int_{x\in
V_{i}}||x-z_{i}||^{2}\rho(x)dx$ s.t, $\displaystyle\hskip 28.45274pt\sum_{i\in
I_{N}}z_{i}=r$ (5)

Comparing with the resource allocation problem (1) we see that the individual
objective functions from (5) are $\\{\int_{x\in
V_{i}}\rho(x)||x-z_{i}||^{2}dx\\}_{i\in I_{N}}$. While the objective function
in (5) is separable, a global distribution $\rho(.)$ governs all the
individual objective functions. This enables embedding of a desired aggregate
behavior in the team through such distributions; the desired aggregate
behavior can arise from modeling individual preferences or from an external
global trendsetting factor depending on the application at hand.

###  II-B Computation of CVT

Given $\Omega,N$ and $\rho(.)$, there are various iterative algorithms to
compute a CVT in $\Omega$. Under the same conditions, CVT need not be unique
for any dimensional region unless certain conditions are imposed on the
density function. In 1-D regions, the CVT is unique for log-concave density
functions with finite second moment [20]. For higher dimensions, finding the
conditions on the uniqueness for the general case, without assumptions on the
region, density or the number of generators $N$, remains an open area of
research. However, it is proved in [21] that for $N=2$, there does not exist a
unique CVT for any density for dimensions greater than one.

Accordingly, the solutions rendered by various algorithms to compute the CVT
need not be the unique global minimizers, and can converge to a local minima.
A deterministic, popular algorithm to obtain a CVT is the Lloyd’s algorithm.
Introduced in [22] to find the optimal quantization in pulse-code modulation,
Lloyd’s algorithm has been modified or adapted in various fields. At the core
of it, Lloyd’s algorithm is an iteration between constructing Voronoi
tessellations and their centroids:

Given: $\Omega\subset\mathbb{R}^{n}$, $N$, $\rho(x)$
Initialize: Generators $\textbf{z}=\\{z_{i}\\}_{i\in I}$, where each
$z_{i}\in\Omega$

  * 1

Construct the Voronoi tessellation $\textbf{V}_{\textbf{z}}$.

  * 2

Compute the mass centroids $z^{c}_{\textbf{V}_{\textbf{z},\rho}}$ of
$\textbf{V}_{\textbf{z}}$.

  * 3

If the computed centroids meet certain stopping criteria then terminate. If
not, then set $\textbf{z}=z^{c}_{\textbf{V}_{\textbf{z},\rho}}$, and return to
Step 1.

Even though Lloyd’s algorithm is iterative and approximate, it has certain
desirable convergence properties. Various global convergence properties of the
Lloyd’s algorithm are rigorously proved in [23]. Specifically for one-
dimensional spaces with log-concave density function, the local convergence
using the Lloyd’s algorithm has been proved in [24]. Depending on the
application at hand, various algorithms that have faster convergence than
Lloyd’s have been proposed, [25], [26], [27]. Taking the probabilistic
approach, [28] offers a Monte Carlo sampling based method for the computation
of CVTs in any dimension. However, the MacQueen’s method only results in the
centroids of the CVT and not their Voronoi partitions.

In one-dimensional spaces, we can obtain the entire tessellation analytically
using a System of Non-linear Equations (SNLE). The core idea is to
parameterize the Voronoi regions in terms of their centroids. In
$\Omega=[a,b]\subset\mathbb{R}$, without loss of generality, let the $N$
generators be $z_{1}<z_{2}<\ldots<z_{N}\in\Omega$. Following Section II-A, by
definition, the Voronoi regions are given as:

$\displaystyle V_{1}$ $\displaystyle=[a,\frac{z_{1}+z_{2}}{2}]$
$\displaystyle\vdots$ $\displaystyle V_{i}$
$\displaystyle=[\frac{z_{i-1}+z_{i}}{2},\frac{z_{i}+z_{i+1}}{2}]$
$\displaystyle\vdots$ $\displaystyle V_{N}$
$\displaystyle=[\frac{z_{N-1}+z_{N}}{2},b]$ (6)

Additionally, by the definition of CVT the Voronoi generators must be the mass
centroids (2). Rewriting the centroids in terms of the parameterized Voronoi
regions from (6), we have $\forall i\in I_{N}$:

$\displaystyle z_{i}^{c}=$
$\displaystyle\frac{\int_{\frac{z_{i-1}^{c}+z_{i}^{c}}{2}}^{\frac{z_{i}^{c}+z_{i+1}^{c}}{2}}x\rho(x)dx}{\int_{\frac{z_{i-1}^{c}+z_{i}^{c}}{2}}^{\frac{z_{i}^{c}+z_{i+1}^{c}}{2}}\rho(x)dx}$
(7)

where $z_{V_{i},\rho}^{c}$ from (2) is denoted as $z_{i}^{c}$ for ease of
notation. In (7), there are $N$ number of unknowns: $\\{z_{i}^{c}\\}_{i\in
I_{N}}$, and $N$ equations. Therefore, solving this system of nonlinear
equations will result in the centroids of the CVT with which the Voronoi
regions can be computed as in (6). This is the exact solution of the
functional $\mathcal{K}$ from (4).

To employ CVTs to solve the resource allocation problem (1), we treat the
Voronoi generators as the agents’ resources. In the next Section, we introduce
the notations for resource and communication graph for the team which
highlights the advantage of employing one-dimensional CVTs in the resource
allocation problem.

###  II-C Notations

Let $\textbf{z}=\\{z_{i}\\}_{i\in I_{N}}$ be the agents’ resources with each
$z_{i}\in\Omega\subset\mathbb{R}$. Let $\rho(.)$ denote a measure of
information or the probability density over $\Omega$.

Let $\mathcal{Z}$ denote the undirected resource graph and $\mathcal{C}$
denote the undirected communication network of all the agents $i\in I_{N}$.
Their vertex and edge sets are $\\{\textbf{z},\mathcal{E}_{Z}\\}$ and
$\\{I_{N},\mathcal{E}_{C}\\}$, respectively. Denote the set of neighbors of
agent $i\in I_{N}$ according to resource and communication graphs as
$\mathcal{N}_{\mathcal{Z}_{i}}$ and $\mathcal{N}_{\mathcal{C}_{i}}$,
respectively.

The set of resource neighbors $\mathcal{N}_{\mathcal{Z}_{i}}$ of each agent
$i\in I_{N}$ is given by [29]:

$\displaystyle\mathcal{N}_{\mathcal{Z}_{i}}=\\{j\in I_{N}:z_{k}<z_{j}<z_{i},\
\forall j,k\in I_{N}\\}\ \ \cup$ $\displaystyle\\{j\in
I_{N}:z_{i}<z_{j}<z_{k},\ \forall j,k\in I_{N}\\}$ (8) $\displaystyle\implies\
\ j\in\mathcal{N}_{\mathcal{Z}_{i}}\iff\\{z_{i},z_{j}\\}\in\mathcal{E}_{Z}\iff
i\in\mathcal{N}_{\mathcal{Z}_{j}}$

Parallelly, since $\mathcal{C}$ is undirected,
$\\{i,j\\}\in\mathcal{E}_{C}\iff i\in\mathcal{N}_{\mathcal{C}_{j}}$ and
$j\in\mathcal{N}_{\mathcal{C}_{i}}$, where $\mathcal{N}_{\mathcal{C}_{i}}$ is
the set of communication neighbors of the agent $i$.

Since the resources $z_{i}\in\mathbb{R}$, the resource network is always a
line graph: each agent can have at most $2$ resource neighbors. While the
communication network $\mathcal{C}$ can be as complex as a full graph, we set
it to be the simplest connected graph in 1-D,
$\mathcal{E}_{C}=\mathcal{E}_{Z}$. That is, the agents are aware of the
resource positions of their neighbors only.

Any agent being informed of the resource positions of its non-neighbor agent
is redundant, since it is not used in the iterative CVT computation methods
like Lloyd’s. Therefore, if each agent were to communicate only with its
resource neighbors, then all the agents would converge to the CVT through
Lloyd’s algorithm with minimal communication in a decentralized manner.

In this Section, we looked into how CVTs provide a natural way of embedding a
desired distribution in the solution, along with obtaining the solution in a
straight-forward decentralized approach with minimal requirements on the team
information structure. We studied that one may obtain 1-D CVTs in a
decentralized manner using one of the simplest communication graphs: a line
graph that is also the same as the resource graph. In the next Section, we
take up the resource allocation problem (1) which is a constrained CVT
minimization problem, and solve it centrally using the analytical CVT
computation method SNLE.

##  III Static resource allocation

The underlying idea in employing CVTs to solve the resource allocation problem
is quite straightforward: the CVT centroids are the resources allocated to the
agents, and accordingly, they must sum up to the available amount of resource
$r$. Comparing the resource allocation constrained CVT minimization problem
(5) with (1), we can observe that the individual agent cost functions are:

$f_{i}(z_{i})=\int_{x\in V_{i}}||x-z_{i}||^{2}\rho(x)dx$ (9)

The main similarity between the constrained CVT minimization problem (5) and
the resource allocation problem (1) is that the objective functions are
separable. That is, the objective functions are decomposed into individual
(agent) objective functions that are convex and differentiable. Additionally,
the integral equation (9) is known as the Fredholm integral equation of the
first kind, [30].

Given the separability of (5) along with the convexity and differentiability
of the agent cost functions, the asymptotic convergence properties developed
in [2] for the resource allocation problem apply to the resource allocation
constrained CVT minimization problem (5). The core idea of our solution to the
problem of resource allocation through CVTs is to embed the resource
allocation constraint within the objective function through the density
$\rho(.)$.

Suppose $\rho(.)$ is defined by $n_{\rho}$ number of parameters:
$v=(v_{1},v_{2},\ldots,v_{n_{\rho}})$. Let $v_{k}\in v$ for some $k\in
I_{N_{\rho}}$ be an unknown or the “free” design parameter, and all the other
parameters defining the density function be known and fixed. To highlight the
dependence of the density function on the free parameter $v_{k}$, denote the
density function as $\rho(x,v_{k})$, where $x$ is it’s support.

The optimal solution of the unconstrained CVT minimization problem (4) is the
set of centroids of the Voronoi regions for every agent. Using the definition
of centroids (2) for $\\{z_{i}\\}_{i\in I_{N}}$ and embedding the resource
allocation constraint transforms the constrained CVT minimization problem into
the following system of nonlinear equations with $N+1$ unknowns –
$(z_{1}^{c},z_{2}^{c},\ldots,z_{N}^{c},v_{k})$:

$\displaystyle z_{i}^{c}=$
$\displaystyle\frac{\int_{\frac{z_{i-1}^{c}+z_{i}^{c}}{2}}^{\frac{z_{i}^{c}+z_{i+1}^{c}}{2}}x\rho(x,v_{k})dx}{\int_{\frac{z_{i-1}^{c}+z_{i}^{c}}{2}}^{\frac{z_{i}^{c}+z_{i+1}^{c}}{2}}\rho(x,v_{k})dx}\
\ \ \forall i\in I_{N}$ $\displaystyle\sum_{i=1}^{N}z_{i}^{c}=r$ (10)

The solution of this system of nonlinear equations,
$(z_{1}^{c},z_{2}^{c},\ldots,z_{N}^{c},v_{k})$, satisfies the following:

  * •

$(z_{1}^{c},z_{2}^{c},\ldots,z_{N}^{c})$ are the $N$ centroids of the CVT in
$\Omega=[a,b]$ under the density function $\rho(x,v_{k})$.

  * •

The centroids sum up to $r$, satisfying the resource allocation constraint in
(1).

The main solution of interest here is the solved design parameter $v_{k}$
which is fed to the Lloyd’s algorithm in its initialization step. In that
case, all the agents can still maintain communication only with their resource
neighbors, and since they are all initialized with the same design parameters,
all the agents are minimizing the cost function (4) under the same
specifications, and obtain the CVT.

We now demonstrate the method with different simulation cases. In Fig. 2, the
region $\Omega=[0,100],\ N=50$, and the density is Gaussian. Out of the three
examples therein, the top two have the same variance but are required to
allocate different amounts of resources – $2500$ in the first and $1500$ in
the second – among the same number of agents. Accordingly, we can observe the
resources allocated among all the agents are lower in the second case than the
first. Moving from the second example to the third (the bottom graph in Fig.
2), the variance is increased while keeping all other parameters the same. In
all these three cases, the free design parameter $v_{k}$ is $\mu$ – the mean
of the Gaussian distribution. The solution of the free parameter obtained from
solving the $N+1$ equations from (10), is shown in the figures and is used to
initialize the Lloyd’s algorithm. The generators obtained from the Lloyd’s
algorithm and the generators from solving (10) are plotted together. We can
observe that the two solutions are very close to each other. Additionally,
both the solutions sum up to the resource to be allocated – $r$, with an
acceptable error.

Figure 2: Allocation of $r$ amount of resource among 50 agents in
$\Omega=[0,100]$ under Gaussian distribution for specified variances – $4$
(top and middle) and $8$ (bottom). The mean of the distribution $\mu$ is the
solution $v_{k}$ from (10).

Similarly, we present another set of simulations in Fig. 3. In the three cases
therein, $\Omega,N$ and $r$ are the same. The difference in the three cases is
the underlying distributions – Gamma distribution in the top figure,
Exponential in the middle, and Gaussian distribution in the bottom figure.
Like in Fig. 2, the solutions from the two approaches are close to each other
and also sum up to $r$.

Figure 3: Allocation of $r$ amount of resource among 50 agents in
$\Omega=[0,300]$ under three different distributions. Top: Gamma distribution
with the free parameter $v_{k}$ being $k$. Middle: Exponential distribution
with the free parameter $v_{k}$ being $\lambda$. Bottom: Gaussian distribution
with the free parameter $v_{k}$ being $\mu$.

Even though in this approach we obtain the solution of (5), the SNLE method is
centralized and its grows with $N$. However, it is worth pointing out that
once the Lloyd’s algorithm is initialized with a design parameter $v_{k}$, CVT
obtained using Lloyd’s algorithm is scalable to any $N$ because regardless of
the total number of agents $N$, each agent can have at most two neighbors.

Like most of the solutions to the resource allocation problem, in this Section
we considered a fixed amount of resource to be allocated among al the agents.
Using the developed static allocation method as the initialization step, in
the next Section we consider the problem of dynamic resource allocation
problem where the amount of resource to be allocated is time-varying and all
the agents are aware of the quantity.

##  IV Dynamic Resource Allocation

In the previous section, we solved the static resource allocation problem by
using the centralized system of nonlinear equations. However, extending the
same approach to varying amount of resource-to-be-allocated results in a
centralized approach. Therefore, in this Section we focus on developing a
decentralized approach to the dynamic resource allocation problem.

Employing the static resource allocation problem as the initialization step,
our solution approach to the dynamic resource allocation problem under Normal
distribution involves a one-step update that maintains the dynamic resource
allocation constraint while preserving the CVT. We employ the following Lemma
1 to obtain such one-step update in Theorem 1. Through the design process, we
assume that the amount of resource to be allocated among all the agents over
the considered time duration is known to all the agents.

Suppose $\rho(.)=\mathcal{N}(\mu,\sigma^{2})$ and the “free” parameter $v_{k}$
is $\mu$. Then we have:

###### Lemma 1

Suppose at time $k$, $\\{z_{i}(k)\\}_{i\in I_{N}}$ are the centroids of the
CVT in $\Omega\subset\mathbb{R}$ with density
$\rho(.)=\mathcal{N}(\mu(k),\sigma^{2})$. Then the following relationship
holds between the time-updated centroids:

$\displaystyle z_{i}(k+1)-z_{i}(k)$ $\displaystyle=z_{j}(k+1)-z_{j}(k)$
$\displaystyle=\mu(k+1)-\mu(k)=-\delta$ (11)

Proof: Let $\mu(k+1)=\mu(k)-\delta$. Because $\\{z_{i}(k)\\}_{i\in I_{N}}$ are
the centroids with normal distribution, we have by definition:

$\displaystyle
z_{i}(k)=\frac{\int_{V_{i}(k)}xe^{\frac{(x-\mu(k))^{2}}{2\sigma^{2}}}dx}{\int_{V_{i}(k)}e^{\frac{(x-\mu(k))^{2}}{2\sigma^{2}}}dx}$

Similarly, writing out the mass centroid for the next time instant $k+1$ using
$\mu(k+1)=\mu(k)-\delta$, we have:

$z_{i}(k+1)=\frac{\int_{V_{i}(k+1)}xe^{\frac{(x-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dx}{\int_{V_{i}(k)}e^{\frac{(x-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dx}$
(12)

Suppose $V_{i}(k)=[a,b]\subset\Omega$. Consider the change of variables
$y=x-\delta$. Then the mass centroids transform as:

$\displaystyle z_{i}(k)$
$\displaystyle=\frac{\int_{a}^{b}xe^{\frac{(x-\mu(k))^{2}}{2\sigma^{2}}}dx}{\int_{a}^{b}e^{\frac{(x-\mu(k))^{2}}{2\sigma^{2}}}dx}$
$\displaystyle=\frac{\int_{a-\delta}^{b-\delta}(y+\delta)e^{\frac{(y+\delta-\mu(k))^{2}}{2\sigma^{2}}}dy}{\int_{a-\delta}^{b-\delta}e^{\frac{(y+\delta-\mu(k))^{2}}{2\sigma^{2}}}dy}$
$\displaystyle=\frac{\int_{a-\delta}^{b-\delta}(y+\delta)e^{\frac{(y-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dy}{\int_{a-\delta}^{b-\delta}e^{\frac{(y-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dy}$
$\displaystyle=\frac{\int_{a-\delta}^{b-\delta}ye^{\frac{(y-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dy+\int_{a-\delta}^{b-\delta}\delta
e^{\frac{(y-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dy}{\int_{a-\delta}^{b-\delta}e^{\frac{(y-(\mu(k)-\delta))^{2}}{2\sigma^{2}}}dy}$
$\displaystyle=\frac{\int_{V_{i}(k+1)}ye^{\frac{(y-\mu(k+1))^{2}}{2\sigma^{2}}}dy}{\int_{V_{i}(k+1)}e^{\frac{(y-\mu(k+1))^{2}}{2\sigma^{2}}}dy}+\delta\frac{\int_{V_{i}(k+1)}e^{\frac{(y-\mu(k+1))^{2}}{2\sigma^{2}}}dy}{\int_{V_{i}(k+1)}e^{\frac{(y-\mu(k+1))^{2}}{2\sigma^{2}}}dy}$
$\displaystyle=z_{i}(k+1)+\delta$ $\displaystyle\implies$ $\displaystyle
z_{i}(k+1)-z_{i}(k)=-\delta$ (13)

Since (13) holds for all $i\in I_{N}$ and $\mu(k+1)=\mu(k)-\delta$, we have
(11) proved.

$\hfill\square$

###### Theorem 1

Suppose we are initialized with static resource allocation solution at
discrete-time $k$ and are at the following conditions: $\\{z_{i}(k)\\}_{i\in
I_{N}}\ s.t\ \sum_{i\in I_{N}}z_{i}(k)=r(k),\ \ \\{z_{i}(k)\\}_{i\in
I_{N}}\sim\mathcal{N}(\mu(k),\sigma^{2})$. Suppose the resource to be
allocated at the next time instant is $r(k+1)$. If agents update their
resources as

$z_{i}(k+1)=z_{i}(k)+\frac{1}{N}(r(k+1)-r(k))$ (14)

then the resulting solution satisfies the following:

  * 1

$\sum_{i\in I_{N}}z_{i}(k+1)=r(k+1)$

  * 2

$\\{z_{i}(k+1)\\}_{i\in I_{N}}\sim\mathcal{N}(\mu(k+1),\sigma^{2})$

where $\mu(k+1)$ is a solution of the $(N+1)$ SNLE (10).

Proof: Obtain the time-difference of the summation of the resources:

$\displaystyle\sum_{i\in I_{N}}z_{i}(k+1)-\sum_{i\in I_{N}}z_{i}(k)$
$\displaystyle=\sum_{i\in I_{N}}z_{i}(k+1)-z_{i}(k)$ $\displaystyle=-N\delta$
$\displaystyle=N(\mu(k+1)-\mu(k))$ $\displaystyle=r(k+1)-r(k)$ (15)

$\hfill\square$

Following Theorem 1 we obtain the CVT that satisfies the dynamic resource
allocation constraint for the desired Normal distribution in a decentralized
manner. While this fulfills our objective, it can be observed that the
approach is quite rigid. In practical applications where the agents have their
own set of dynamics and are trying to navigate around certain local objectives
as well, this approach can be restrictive. Therefore, to extend its
applicability we introduce flexibility in the design by allowing for (local)
negotiations between neighbors through what we call a “civility model”.

Before detailing the civility model, let us introduce some new notations. For
each agent $i\in I_{N}$, denote its desired resource amount at time $k$ that
meets its local objective as $u_{i}(k)$. For example, if the agent $i$ is
responsible for the control of a certain system modeled as a state-space, such
$u_{i}(k)$ could be the control input from a state-feedback controller, from
an LQR or from any such local controller. Since we are operating in 1-D
spaces, recall from Section II-C that the resource and communication graphs
are the same. Following the same notation therein, denote the communication
graph at time $k$ as $\mathcal{C}^{k}$, and the neighbors of agent $i$ at time
$k$ as $\mathcal{N}_{\mathcal{C}^{k}_{i}}$.

Initialization: All agents are aware of the total resources $r(k),\forall k\in
T$ and the initial communication network $\mathcal{C}^{k-1}$. Solve the static
allocation problem for the resource $r(k-1)$.

Following the initialization, the civility model for local negotiations is
developed as follows.

Civility model for local negotiations

For every agent $i\in I_{N}$, at every time $k\in T$, do:

  * 1

Compute the resource update $z_{i}(k)$ from (14). Compute $u_{i}(k)$ based on
the local requirements, possibly from the local controller.

  * 2

Compute the neighbor of interest as
$\hat{j}=\\{j\in\mathcal{N}_{\mathcal{C}^{k}_{i}}\cup i\text{ such that
}||u_{i}(k)-z_{\hat{j}}(k)||<||u_{i}(k)-z_{j}(k)||\\}$.

  * 3

Swap resources with the neighbor of interest $\hat{j}$ from the previous step,
if $\hat{j}$ indicates it has not already been taken. This results in
$z_{i}(k)=z_{\hat{j}}(k)$. If $\hat{j}$ has already negotiated with its other
neighbor and is hence taken, or if $\hat{j}=i$, then implement the resource
update $z_{i}(k)$ from Step 1.

It is worth noting that the communication network is dynamically updated in a
decentralized manner, and that such an update naturally follows from the
resource swap during the local negotiations. We call this approach the
civility model because if a neighbor asks to swap, the agent complies with it
regardless of its own local requirement. And hence, since all the agents
follow the same model, no agent is at a disadvantage in following such an
approach.

To demonstrate the clarity and effectiveness of the proposed method to
dynamically allocate resources in a decentralized manner, we consider the
application of demand-response in smart grids. Specifically, we consider a
group of Heating, Ventilation, and Air Conditioning (HVAC) units that have
their local objectives to maintain their indoor air temperatures according to
certain desired setpoints, but are also required to respond to certain demand
(power) curve by consuming the available power as a team of agents.

##  V Application to Demand Response

To demonstrate the developed method, we consider power allocation in a group
of building HVACs. In this application of demand-response, the agents are the
building HVACs. The resources to be allocated to all the agents are the powers
consumed by the HVACs to maintain the local indoor air temperatures. We adapt
the state-space model from [31] to simulate the indoor air temperatures for
each agent $i$ as:

$\displaystyle\dot{x}_{i}(t)=A_{i}x_{i}(t)+B_{i}u_{i}(t)+G_{i}w_{i}(t)$
$\displaystyle y_{i}(t)=C_{i}x_{i}(t)+D_{i}u_{i}(t)$ (16)

The input $u_{i}$ is the power consumption of the HVAC (agent $i$), the output
$y_{i}$ is the indoor air temperature, and $w_{i}$ is the vector of
disturbances – outdoor air temperature and solar radiation. The system
matrices for each agent are given by:

$A_{i}=\begin{bmatrix}\frac{-(K_{1}^{i}+K_{2}^{i}+K_{3}^{i}+K_{5}^{i})}{C_{1}^{i}}&\frac{(K_{1}^{i}+K_{2}^{i})}{C_{1}^{i}}&\frac{K_{5}^{i}}{C_{1}^{i}}\\\
\frac{K_{1}^{i}+K_{2}^{i}}{C_{2}^{i}}&\frac{-(K_{1}^{i}+K_{2}^{i})}{C_{2}^{i}}&0\\\
\frac{K_{1}^{i}}{C_{3}^{i}}&0&\frac{-(K_{4}^{i}+K_{5}^{i})}{C_{3}^{i}}\end{bmatrix}$

$\displaystyle B_{i}=\begin{bmatrix}\frac{1}{C_{1}^{i}}+\frac{1}{C_{2}^{i}}\\\
0\\\
0\end{bmatrix}G_{i}=\begin{bmatrix}\frac{K_{3}^{i}}{C_{1}^{i}}&\frac{1}{C_{1}^{i}}\\\
0&\frac{1}{C_{2}^{i}}\\\
\frac{K_{4}^{i}}{C_{3}^{i}}&0\end{bmatrix}C_{i}=\begin{bmatrix}1&0&0\end{bmatrix}$

with $D_{i}$ being a zero matrix. The system parameters, which are resistances
and capacitances in the thermal dynamics of the building model, for each agent
$i$ are obtained as realizations of the following normal distributions:

$\displaystyle K_{1}\sim\mathcal{N}(16.48,0.1)\hfill$ $\displaystyle
K_{5}\sim\mathcal{N}(23.04,0.1)$ $\displaystyle
K_{2}\sim\mathcal{N}(108.5,0.1)\hfill$ $\displaystyle
C_{1}\sim\mathcal{N}(9.36\times 10^{5},1)$ $\displaystyle
K_{3}\sim\mathcal{N}(5,0.1)\hfill$ $\displaystyle
C_{2}\sim\mathcal{N}(2.97\times 10^{6},1)$ $\displaystyle
K_{4}\sim\mathcal{N}(30.5,0.1)\hfill$ $\displaystyle
C_{3}\sim\mathcal{N}(6.695\times 10^{5},1)$

We implement the agent’s model by discretizing the state-space model (16) with
a sampling time of 10 minutes. In the HVAC model, the input $u_{i}$
corresponds to cooling when negative and to heating when positive. Regardless,
its absolute value is the power consumed, and therefore we use that for local
negotiations and let the individual agent decide whether to use the allocated
power for heating or cooling based on its local control. To maintain the
indoor air temperatures from a local control, we employ a state-feedback
controller for pole placement for every agent to determine its $u_{i}(k)$. We
consider the same disturbances for all the agents; the outdoor air temperature
and the solar radiation, [32], we use for our simulations are shown in Fig. 4.

Figure 4: Disturbances in the HVAC model (16)

To begin the dynamic resource allocation we initialize $\rho(.)$ as
$\mathcal{N}(\mu,\sigma^{2})$, and following Section III, solve the first-time
allocation (initialization) as a static allocation problem. Communicating to
all the agents the resulting mean $\mu$, we begin the decentralized dynamic
allocation as laid out in Section IV.

Even though the performance of the developed approach depends on the total
available resource and the local requirements, the civility model allows for
flexibility, and that could be necessary to compensate for local disturbances
or for improper selection of the (desired) distribution in the tessellation.
To explain the graphical setup of our results, we begin with $N=5$ in Fig 5.
The top figure shows the power consumption of all the agents at every time
instant, and the bottom figure shows their total power consumption versus the
available power. Augmenting, Fig 6 shows the individual indoor air
temperatures when the agents implement the allocated power from Fig 5.

Figure 5: Baseline power consumptions: Agents acting based on the resource
allocation constraint without the civility model. Left: Individual power
consumption. Right: Total power consumption. Figure 6: Baseline indoor air
temperatures: Agents acting based on the resource allocation constraint
without the civility model.

Next, we demonstrate the civility model from Section IV by allowing for
swapping through local negotiations. Continuing the previous case we first
consider only $5$ agents in the team in Fig. 7 and then demonstrate for $15$
agents. For every agent, the power consumptions and the indoor air
temperatures are shown in the same color throughout the simulation duration.
For example, agent $2$ is shown in red. Thus one can follow the agents’
negotiations and the resulting swaps and communication network by following
the individual power consumption of the agents through their colors. In the
subsequent cases, we do not show the satisfaction of the resource allocation
constraint through a dedicated figure since we can concisely express it
numerically as the error between total power consumption of all the agents and
the available power; we use $l_{2}$ norm to compute the power consumption
error.

Figure 7: Civility model with local state feedback controller for $5$ agents.

Figure 8: Civility model with local state feedback controller for $15$ agents.
The temperature setpoints for all the HVACs are at $72\degree F$.

The strengths of the developed method lie in its robustness in maintaining the
resource allocation constraint while accounting for local preferences in a
truly decentralized manner. To demonstrate the same, we perturb the setpoints
of certain agents and observe the corresponding resource negotiations and the
air temperatures in Fig 9. We can observe the increased amount of negotiations
in the increased number of swaps spreading throughout the team to correct for
the disturbances for some of the agents. Quantifying the swaps, we have that
out of $144$ time-steps in the simulation, each agent swapped $126.9$ times on
average and that every agent has been a neighbor of almost every other agent.
This suggests a high degree of variation in the communication network, further
suggesting that the amount of information is so fragmented among all the
agents that it is sufficient to meet the resource allocation constraint while
following the desired distribution in the tessellation but not enough for any
agent to recreate the behavior of any other agent.

Figure 9: With swaps and local state feedback controller of 15 agents under
disturbed setpoints

The decentralized dynamic resource allocation solution proposed in this work
follows the idea of “Global trendsetting, local negotiations”. Here, the the
global trend is for the agents’ resources to be Gaussian distributed while
summing up to the available power, and the local negotiations happen to
maintain the balance between following such global trend and accounting for
the local requirements simultaneously.

##  VI Conclusions

CVTs in one-dimensional spaces are desirable due to their inherent line
structure and ease of computation of the entire tessellations, allowing for
verification of the quality of the solution. Employing them in the resource
allocation problem brings forth additional advantages as embedding desired
global trends in the team through probability distributions.

For a fixed amount of resource, the static resource allocation method offers
an analytical, although central, solution by posing the constrained CVT
minimization problem as a system of non-linear equations. The generalizability
of the developed static allocation framework is worth remarking. Instead of
the summation constraint, one can have any constraint from
$\mathbb{R}^{N}\to\mathbb{R}$, and one can also have as many constraints as
the number of parameters defining the desired distribution.

The developed decentralized dynamic allocation solution using CVTs provides a
natural way to embed the aggregate team behavior or to set the desired global
trend through the distribution of the tessellations. The developed method
through the civility model allows for flexibility on the local end by
absorbing and distributing disturbances throughout the team. We observed a
hint of inherent privacy in the architecture through the highly dynamic
communication network that nevertheless was a simple line graph at all times,
demonstrating the scalability of the method.

Building on this work, we aim to generalize the develop decentralized resource
allocation method to global trends that are described by more distributions,
and not just Gaussian. We aim to verify the robustness of the architecture to
a changing number of agents in the team, possibly due to communication
failures.

## References

  * [1] K. Arrow and F. H. Hahn, General Competitive Analysis.  North Holland, December 1983.
  * [2] H. Lakshmanan and D. P. de Farias, “Decentralized resource allocation in dynamic networks of agents,” SIAM Journal on Optimization, vol. 19, no. 2, pp. 911–940, 2008.
  * [3] L. Xiao and S. Boyd, “Optimal scaling of a gradient method for distributed resource allocation,” Journal of Optimization Theory and Applications, vol. 129, no. 3, pp. 469–488, 2006.
  * [4] G. Banjac, F. Rey, P. Goulart, and J. Lygeros, “Decentralized resource allocation via dual consensus ADMM,” in 2019 American Control Conference (ACC), pp. 2789–2794, 2019.
  * [5] T. T. Doan and C. L. Beck, “Distributed resource allocation over dynamic networks with uncertainty,” IEEE Transactions on Automatic Control, vol. 66, no. 9, pp. 4378–4384, 2021.
  * [6] H. Cai, “Decentralized control of stochastic dynamic systems with applications to resource allocation and portfolio management,” 2012.
  * [7] J. Snow, On the Mode of Communication of Cholera.  1855\.
  * [8] L. Ju, T. Ringler, and M. Gunzburger, Voronoi Tessellations and Their Application to Climate and Global Modeling, pp. 313–342.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2011.
  * [9] J. Cortes, S. Martinez, T. Karatas, and F. Bullo, “Coverage control for mobile sensing networks,” IEEE Transactions on Robotics and Automation, vol. 20, no. 2, pp. 243–255, 2004.
  * [10] D. Inoue, Y. Ito, and H. Yoshida, “Optimal transport-based coverage control for swarm robot systems: Generalization of the Voronoi tessellation-based method,” IEEE Control Systems Letters, vol. 5, no. 4, pp. 1483–1488, 2021\.
  * [11] Y. Bai, Y. Wang, M. Svinin, E. Magid, and R. Sun, “Adaptive multi-agent coverage control with obstacle avoidance,” IEEE Control Systems Letters, vol. 6, pp. 944–949, 2022.
  * [12] A. Davydov and Y. Diaz-Mercado, “Sparsity structure and optimality of multi-robot coverage control,” IEEE Control Systems Letters, vol. 4, no. 1, pp. 13–18, 2020.
  * [13] T. Chen, H. Pourbabak, and W. Su, “Electricity market reform,” The Energy Internet, pp. 97–121, 2019.
  * [14] S. Williams and M. Short, “Electricity demand forecasting for decentralised energy management,” Energy and Built Environment, vol. 1, no. 2, pp. 178 – 186, 2020.
  * [15] C. Park and T. Yong, “Comparative review and discussion on P2P electricity trading,” Energy Procedia, vol. 128, pp. 3–9, 2017.
  * [16] C. Zhang, J. Wu, C. Long, and M. Cheng, “Review of existing peer-to-peer energy trading projects,” Energy Procedia, vol. 105, pp. 2563–2568, 2017\.
  * [17] A. Malizou, “Electricity aggregators: Starting off on the right foot with consumers,” BEUC, The European Consumer Organization, 2018.
  * [18] S. Burger, J. Chaves-Ávila, C. Batlle, and I. Pérez-Arriaga, “A review of the value of aggregators in electricity systems,” Renewable and Sustainable Energy Reviews, vol. 77, pp. 395–405, 2017.
  * [19] Q. Du, V. Faber, and M. Gunzburger, “Centroidal voronoi tessellations: Applications and algorithms,” SIAM Review, vol. 41, no. 4, pp. 637–676, 1999.
  * [20] P. Fleischer, “Sufficient conditions for achieving minimum distortion in a quantizer,” IEEE International Convention Record, Pt I, pp. 104–111, 1964\.
  * [21] J. C. Urschel, “On the characterization and uniqueness of centroidal Voronoi tessellations,” SIAM Journal on Numerical Analysis, vol. 55, no. 3, pp. 1525–1547, 2017.
  * [22] S. Lloyd, “Least squares quantization in pcm,” IEEE Transactions on Information Theory, vol. 28, no. 2, pp. 129–137, 1982.
  * [23] Q. Du, M. Emelianenko, and L. Ju, “Convergence of the lloyd algorithm for computing centroidal voronoi tessellations,” SIAM Journal on Numerical Analysis, vol. 44, no. 1, pp. 102–119, 2006.
  * [24] J. Kieffer, “Uniqueness of locally optimal quantizer for log-concave density and convex error weighting function,” IEEE Transactions on Information Theory, vol. 29, no. 1, pp. 42–47, 1983.
  * [25] Y. Liu, W. Wang, B. Lévy, F. Sun, D.-M. Yan, L. Lu, and C. Yang, “On centroidal Voronoi tessellation—energy smoothness and fast computation,” ACM Trans. Graph., vol. 28, sep 2009.
  * [26] X. Wang, X. Ying, Y.-J. Liu, S.-Q. Xin, W. Wang, X. Gu, W. Mueller-Wittig, and Y. He, “Intrinsic computation of centroidal voronoi tessellation (cvt) on meshes,” Computer-Aided Design, vol. 58, pp. 51–61, 2015.  Solid and Physical Modeling 2014.
  * [27] J. Hatless, H. Wei, and C. L, “Fast methods for computing centroidal voronoi tessellations,” Journal of Scientific Computing, vol. 63, pp. 185–212, 2015.
  * [28] J. Macqueen, “Some methods for classification and analysis of multivariate observations,” In 5-th Berkeley Symposium on Mathematical Statistics and Probability, pp. 281–297, 1967.
  * [29] M. Loebl, Introduction to Graph Theory, pp. 13–49.  Wiesbaden: Vieweg+Teubner, 2010.
  * [30] A.-M. Wazwaz, Fredholm Integral Equations, pp. 119–173.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2011.
  * [31] X. Ma, J. Dong, S. M. Djouadi, J. J. Nutaro, and T. Kuruganti, “Stochastic control of energy efficient buildings: A semidefinite programming approach,” in 2015 IEEE International Conference on Smart Grid Communications (SmartGridComm), pp. 780–785, 2015.
  * [32] K. Amasyali, Y. Chen, B. Telsang, M. Olama, and S. M. Djouadi, “Hierarchical model-free transactional control of building loads to support grid services,” IEEE Access, vol. 8, pp. 219367–219377, 2020.


# Exploring Perceived Vulnerability of Pedestrians: Insights from a Forced-
Choice Experiment

Paul Geoerg
Vereinigung zur Förderung des Deutschen Brandschutzes e.V.
Münster
Germany
<EMAIL_ADDRESS>
Ann Katrin Boomers
Forschungszentrum Jülich GmbH
Jülich
Germany
<EMAIL_ADDRESS>
Maxine Berthiaume
University of Ottawa
Ottawa
Canada
<EMAIL_ADDRESS>
Maik Boltes
Forschungszentrum Jülich GmbH
Jülich
Germany
<EMAIL_ADDRESS>
Max Kinateder
National Research Council Canada
Ottawa
Canada
<EMAIL_ADDRESS>

###### Abstract

Individual differences in mobility (e.g., due to wheelchair use) during crowd
movement are not well understood. Perceived vulnerability of neighbors in a
crowd could affect, for example, how much space is given to them by others. To
explore how pedestrians perceive people moving in front of them, in
particular, how vulnerable they believe them to be, we asked $51\text{\,}$
participants to complete a Two-Alternatives-Forced Choice task (2AFC) in an
internet browser. Participants were shown pairs of images each showing a
person and then asked to select the person who appeared more vulnerable to
them. For example, participants would choose between a male person in a
wheelchair and a female person carrying a suitcase. In total $16\text{\,}$
different stimuli (male vs female; no item/device, 1 suitcase, 2 suitcases,
small backpack, large backpack, stroller, cane, and wheelchair), yielding
$n(n-1)/2=120$ potential pairwise comparisons per participant. Results showed
that wheelchair users appeared the most vulnerable and persons without any
items/devices the least vulnerable. Persons carrying two suitcases were in the
middle. These results informed the design of a main behavioral study (not
reported here).

_Keywords_ Pedestrian dynamics $\cdot$ Accessibility $\cdot$ Heterogeneous
crowds $\cdot$ Demographic change $\cdot$ 2AFC $\cdot$ Online study $\cdot$
Perceived vulnerability $\cdot$ Wheelchair user

##  1 Introduction

This manuscript documents the background, rationale, procedure, and results of
a pilot study designed to test the perceived vulnerability of pedestrians
walking ahead of another person. The results of the pilot study informed the
design of a larger behavioral study on pedestrian movement (not reported here,
see [1])

Individual differences in mobility (e.g., due to wheelchair use) during crowd
movement are not well understood. This has been recognized in several review
publications, e.g. [2, 3, 4]. While these have pointed out how, for example,
certain functional limitations could affect pedestrian egress movement, they
also report a lack of empirical data. Likely as a consequence, many
engineering tools that aim to predict the evacuation performance of crowds are
based exclusively on data from young adults without disabilities. However,
previous work has found anecdotal evidence that pedestrians keep a larger
distance from wheelchair users when walking in a crowd [5].

The perceived vulnerability of people walking next to each other (also called
neighbors) in a crowd could potentially explain this effect on microscopic
movement parameters (e.g., movement speed). More specifically, a relevant
question that has not been answered is how the visibility of a disability
(e.g., through recognition of an assistive device) or mobility-relevant
properties (e.g., carrying heavy items, pushing strollers, etc.) shape how
people see neighbors in a crowd.

Individual pedestrians’ behavioral reactions to visible mobility attributes
may have cascading effects on micro- and macroscopic movement patterns in the
crowd. Here, we explore two potential explanatory mechanisms: perceived
vulnerability and perceived required space. That is, will people increase
their interpersonal distance to wheelchair users because there is a social
norm to be mindful of people who appear vulnerable or because wheelchair users
simply appear to take up more space?

The purpose of this pilot study was to collect subjective impressions from
participants via an online study that informed the main behavioral study (a
group of participants moving through a bottleneck). Participants reported
ratings of perceived vulnerability of pedestrians and wheelchair users with
different mobility attributes using a Two-Alternatives-Forced-Choice (2AFC)
paradigm [6].

In a 2AFC paradigm, participants were shown two stimuli (e.g., images) and
then need to select one. The 2AFC procedures allowed for generating unbiased
responses and avoid certain response patterns compared to other procedures
such as simple rating scales (e.g., anchoring) [6, 7, 8].

###  1.1 Research Questions

The main research questions addressed in this pilot study were:

  1. 1.

Are wheelchair users generally perceived to be more or less vulnerable than
others?

  2. 2.

Are there differences/similarities compared to other mobility-related
attributes?

##  2 Methods

###  2.1 Design

In each trial, two images showing a person were presented. The participant’s
task was to select the person that appeared more vulnerable to them. More
specifically, they were given the following instructions: "Imagine that you
are walking behind each person while walking in a crowd. Please click on the
person you would be more cautious around" (Figure 1(c)).111We queried feedback
on optimal item formulation from an expert in ecological psychology and
experimental design Participants then chose between images in which the
following attributes were manipulated:

  1. 1.

Mobility attributes (8 levels): baseline (no mobility aid or luggage), 1
suitcase, 2 suitcases, small backpack, large backpack, stroller, cane,
wheelchair

  2. 2.

Visible gender attributes of the person shown (male or female)

This yielded a total of 16 unique stimuli (see 2.2). Stimuli were placed
randomly either on the left or right side of the screen (see Figure 1(c)).
Conditions in which two identical images were shown on the left and right side
were excluded. In a 2AFC design, this translates into $n(n-1)/2=120$ potential
pairwise comparisons.

###  2.2 Stimuli

Sixteen stimuli were generated for the main trials of this study (Figure 1(b))
as well as three stimuli for the practice trials (Figure 1(a)). Each stimulus
was generated from a vector graphic based on two baseline models purchased
from a commercial stock image database (shutterstock.com) and then loaded into
an image editing software (Gimp v.2.10.32, www.gimp.org). The baseline images
showed stylized views of an adult person, with either stereotypical male or
female visual attributes. The "male" character had short hair, a jacket and
pants. The "female" character had long hair and was wearing a skirt. Both
characters were Caucasian. The baseline models were then modified for varying
mobility profiles. Care was taken so that the basic appearance (e.g., posture,
size) of the baseline characters did not change, except for the changes
associated with the mobility profiles. The three practice trial stimuli showed
similar characteristics, however, their design differed so that no feature
(e.g., assistive device) would match those of the main trials. They included a
"female" character with crutches, a female character with grey hair, and an
adult male character (Figure 1(a)).

(a)

(b)

(c)

Figure 1: Pilot study: (a) stimuli practice trials; (b) stimuli main trials;
and (c) example screenshot of 2AFC task.

###  2.3 Procedure

The study was implemented using PsychoPy (v2021.1.4) and hosted online using
Pavlovia.org (an online platform to run, share, and explore psychometric
studies). Participants needed to use an internet-connected device (either
computer, tablet, or phone) to access the study. The study procedure itself is
structured in a sequential order (Fig 2. Participants had to give informed
consent before starting the study procedure, were instructed, practiced three
exemplary trials to get familiar with the study design, and then started the
120 decision trials. Last, participants were asked for demographic information
(age group, gender identity, prior experience with people with disabilities,
and movement in crowds). The procedure was approved by the NRC Research Ethics
Board (REB 2021-88).

Figure 2: Pilot study structure and procedure.

###  2.4 Sample

Participants were recruited via e-mail invitation and posting on social media.
$51\text{\,}$ participants completed at least some of the main trials of the
study. $64.7\text{\,}\mathrm{\char 37\relax}$ ($n=$33\text{\,}$$) of the
participants provided a complete data set and provided demographic
information. Two participants ($\approx$6.1\text{\,}\mathrm{\char 37\relax}$$)
stated that they had a disability. Figure 3 shows the distribution of age,
gender, and the amount of experience that participants had interacting with
people with disabilities in their daily lives.

(a)  (b)  (c)

Figure 3: Distribution of (a) age, (b) gender, and (c) reported frequency of
interaction with persons with disabilities.

##  3 Results and Discussion

In this study, we investigated how visible mobility attributes, such as using
a wheelchair, carrying luggage, or pushing a stroller affect the perceived
vulnerability of pedestrians. 2AFC designs provide two metrics: Which image
was selected (choice) and how long it took to make that decision. The former
can provide direct answers to the research question (see section 1.1) and test
for potential biases and stereotypical response patterns (e.g., a participant
always clicking on the right stimulus). The latter can provide information on
how difficult a decision was (the longer participants needed to decide between
two stimuli, the harder the decision). In addition, response times can be used
to filter non-credible responses (either too fast or too long, for a
discussion in an adjacent area see [9]).

###  3.1 Data processing

In order to exclude data with non-credible responses, we first filtered the
data for response time by removing values above and below two standard
deviations of the average data (response time
$\bar{RT}=$1.67\text{\,}\mathrm{s}$\pm$7.62\text{\,}\mathrm{s}$$). Note that
this approach has limitations (see [10] for a discussion). As a result,
$4.45\text{\,}\mathrm{\char 37\relax}$ ($265\text{\,}$) from a total of
$5949\text{\,}$ responses were removed from the original data set.

Next, we plotted histograms of the left and right clicks of participants to
identify potential biases (see Figure 6(a)). None of the participants appeared
to systematically prefer left to right and consequently, thus no participants
were excluded.

###  3.2 Choice

Figure 4 shows the absolute frequencies that each stimulus was selected. The
following observations were made: Stimuli showing …

  1. 1.

…a person without any assistive device were consistently rated as the least
vulnerable;

  2. 2.

…a person in a wheelchair appeared to be the most vulnerable, followed by
those with a cane and stroller;

  3. 3.

…persons with either one or two suitcases appeared to be in the middle;

  4. 4.

…persons carrying larger and smaller backpacks appeared less vulnerable than
those carrying suitcases but more than a person without any assistive device
or travel item;

  5. 5.

…persons with stereotypically female attributes consistently appeared to be
more vulnerable than those with male attributes. These differences were most
pronounced for the stimuli showing cane users, and smallest for stimuli
showing wheelchair and stroller users.

This indicates that the perception of vulnerability scales with the space
required by the item, i.e., a small backpack is perceived as less vulnerable
than a large backpack, which is perceived as less vulnerable than a stroller.
Similarly, the larger an item, the more vulnerable a person appears to be.
Interestingly, this effect does not seem to be linear. For example, the
differences in perception from small to big backpacks and from one to two
suitcases were larger compared to the differences between other conditions.

Figure 4: Absolute frequencies of each stimulus being selected.

###  3.3 Response times

Figure 6(b) shows the response times for each stimulus, when selected. The
average response time across all stimuli was
$$1.40\text{\,}$\pm$1.79\text{\,}\mathrm{s}$$ (median
$=$0.81\text{\,}\mathrm{s}$$). Medians were consistently higher than averages.
There were no noticeable differences in response times across stimuli
(regardless of displayed gender attributes and items). Note that the number of
data points varied notably because we only report the response times for the
selected stimulus (see Figure 4).

Figure 5 shows a heatmap of the average response time for each pairwise
comparison. These data could indicate comparisons that were harder for
participants to differentiate. The longest response times were reported for
comparisons between male and female cane users ($2.83\text{\,}\mathrm{s}$),
male with one suitcase and female with a large backpack
($2.79\text{\,}\mathrm{s}$), and male and female wheelchair users
($2.74\text{\,}\mathrm{s}$). The fastest response times were reported for
comparisons between females with no items and female wheelchair users
($0.47\text{\,}\mathrm{s}$). This suggests the following patterns:

  * •

Decisions between male and female figures were harder when they were otherwise
similar in terms of mobility attributes.

  * •

Decisions became easier when the contrast between items was larger (see, for
example, the column for the male wheelchair user ($m_{whee}$) in Figure 5).

Figure 5: Heatmap of average response times as a function of pairwise
comparisons.

(a)  (b)

Figure 6: (a) Histograms of left and right clicks for each participant ID, and
(b) boxplots of response times for each stimulus, when selected.

##  4 Conclusions

This study reports on a pilot experiment designed to provide insights and
guidance for a larger behavioral study [1]. The goal was to identify mobility
profiles (i.e., attributes of people moving in a crowd) that differ in
perceived vulnerability and space requirements. To this end, we asked
participants to compare images showing characters that differed in their
appearance, such as their stereotypical gender attributes (male vs female) or
the kind of travel items (suitcases and backpacks) and assistive devices (cane
and wheelchair) with them.

We chose a 2AFC online study for this purpose and while this approach has
clear limitations, a pattern could clearly be established, which answered the
research questions and informed the design of the behavioral study:

  * •

Pedestrians in wheelchairs consistently appeared to be more vulnerable than
any other mobility profile (research question 1). The inverse was true for
pedestrians without any other item/mobility device (research question 2).

  * •

Other mobility-related attributes appeared to influence the perceived
vulnerability as a function of their size, weight, and ease with which the
items could be moved (e.g., characters with bag packs, in general, were seen
as less vulnerable as those with suitcases and strollers). The only exception
was the cane, which was likely interpreted as a clear sign of vulnerability,
but to a lesser degree than the wheelchair (research question 2).

  * •

Within each category of items, characters that appeared to be female were
consistently rated to be more vulnerable.

  * •

The more similar the displayed stimuli were, the longer participants needed to
decide.

The present work has a number of limitations that should be considered.
Firstly, we only reported on a selection of mobility profiles. The profiles
were selected with the future behavioral study in mind, but the list is
certainly not exhaustive. For instance, we did not vary the appearance of age
(neither children nor seniors were displayed). Secondly, we did not measure
perceived vulnerability directly, but asked participants to indicate the
person they would be "more cautious" around. This approach was chosen given
that the term ’vulnerable’ may be interpreted differently based on the
participants’ own backgrounds. However, prior to testing, we solicited
feedback on this formulation from native speakers of English who were experts
in environmental psychology and experimental design. Finally, a large number
of participants did not complete the full study; Consequently, the data set
might be imbalanced given the amount of missing demographic information. In
addition, the influence of the participant’s own mobility profile (e.g.,
gender, age, or living with a disability) could not be investigated.

However, we believe that the limitations did not prevent this work from
achieving its goal of informing the design of a behavioral study. In the
behavioral study, groups of participants will be asked to move together
through a bottleneck; critically, the mobility attributes of two participants
at the center of the group will be manipulated. In wheelchair conditions, two
participants will be wheelchair users, in the luggage condition, two
participants will be carrying two suitcases each, while in control conditions
all participants will have similar mobility attributes.

## References

  * [1] Paul Geoerg et al.  “Together apart: the influence of increased crowd heterogeneity on crowd dynamics at bottlenecks”  In _Proceedings of Pedestrian and Evacuation Dynamics 2023_ , 2023
  * [2] Olivera Bukvic et al.  “A review on the role of functional limitations on evacuation performance using the International Classification of Functioning, Disability and Health”  In _Fire Technology_ 57.2, 2021, pp. 507–528  DOI: 10.1007/s10694-020-01034-5
  * [3] G. Carlsson et al.  “A scoping review of public building accessibility”  In _Disability and health journal_ 15.2, 2022, pp. 101227  DOI: 10.1016/j.dhjo.2021.101227
  * [4] Haley Hostetter and M.. Naser  “Characterizing disability in fire evacuation: A progressive review”  In _Journal of Building Engineering_ 53, 2022, pp. 104573  DOI: 10.1016/j.jobe.2022.104573
  * [5] Paul Geoerg, Jette Schumann, Maik Boltes and Max Kinateder  “How people with disabilities influence crowd dynamics of pedestrian movement through bottlenecks”  In _Scientific Reports_ 12.1, 2022  DOI: 10.1038/s41598-022-18142-7
  * [6] Gustav Theodor Fechner  “Elements of psychophysics, 1860”  In _Readings in the history of psychology_ East Norwalk: Appleton-Century-Crofts, 1948, pp. 206–213  DOI: 10.1037/11304-026
  * [7] Rafal Bogacz et al.  “The physics of optimal decision making: a formal analysis of models of performance in two-alternative forced-choice tasks”  In _Psychological review_ 113.4, 2006, pp. 700–765  DOI: 10.1037/0033-295X.113.4.700
  * [8] Carl Auerbach  “Correcting two-alternative forced-choice data for response bias”  In _Perceptual and motor skills_ 32.2, 1971, pp. 533–534  DOI: 10.2466/pms.1971.32.2.533
  * [9] Yoram Braw  “Response Time Measures as Supplementary Validity Indicators in Forced-Choice Recognition Memory Performance Validity Tests: A Systematic Review”  In _Neuropsychology review_ 32.1, 2022, pp. 71–98  DOI: 10.1007/s11065-021-09499-z
  * [10] ChristianMichael Lachaud and Oliver Renaud  “A tutorial for analyzing human reaction times: How to filter data, manage missing values, and choose a statistical model”  In _Applied Psycholinguistics_ 32.2, 2011, pp. 389–416  DOI: 10.1017/S0142716410000457