# Efficient Large-Scale Simulation of Fish Schooling Behavior Using Voronoi
Tessellations and Fuzzy Clustering

Salah Alrabeei
Department of Computer Science
Western Norway Univ. of Applied Sciences
Bergen, Norway
<EMAIL_ADDRESS>
& Sam Subbey
Inst. of Marine Research
PO Box 1870
5817 Bergen, Norway
<EMAIL_ADDRESS>
&Talal Rahman
Department of Computer Science
Western Norway Univ. of Applied Sciences
Bergen, Norway
<EMAIL_ADDRESS>
Also at Department of Computer Science, Western Norway Univ. of Applied
Sciences, PO Box 7030, 5020 Bergen, Norway

###### Abstract

This paper introduces an efficient approach to reduce the computational cost
of simulating collective behaviors, such as fish schooling, using Individual-
Based Models (IBMs). The proposed technique employs adaptive and dynamic load-
balancing domain partitioning, which utilizes unsupervised machine-learning
models to cluster a large number of simulated individuals into sub-schools
based on their spatial-temporal locations. It also utilizes Voronoi
tessellations to construct non-overlapping simulation subdomains. This
approach minimizes agent-to-agent communication and balances the load both
spatially and temporally, ultimately resulting in reduced computational
complexity.

Experimental simulations demonstrate that this partitioning approach
outperforms the standard regular grid-based domain decomposition, achieving a
reduction in computational cost while maintaining spatial and temporal load
balance. The approach presented in this paper has the potential to be applied
to other collective behavior simulations requiring large-scale simulations
with a substantial number of individuals.

##  1 Introduction

Collective behavior is a natural phenomenon observed in various dynamic
complex living systems, where large groups of autonomous individuals
participate based on limited local information transferred between each other
through local interaction or their reaction to the surrounding environment.
Such individual behaviors spread through the system, producing a collective
pattern. Collective behaviors can be observed in biological and ecological
systems such as schooling in fish [1], flocking in birds [2], herding mammals
[3], and human crowding [4]. They are also observed in microorganisms such as
cell populations [5], physics [6], and robotics [7].

Agent-based modeling (ABM) is a powerful tool used to simulate and capture
such emergent collective behaviors in complex systems [8]. Continuum-based
models, in which system variables evolve according to continuous functions
over time, have limitations in modeling and simulating such collective
behaviors and capturing emerging patterns [9]. Simulating realistic collective
behaviors often requires a large number of agents. However, running large-
scale agent-based simulations has been a major challenge in ABM. This
challenge is being addressed through the use of high-performance parallel
computing and domain decomposition algorithms, which can be efficiently
utilized to run such simulations. In particular, domain decomposition
techniques have been commonly used by the research community in this field.
Large-scale simulation workloads can be implemented by allocating each
computing resource (processor, cluster, etc.) to a part of the simulation
domain.

In domain decomposition-based simulation, selecting the proper partitioning
mechanism for the simulated agent system is the main challenge. Two main
approaches are used in the literature: grid-based and cluster-based [10]. The
grid-based approach consists of dividing the simulation domain into non-
overlapping and uniform subdomains, each with a set of individuals residing
temporarily in that subdomain. The cluster-based approach groups simulated
individuals based on specific criteria and solves each subdomain and its
simulated individuals or cluster separately in various modes, including
sequential, parallel, or disturbed versions [11]. Both approaches can be
either static or dynamic, but the static method fails to maintain load
balancing throughout the simulation evolution [9, 12].

Different types of partitioning have been implemented in grid-based domain
decomposition, with most techniques tailored to the specific needs of a
particular application area [13, 14, 15]. The simplest type of domain
decomposition is to divide the simulation domain into equal rectangular
subdomains, but other techniques have been developed, such as regular micro-
cells grouped into uniform subdomains, striped decomposition, block-striped
decomposition, or even irregular grid cells [16, 17, 18, 19, 20]. However,
using a partitioning method without an explicit load-balancing process can
lead to increased computational time due to the extensive information exchange
between agents at subdomain boundaries [20].

K-means is a widely used unsupervised and non-deterministic learning algorithm
for clustering large datasets. It has been particularly useful in low-
dimensional data and as a partitioning technique in agent-based simulations
[21, 22]. For instance, it has been adapted to partition agent-based crowd and
vortex particle simulations [23, 24]. Several improved versions of k-means-
based partitioning techniques have also been developed to efficiently
distribute simulation workloads across computing resources [25, 26, 27].

However, the k-means algorithm becomes less efficient when it needs to account
for individual movement in the simulation model [12]. Furthermore, it does not
consider cluster cardinalities, resulting in uneven cluster sizes when the
dataset has varying densities, such as in particle interactions where
aggregation and repulsion are present [28]. Consequently, the simulation
workload may become unbalanced, leading to expensive and inefficient large-
scale simulators.

This study introduces a hybrid cluster- and region-based partitioning approach
in a fish schooling simulation model. Our method is based on fuzzy (soft)
clustering and Voronoi tessellations, ensuring spatial and temporal load
balancing. The algorithm can be implemented sequentially or in parallel and
can run on a local or remotely distributed computing resource. In Section 2,
we introduce the biological fish schooling behavior model. In Section 3, we
describe our simulation domain partitioning and population cluster algorithm.
We present simulation results and a comparison between the static rectangular
decomposition grid-based partitioning method and our proposed algorithm in
Section 4. Finally, we conclude by discussing future work in Section 5.

##  2 Description of fish schooling model

Collective behaviors, such as fish schooling, bird flocking, and mammal
crowding, have been widely studied using an agent-based modeling approach.
This remains an active area of research in computational simulations. One of
the earliest models for studying collective behavior is the Boid (or Reynolds)
model [29]. In this model, each individual aligns its movement direction with
its neighbors, moves away from those neighbors who are too close to avoid
collisions, and then flocks together. Subsequent extensions of the Boid model
have incorporated additional features, such as obstacle avoidance [30, 31],
predator avoidance [32, 33], leadership [34, 35], field preference attraction
[36, 37], and multiple interaction zones [38]. Other variations include the
type of neighboring measure used [39].

Here, we adopt the Boid model [29] to capture the collective behaviors of
fish. Within a local range, each fish follows several interacting rules as
illustrated in Figure 1. The cohesion rule is composed of two parts: (i)
collision avoidance, which generates a repulsive force to ensure a minimum
distance between schoolmates (see Figure 2(a)), and (ii) school attraction,
which produces a positive force to maintain school formation within the region
of interaction (see Figure 2(c)). By combining these two rules, the fish
school within a mutual range, achieving a cohesive group behavior. The
cohesion force vector is denoted by $F_{u_{c}}$, and its mathematical
expression is given in Eq.(1).

$F_{u_{c}}=\frac{u_{c}}{\left\|u_{c}\right\|}\cdot\begin{cases}\frac{\left\|u_{c}\right\|-\text{r}_{\mathrm{keep}}}{\text{r}_{\mathrm{keep}}}&\text{if}\
\left\|u_{c}\right\|\leq\text{r}_{\mathrm{interact}},\\\
0&\text{Elsewhere}.\\\ \end{cases}$ (1)

where $u_{c}$ is the vector from the reference individual to the nearest
schoolmate, $\text{r}_{\mathrm{keep}}$ is the minimum safe distance to avoid
collision with a schoolmate, and $\text{r}_{\mathrm{interact}}$ is the maximum
distance to interact (range of the region of interaction).

Figure 1: Regions of interaction, $r_{interact}$, and repulsion, $r_{keep}$,
are centered on reference fish, $rf$ (black dot).

(a) Direction change (black to green arrow) to avoid too close schoolmates
(red dots).

(b) Movement towards the centroid of positions where schoolmates are located.

(c) Direction change (black to blue arrow) to align movement with schoolmates.

Figure 2: Fish local interaction rules – the reference individual fish $rf$
(black dot) whose original direction (black arrow) is influenced by
neighboring fish (red and blue dots) whose movement directions are indicated
by red and blue arrows.

The alignment rule in the fish school model is described by the force vector
$F_{u_{a}}$ which represents the tendency of an individual fish to align its
swimming velocity with its nearby schoolmates. The mathematical formula for
$F_{u_{a}}$ is given by Eq.(2):

$F_{u_{a}}=\frac{1}{n}\sum_{i\in R_{n}}u_{i}-u_{rf},$ (2)

where $u_{i}$ is the velocity vector of individual $i$, $R_{n}$ is the set of
$n$ individuals located within the interaction region of the reference
individual $rf$, and $u_{rf}$ is the velocity vector of the reference
individual.

In addition to the cohesion and alignment rules, there are also attraction and
repulsion rules that are based on stimuli and risks in the environment,
respectively [36]. The force vectors $F_{u_{f}}$ and $F_{u_{r}}$ for
attraction and repulsion, respectively, are given by Eq.(3) and Eq.(4), where
$u_{f}$ is the vector from the reference individual to the nearest stimuli
source, and $u_{r}$ is the vector from the reference individual to the nearest
risk source.

$\displaystyle F_{u_{f}}$ $\displaystyle=$
$\displaystyle\frac{u_{f}}{\left\|u_{f}\right\|},$ (3) $\displaystyle
F_{u_{r}}$ $\displaystyle=$
$\displaystyle\frac{u_{r}}{\left\|u_{r}\right\|}\cdot\begin{cases}\frac{\left\|u_{r}\right\|-\text{r}_{\mathrm{interact}}}{\text{r}_{\mathrm{interact}}}&\text{if}\
\left\|u_{r}\right\|\leq\text{r}_{\mathrm{interact}},\\\
0&\text{Elsewhere}.\\\ \end{cases}$ (4)

The ultimate strategy movement vector $V_{rf}$ of the $rf$ individual fish is
defined as a weighted combination of all the influencing forces, which is
expressed by Eq. 5,

$\displaystyle V_{rf}$ $\displaystyle=$
$\displaystyle\omega_{1}F_{u_{c}}+\omega_{2}F_{u_{a}}+\omega_{3}F_{u_{f}}++\omega_{4}F_{u_{r}}$
(5)

where $\left\\{\omega_{i}\right\\}_{i=1}^{4}$, are the controlling weights
that are based on the Spatiotemporal state of the individual surroundings.

##  3 Domain Decomposition Interaction Simulation Model

Simulating collective behavior can be an extremely time-consuming process,
often taking days, weeks, or even months. To reduce computational time, one
effective solution is to partition the simulation domain into sub-domains.
While domain decomposition may sound simple, its implementation requires
careful consideration of issues that arise at inter-regional boundaries [40].

In this work, we introduce an adaptive and load-balancing domain decomposition
technique to efficiently simulate fish schooling behavior. Our approach
utilizes an unsupervised machine learning algorithm to cluster the simulated
fish into sub-schools based on their spatial-temporal locations. We then
employ Voronoi tessellations to create bounded simulation subdomains that
house the clustered population separately.

To mitigate communication overload between agents on the boundaries of
subdomains and their neighboring subdomains, we replace subdomain-to-subdomain
communication with area-to-area communication. This involves dividing each
subdomain (polytope) into smaller areas based on their vertices, allowing
agents within each area to communicate and interact solely with those in
neighboring areas, rather than with the entire subdomain.

Our implementation approach incorporates several innovative ideas for
leveraging this new algorithm to reduce complexity and enhance the efficient
utilization of computational resources.

###  3.1 Fuzzy c-Means Clustering

Fuzzy c-Means (FCM) is a very popular soft (overlapping) clustering technique.
FCM [41] is a clustering technique allowing each data point to belong to more
than one cluster with association probabilities between 0 and 1 (see Fig. 3).

Figure 3: In this FCM clustering example, each cluster is centered at the
black cross, with members of the same color (red, blue, or green). The
intensity of the cluster color represents the association probability of
members with the cluster. Members closer to the center have a higher
association probability, while those near the boundary have a lower
probability.

The objective is to cluster a set of $N$ data point locations, denoted as
$P=\\{p_{1},x_{2},\ldots,x_{N}\\}$, into $k$ clusters with centroids
$C=\\{c_{1},c_{2},\ldots,c_{k}\\}$. The aim is to minimize the distances
between cluster members while maximizing the distances between clusters. In
other words, FCM formulates this as a minimization problem with the following
objective function:

$J_{m}=\sum_{i=1}^{N}\sum_{j=1}^{k}\mu_{ij}^{m}\|p_{i}-c_{j}\|^{2},$ (6)

where $m>1$ is the fuzzy partition matrix exponent, which controls the degree
of fuzzy overlap, and $\mu_{ij}$ is the probability of associating the point
$x_{i}$ with the $j^{\mbox{th}}$ cluster.

The FCM algorithm [42] works by initially generating random association
probabilities $\mu_{ij}$ of each point, $x_{i}$, and each cluster, $c_{j}$.
The next step involves calculating the centroid of each data partition group
(cluster) as in Eq. (7),

$c_{j}=\frac{\sum_{i=1}^{N}\mu_{ij}^{m}x_{i}}{\sum_{i=1}^{N}\mu_{ij}^{m}},$
(7)

followed by an update of the association probabilities according to Eq. (8),

$\mu_{ij}=\frac{1}{\sum_{r=1}^{k}\bigg{(}\frac{\|x_{i}-c_{j}\|}{\|x_{i}-c_{r}\|}\bigg{)}^{\frac{2}{m-1}}},$
(8)

and calculation of the object function, $J_{m}$. This iterative process
continues until a stopping criterion is met, with the objective function
$J_{m}$ consistently improving. In summary, the FCM algorithm proceeds as
follows: during each iteration, it calculates association probabilities and
cluster centers as per Eq. (7) and Eq. (8), respectively. If the objective
function $J_{m}$ falls below the specified threshold, the process concludes.
Otherwise, the algorithm continues to update $\mu_{i,j}$, $c_{i}$, and $J_{m}$
until the stopping criterion is satisfied.

###  3.2 Voronoi Diagram

Voronoi Diagram (VD) [43, 44] is defined as a set of Voronoi cells
$VD=\bigg{\\{}V(p_{1}),V(p_{2}),...,V(p_{n})\bigg{\\}}$, generated by the
points $p_{1},p_{2},...,p_{n}$, where, given $d$ as the Euclidean distance,
$V(p_{i})=\bigg{\\{}x|\;d(x,p_{i})\leq d(x,p_{j})\;\;\text{for}j\neq
i\bigg{\\}}$
Voronoi cells are non-overlapping, space-filling and the size (volume, area)
of each cell is inversely proportional to the density of the generating points
(see Fig. 4). These properties of Voronoi cells are attractive in high-
dimensional computations [45].

|  ⇨ |
---|---|---
Random points |  | Voronoi diagram
Figure 4: Illustration of Voronoi tessellation

###  3.3 Model Architecture

Our model architecture comprises three key classes: Agent, Region, and Domain.

  * •

The Agent class represents the smallest unit in the simulation and includes
attributes related to its spatial-temporal state.

  * •

The Region class defines a bounded area within the simulation domain that
encompasses agents within its boundaries and neighboring regions. It also
provides functions for managing communication between agents within its range
and those in neighboring regions.

  * •

The Domain class represents a larger section of the simulation composed of
multiple non-overlapping regions. It oversees the allocation of agents to
regions and manages communication between these regions.

#### Region class

In our simulation, we adopt an adaptive approach to region architecture, as
opposed to a predefined "divide and conquer" method. The structure of regions
is determined by the spatial-temporal distribution of agents. Our design for
the region class aims to host an approximately equal number of agents,
mitigating communication overhead in some regions while preventing others from
remaining idle and causing redundant computation. This region construction
process occurs in three steps: clustering, balancing, and partitioning.

  1. 1.

Clustering: We utilize the fuzzy c-Means (FCM) algorithm to group all fish
agents into sub-schools, which occupy temporally non-overlapping regions
within the simulation domain.

  2. 2.

Balancing: We ensure that each cluster has a balanced cardinality, maintaining
an even distribution of agents.

  3. 3.

Partitioning: Finally, we employ the centroids of these balanced clusters to
create non-overlapping Voronoi bounded regions that encompass the clustered
agents (as shown in Fig. 5).

|  ⇨ |  |  ⇨ |
---|---|---|---|---
A large group agents |  | Clustering the agents |  | Partitioning the domain
Figure 5: Illustration of the CP algorithm

#### Communication between agents

Effective communication between agents to share information about their state
and the surrounding environment is crucial in individual-based modeling.
However, it can also significantly contribute to computational costs. Our CP
algorithm and the internal architecture of our region class enable us to
minimize redundant agent-to-agent communication.

We implement two scales of communication per simulation time step:

  1. 1.

Region Scale (School-Scale): Each school within a region informs other schools
in the entire fish population about member agents entering or leaving. Our
adaptive region construction and communication are performed simultaneously
using the CP algorithm.

  2. 2.

Individual Scale: Agents within each region communicate with those in their
neighboring regions. Based on our region class’s internal architecture, each
agent only communicates with agents in the same region and those assigned to
the nearest vertex of the neighboring region. This optimized communication
strategy is illustrated in Fig. 6.

Figure 6: Agent-to-agent communication occurs at two levels within the region.
Internal Communication: Agents within the same region (light blue) exclusively
communicate with one another (blue communication graph); External
Communication: External agents extend their communication to agents in
neighboring regions who share the same vertices.

#### Assign Agent to Region

Our Cluster-Partitioning (CP) algorithm offers several advantages beyond
adaptive region definition. It not only partitions agents into regions but
also identifies inter-region neighbors during the clustering step, eliminating
the need for additional assignment and organization steps.

To minimize communication overhead between agents at region boundaries, we
implement vertex-based partitioning within each region (polytope). This
reduces communication from region-to-region to area-to-area. The
classification of agents as internal or external is based on their interaction
range. Agents with an interaction range confined within their region are
classified as internal, while those reaching neighboring regions are
classified as external (refer to Fig. 7).

This classification is determined by constructing a triangle using the agent’s
location and the nearest two vertices of the region’s corners. We then check
whether the height of the triangle exceeds the interaction radius (as shown in
Fig. 8).

Figure 7: Internal agents are represented by blue filled circles, while
external agents are denoted by red filled circles. The interaction radii of
internal, and external agents are indicated by blue and red circles,
respectively.  Figure 8: The vertices of triangle ABC are such that A
represents the agent’s location, while B and C correspond to the agent’s
nearest two region corners. The classification of an agent as internal or
external, is based on the triangle’s height, $\overline{AD}$.

Equation (9) is mathematical and functional representation of a classifier:

$f(\textit{Area},b)=\begin{cases}Internal&\text{if }\hskip
4.28003pt\text{{Area}}/b\leq R\\\ External&\text{if}\hskip
4.28003pt\text{{Area}}/b>R\\\ \end{cases}$ (9)

where Area and b are the area and base of the triangle ABC, and R is the
radius of interaction.

The second classifier categorizes each external agent based on its nearest
region corner. This classification serves to eliminate redundant communication
overhead from agents in distant neighboring regions, ensuring that agents
communicate only with the essential regions (refer to Fig. 9).

Figure 9: Two agents, represented by red and blue filled circles, are
classified and assigned based on their proximity to region vertices (1 for the
red agent and 3 for the blue agent). These external agents have interaction
ranges extending to neighboring regions, indicated by the light red and purple
regions, respectively

##  4 Simulation experiments

Domain decomposition-based simulation has the potential to significantly
reduce computational costs, and researchers in this field have devoted
substantial efforts to identifying efficient techniques for adaptive workload
redistribution throughout the simulation. Workload can be quantified based on
the number of agents in each subdomain or the computational time spent on each
subdomain, as discussed in [20].

In our simulations, we measured workload in three ways: across the entire
simulation’s duration to assess overall temporal load balancing, within each
subdomain to evaluate spatial load balancing, and at each simulation time step
for temporal load balancing.

Our partitioning algorithm demonstrates superior performance compared to
regular grid-based partitioning. It excels in terms of total running time and
the preservation of both spatial and temporal load balancing. Notably, our
algorithm is versatile and suitable for parallel simulation or distributed
cloud computing, although we exclusively utilized sequential implementations
in our simulations.

All algorithms and simulations were executed on the MATLAB 2020b platform,
operating on a PC equipped with an Intel Core i7, 2.6 GHz CPU, and 16 GB RAM.

In contrast to the static rectangular grid-based approach [10], our proposed
method ensures spatial and temporal load balancing, regardless of the initial
agent distribution and simulation duration (as depicted in Fig. 10).

(a) Spatial distribution of simulated fish

(b) Rectangle-based partitioning

(c) Cluster-based partitioning

Figure 10: Illustration of how ten thousand agents (in a) can be distributed
over 36 subdomains constructed based on grid-based (in b) or cluster-based (in
c) partitioning.

For our initial simulation experiment, we assessed the computational
performance of three implementation methods: no partitioning, regular grid-
based partitioning (square), and partitioning using our CP algorithm. We
conducted four simulations for each method, each involving a different number
of simulated fish.

Table 1 illustrates the results, demonstrating that when the number of
simulated fish is small, partitioning is unnecessary. However, as the number
of simulated fish increases, the non-partitioning method becomes significantly
more time-consuming than the partitioned methods. This disparity arises
because the communication between simulated fish exhibits an order of
$O(N^{2})$, indicating that with a higher number of fish, more time is
required for communication.  N/Type NDD RDD VDD 1e+02 0.17 0.27 2.4 1e+03
13.08 7.85 6.60 1e+04 1.4 e3 718.38 258.48 2e+04 5.51 e+03 2.14e+03 1.02e+03
Table 1:  Computational CPU time costs, measured in seconds, for fish
simulations employing various implementation methods. These simulations
involve different numbers of simulated fish. Specifically, ’NDD’ represents no
domain decomposition, ’RDD’ stands for regular domain decomposition, and ’VDD’
denotes Voronoi domain decomposition, facilitated by the CP algorithm.  Figure
11: A computation cost comparison between regular grids (RDD) and the Voronoi
irregular grids (VDD) domain decomposition implementation. The figures show
the CPU time as a function of the number of subdomains and the size of the
simulated fish school

###  4.1 Load balancing performance comparison: Regular Cells vs. Voronoi
Cells

We conducted a load balancing performance comparison between our Voronoi-based
domain decomposition (VDD) and the regular grid-based domain decomposition
(RDD) methods. As depicted in Fig. 12, the VDD implementation excels in
achieving a well-balanced workload distribution among subdomains, resulting in
a lower variance of CPU time across subdomains when compared to the RDD
implementation. This advantage is primarily attributed to VDD’s adaptability
in adjusting subdomain size and shape based on the spatial distribution of
fish, which promotes an even workload distribution. In contrast, RDD divides
the simulation domain into equal-sized regular grid cells, which may not
accurately reflect the fish’s spatial distribution and can lead to subdomains
with significantly different fish counts, resulting in load imbalance.

In summary, our results clearly demonstrate that the Voronoi-based domain
decomposition approach proves to be a more effective and efficient method for
achieving load balance in simulations of fish behavior. This improvement holds
the potential to yield significant enhancements in computational performance.

Likewise, the graphs in Fig. 13 illustrate the superior load balancing
performance of our approach compared to the regular grid cells (RDD)
implementation within each subdomain of the simulation. In comparison to the
RDD approach (as shown in Fig. 13a), our approach exhibits significantly
reduced variation in the computational cost of simulation within each
subdomain, aligning more closely with the average workload (as seen in Fig.
13b).

Figure 12: The computation cost profile of the fish simulation is depicted for
both RDD (solid lines) and VDD (dashed lines) implementations. The stippled
lines represent the CPU time for 100 time steps, while the bold lines denote
the average CPU time.

(a)

(b)

Figure 13: The graph shows the mean (represented by the curve) and standard
deviation (indicated by vertical lines) of the computational cost per region
per iteration. This simulation involves 5,000 agents.

##  5 Conclusions and Future Work

We have successfully implemented a large-scale simulation of fish schooling
behavior using Voronoi tessellations generated by fuzzy c-Mean cluster
centroids. Our approach effectively addresses common challenges in agent-based
simulations, including communication overload and workload imbalances. Our
simulations have shown that our partitioning algorithm maintains spatial and
temporal load balance while significantly reducing computational costs
compared to other methods.

Looking ahead, our future work will involve applying this approach to simulate
fish migratory behaviors, incorporating environmental variables like
temperature, salinity, and ocean currents. Given the varying resolutions of
these variables, this simulation will be more complex, necessitating efficient
domain decomposition techniques to mitigate communication overload and high
computational costs.

## References

  * [1] Andreas Huth and Christian Wissel  “The simulation of the movement of fish schools”  In _Journal of theoretical biology_ 156.3 Elsevier, 1992, pp. 365–385
  * [2] CARL-OSCAR Erneholm  “Simulation of the flocking behavior of birds with the boids algorithm”  In _Royal Institute of Technology_ , 2011
  * [3] Jade Alglave, Luc Maranget and Michael Tautschnig  “Herding cats: Modelling, simulation, testing, and data mining for weak memory”  In _ACM Transactions on Programming Languages and Systems (TOPLAS)_ 36.2 ACM New York, NY, USA, 2014, pp. 1–74
  * [4] William H Warren  “Collective motion in human crowds”  In _Current directions in psychological science_ 27.4 SAGE Publications Sage CA: Los Angeles, CA, 2018, pp. 232–240
  * [5] Théo Maire and Hyun Youk  “Molecular-level tuning of cellular autonomy controls the collective behaviors of cell populations”  In _Cell systems_ 1.5 Elsevier, 2015, pp. 349–360
  * [6] John D Treado et al.  “Bridging particle deformability and collective response in soft solids”  In _Physical Review Materials_ 5.5 APS, 2021, pp. 055605
  * [7] Melanie Schranz, Martina Umlauft, Micha Sende and Wilfried Elmenreich  “Swarm robotic behaviors and current applications”  In _Frontiers in Robotics and AI_ 7 Frontiers Media SA, 2020, pp. 36
  * [8] Paul F Egan, Jonathan Cagan, Christian Schunn and Philip R LeDuc  “Design of complex biologically based nanoscale systems using multi-agent simulations and structure–behavior–function representations”  In _Journal of Mechanical Design_ 135.6 American Society of Mechanical Engineers, 2013, pp. 061005
  * [9] Roberto Solar, Remo Suppi and Emilio Luque  “High performance distributed cluster-based individual-oriented fish school simulation”  In _Procedia Computer Science_ 4 Elsevier, 2011, pp. 76–85
  * [10] Guillermo Vigueras, Miguel Lozano, Juan Manuel Orduna and Francisco Grimaldo  “A comparative study of partitioning methods for crowd simulations”  In _Applied Soft Computing_ 10.1 Elsevier, 2010, pp. 225–235
  * [11] Roberto Solar, Remo Suppi and Emilio Luque  “Proximity load balancing for distributed cluster-based individual-oriented fish school simulations”  In _Procedia Computer Science_ 9 Elsevier, 2012, pp. 328–337
  * [12] Chahrazed Labba, Narjes Bellamine Saoud and Julie Dugdale  “Towards a conceptual framework to support adaptative agent-based systems partitioning”  In _2015 IEEE/ACIS 16th International Conference on Software Engineering, Artificial Intelligence, Networking and Parallel/Distributed Computing (SNPD)_ , 2015, pp. 1–5  IEEE
  * [13] Lamia Youseff et al.  “Parallel modeling of fish interaction”  In _2008 11th IEEE International Conference on Computational Science and Engineering_ , 2008, pp. 234–241  IEEE
  * [14] Gennaro Cordasco, Vittorio Scarano and Carmine Spagnuolo  “Distributed mason: A scalable distributed multi-agent simulation environment”  In _Simulation Modelling Practice and Theory_ 89 Elsevier, 2018, pp. 15–34
  * [15] Samuel J Araki and Robert S Martin  “Dynamic load balancing with over decomposition in plasma plume simulations”  In _Journal of Parallel and Distributed Computing_ 163 Elsevier, 2022, pp. 136–146
  * [16] Dongliang Zhang, Changjun Jiang and Shu Li  “A fast adaptive load balancing method for parallel particle-based simulations”  In _Simulation Modelling Practice and Theory_ 17.6 Elsevier, 2009, pp. 1032–1042
  * [17] Roberto Solar, Remo Suppi and Emilio Luque  “High performance individual-oriented simulation using complex models”  In _Procedia Computer Science_ 1.1 Elsevier, 2010, pp. 447–456
  * [18] Bo Zhou and Suiping Zhou  “Parallel simulation of group behaviors”  In _Proceedings of the 2004 Winter Simulation Conference, 2004._ 1, 2004  IEEE
  * [19] MA Stijnman, RH Bisseling and GT Barkema  “Partitioning 3D space for parallel many-particle simulations”  In _Computer physics communications_ 149.3 Elsevier, 2003, pp. 121–134
  * [20] Maria Serg Egorova, Sergey A Dyachkov, Anatoly N Parshikov and VV Zhakhovsky  “Parallel SPH modeling using dynamic domain decomposition and load balancing displacement of Voronoi subdomains”  In _Computer Physics Communications_ 234 Elsevier, 2019, pp. 112–125
  * [21] Ji-Gui Sun, Jie Liu and Lian-Yu Zhao  “Clustering algorithms research”  In _Journal of software_ 19.1, 2008, pp. 48–61
  * [22] Shi Na, Liu Xumin and Guan Yong  “Research on k-means clustering algorithm: An improved k-means clustering algorithm”  In _2010 Third International Symposium on intelligent information technology and security informatics_ , 2010, pp. 63–67  Ieee
  * [23] Yongwei Wang et al.  “Cluster based partitioning for agent-based crowd simulations”  In _Proceedings of the 2009 Winter Simulation Conference (WSC)_ , 2009, pp. 1047–1058  IEEE
  * [24] Youssef M Marzouk and Ahmed F Ghoniem  “K-means clustering for optimal partitioning and dynamic load balancing of parallel hierarchical N-body simulations”  In _Journal of Computational Physics_ 207.2 Elsevier, 2005, pp. 493–528
  * [25] Moritz Looz, Charilaos Tzovas and Henning Meyerhenke  “Balanced k-means for parallel geometric partitioning”  In _Proceedings of the 47th International Conference on Parallel Processing_ , 2018, pp. 1–10
  * [26] Botao Zhu et al.  “Improved soft-k-means clustering algorithm for balancing energy consumption in wireless sensor networks”  In _IEEE Internet of Things Journal_ 8.6 IEEE, 2020, pp. 4868–4881
  * [27] Chunqiong Wu et al.  “k-Means Clustering Algorithm and Its Simulation Based on Distributed Computing Platform”  In _Complexity_ 2021 Hindawi, 2021
  * [28] Pengcheng Shen and Chunguang Li  “Distributed information theoretic clustering”  In _IEEE Transactions on Signal Processing_ 62.13 IEEE, 2014, pp. 3442–3453
  * [29] Craig W Reynolds  “Flocks, herds and schools: A distributed behavioral model”  In _Proceedings of the 14th annual conference on Computer graphics and interactive techniques_ , 1987, pp. 25–34
  * [30] Ugo Erra, Rosario De Chiara, Vittorio Scarano and Maurizio Tatafiore  “Massive simulation using gpu of a distributed behavioral model of a flock with obstacle avoidance”  In _Proceedings of Vision, Modeling and Visualization 2004 (VMV)_ , 2004
  * [31] Vu Phi Tran, Matthew A Garratt and Ian R Petersen  “Switching formation strategy with the directed dynamic topology for collision avoidance of a multi-robot system in uncertain environments”  In _IET Control Theory & Applications_ 14.18 IET, 2020, pp. 2948–2959
  * [32] Martin Barksten and David Rydberg  “Extending Reynolds’ flocking model to asimulation of sheep in the presence of a predator”, 2013
  * [33] Gabriel Chang and Michaela Stjerndal  “Investigating and Modeling the Emergent Flocking Behaviour of Sheep Under Threat with Fear Contagion”, 2019
  * [34] Christopher Hartman and Bedrich Benes  “Autonomous boids”  In _Computer Animation and Virtual Worlds_ 17.3-4 Wiley Online Library, 2006, pp. 199–206
  * [35] Saleh Alaliyat, Harald Yndestad and Filippo Sanfilippo  “Optimisation Of Boids Swarm Model Based On Genetic Algorithm And Particle Swarm Optimisation Algorithm (Comparative Study).”  In _ECMS_ , 2014, pp. 643–650  Citeseer
  * [36] Yen-Wei Chen, Kanami Kobayashi, Xinyin Huang and Zensho Nakao  “Genetic algorithms for optimization of boids model”  In _International Conference on Knowledge-Based and Intelligent Information and Engineering Systems_ , 2006, pp. 55–62  Springer
  * [37] Andriy Dmytruk et al.  “Safe Tightly-Constrained UAV Swarming in GNSS-denied Environments”  In _2021 International Conference on Unmanned Aircraft Systems (ICUAS)_ , 2021, pp. 1391–1399  IEEE
  * [38] Alethea BT Barbaro et al.  “Discrete and continuous models of the dynamics of pelagic fish: application to the capelin”  In _Mathematics and Computers in Simulation_ 79.12 Elsevier, 2009, pp. 3397–3414
  * [39] Michele Ballerini et al.  “Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study”  In _Proceedings of the national academy of sciences_ 105.4 National Acad Sciences, 2008, pp. 1232–1237
  * [40] UR Hanebutte and Adrian Michel Tentner  “Traffic simulations on parallel computers using domain decomposition techniques”, 1995
  * [41] Hui Yang et al.  “Resource assignment based on dynamic fuzzy clustering in elastic optical networks with multi-core fibers”  In _IEEE Transactions on Communications_ 67.5 IEEE, 2019, pp. 3457–3469
  * [42] James C Bezdek  “Pattern recognition with fuzzy objective function algorithms”  Springer Science & Business Media, 2013
  * [43] B Boots, A Okabe and K Sugihara  “Spatial tessellations”  In _Geographical information systems_ 1 John Wiley & Sons New York, NY, 1999, pp. 503–526
  * [44] Wojciech Pokojski and Paulina Pokojska  “Voronoi diagrams–inventor, method, applications”  In _Polish Cartographical Review_ 50.3, 2018, pp. 141–150
  * [45] Sam Subbey, Christie Mike and Malcolm Sambridge  “A strategy for rapid quantification of uncertainty in reservoir performance prediction”  In _SPE Reservoir Simulation Symposium_ , 2003  OnePetro