# How to Implement Continuous-Time Multi-Agent Crowd Simulation

## Table of Contents

## Introduction

In this article, we discuss the mathematics on how to implement a continuous-time multi-agent crowd simulation. The article covers two agent models, circular, unorientable model, and three-circle, orientable model. We include extensive discussion about simulation objects, simulation mechanics, and implementation details and conclude by addressing the shortcoming and challenges of this type of simulation. We intended this article as a reference implementation for crowd simulation. We have experimented with the implementation in the crowddynamics library.

We define *crowd dynamics* as the movement of crowds of people in respect of time. We implemented the simulation model using a continuous-time, Newtonian mechanics with an *agent-based* model, where each agent is a rigid body, in a two-dimensional lattice. In the model, the agents are intelligent self-propelled particles, whose movement is produced by non-physical and physical forces. The non-physical forces are referred to as *social forces*, for example, the desire to move toward a goal or collision avoidance. Hence, we refer to the model as the *social force model*, a term originating from ^{1}.

There also exist other types of models, such as *cellular automata* which models agents as an occupied or unoccupied cells in a discrete lattice, that is, a grid, and *continuum-based* models which model the crowd as a continuous mass. The agent-based model has the advantage of being the most realistic model, but it is also computationally the most challenging model.

There are real-life applications for research on how crowds of people move, for example:

- Improving the safety of buildings, mass events, and crowded places.
- Estimating evacuation times and improving the evacuation process.
- Venue design and improving the flow rate of the crowd in crowded places, such as subways.
- Producing realistic movement in games, animations, and virtual reality.

The authors of the book ^{2} give a comprehensive look at the subject of crowd dynamics.

In the following sections, we are going to explain in-depth the mathematics behind the agent-based model and how to implement the simulation.

## Geometric Primitives

We define two-dimensional geometric primitives as follows:

*Point*is denoted as $๐ฉโโ^2.$*Circle*is defined as a pair $(๐ฉ, r)$ where $๐ฉ$ is the center, and $r$ is the radius.*Line segment*is defined as a pair of points $(๐ฉ_1,๐ฉ_2)$.*Polygonal chain*is a sequence of points $(๐ฉ_1,…,๐ฉ_m),mโฅ2$ where each pair $(๐ฉ_i,๐ฉ_{i+1})$ for $i=1,…,m-1$ is a line segment.*Closed polygonal chain*is a polygonal chain where the endpoints $(๐ฉ_1,๐ฉ_m)$ form a line segment.*Polygon*is a filled shape whose exterior is a closed polygonal chain.

We use these geometric primitives for implementing some of the objects in the simulation.

## Objects

We will refer to the collection of objects that exists in the simulation as the *field*. The field consists of the following objects:

The

*domain*$ฮฉโโ^2$, a subset of the lattice which contains all the other objects.The set of

*agents*$\mathcal{A}_iโฮฉ$ for $iโA$ where $A={1,…,n}$ are the indices. An agent is an impassable, dynamic object. We will denote two distinct agents using the indices as $i,jโA$ where $iโ j$.The set of

*obstacles*$\mathcal{O}โฮฉ$. An obstacle is an impassable, static object. We denote obstacles using the subscript $wโW$, that is, $\mathcal{O}_w โ \mathcal{O}$.The set of

*targets*$\mathcal{T}โฮฉ$. A target is a passable, static object.The set of

*sources*$\mathcal{S}โฮฉ$. A source is a passable, static object. Initially, we place the agents inside a source without overlapping each other. For example, we can sample agent positions inside the source from a uniform distribution and whether the agent overlaps with other agents. If not, place the agent; otherwise, we will resample.

We implement the obstacles using polygonal chains and the domain, targets, and sources using polygons. We describe the implementation of the agent’s shape in the Agents section.

## Agents

In the agent-based model, each agent is a discrete object, modeled as a rigid body. We use two models for modeling the shape of the agent; circular and three-circle, models. In addition to shape, agents have other attributes such as mass and desired velocity.

### Circular Model

In the *circular model*, we model the agent’s shape using a circle. It is the simplest model and rotationally symmetrical. It is easiest to implement and compute and work well for sparse and medium-density crowds. However, it is not a very realistic model of the cross-section of the human torso. In reality, the human body is narrower and can fit through smaller spaces. In dense crowds, the circle model can result in artifacts and unrealistic density measures.

Each circular agent $iโA$ has the following attributes.

- $๐ฑ_i$
*center of mass* - $๐ฏ_i$
*velocity* - $๐_i$
*force* - $\hat{๐}_i$
*desired direction* - $v_i$
*desired velocity* - $r_i$
*radius* - $m_i$
*mass*

### Three-circle Model

The *three-circle model* consists of three circles representing the shoulders and midsection. The three-circle model is a more realistic representation of the human cross-section. However, since it is not rotationally symmetrical, simulation mechanics become more complex. For example, the model must include rotational motion and more complicated collision detection. In dense crowds, agents can move sideways, which requires solving desired body rotation. ^{3} ^{4} ^{5} ^{6}

Each orientable agent has following additional attributes:

- $ฯ_i$
*orientation* - $ฯ_i$
*angular velocity* - $M_i$
*torque* - $ฯ_i^0$
*desired orientation* - $ฯ_i^0$
*desired angular velocity* - $I_i$
*rotational moment*

We constrained the values of the orientations $ฯ$ and $ฯ^0$ into interval $[-ฯ,ฯ].$

In addition to the attributes of the circular and rotatable agent, each three-circle agent has attributes as follows.

- $r_i^{t}$
*torso radius* - $r_i^{s}$
*shoulder radius* - $r_i^{ts}$
*torso-to-shoulder radius*

Identity $r_i^{ts} + r_i^s = r_i$ holds for these attributes.

The following sections will discuss the simulation mechanics, including forces and torques. We cover the mechanics for the circular model in detail. The forces for the circular model can be generalized for the three-circle agent model by deriving the variables, social force, and contact force between three-circle agents and three-circle agent and line obstacle. However, for simplicity reasons, we have left the generalization out from this article for now.

## Variables

We define the *relative position* $\tilde{๐ฑ}$, *relative velocity* $\tilde{๐ฏ}$, *skin-to-skin distance* $h$ (can be negative), *normal (unit) vector* $\hat{๐ง}$ and *tangent (unit) vector* $\hat{๐ญ}$. We use these variables for computing the interaction potentials. We use $R(ฯ)$ to define the rotation matrix in two dimensions where $ฯโ[0ยฐ,360ยฐ]$ is the angle in degrees.

### Between Two Distinct Circular Agents

Between two distinct circular agents $i,jโA,iโ j$ we have the variables as follows: Relative position $\tilde{๐ฑ}=๐ฑ_j-๐ฑ_i$, relative velocity $\tilde{๐ฏ}=๐ฏ_j-๐ฏ_i$, normal vector $\hat{๐ง}=\tilde{๐ฑ}/|\tilde{๐ฑ}|$, tangent vector $\hat{๐ญ}=R(270ยฐ)โ \hat{๐ง}$ and skin-to-skin distance $h=|\tilde{๐ฑ}|-(r_j+r_i).$

### Between Circular Agent and Line-obstacle

Between a circular agent $i$ and static line-obstacle $w$ defined as line segment $(๐ฉ_1,๐ฉ_2)$, we have following auxialiary variables: Line vector $\tilde{๐ฉ} = ๐ฉ_2-๐ฉ_1$, line length $l = |\tilde{๐ฉ}|$, line tangent $๐ญ_w = \tilde{๐ฉ}/l$, line normal $๐ง_w = R(90ยฐ) โ ๐ญ_w$, two vectors from the agent’s center of mass to the line points $๐ช_1 = ๐ฉ_1 - ๐ฑ_i$ and $๐ช_2 = ๐ฉ_2 - ๐ฑ_i$, and tangential distance of agent’s center of mass from line points $l_t = ๐ญ_wโ (๐ช_1 + ๐ช_2).$

Now, we determine the relative position of the center of mass compared to the line segment. We have three distinct cases:

- $l_t > l$: Center of mass $๐ฑ_i$ is on the left of the line perpendicular to $\tilde{๐ฉ}$ intersecting $๐ฉ_1$.
- $l_t < -l$: Center of mass $๐ฑ_i$ is on the right of the line perpendicular to $\tilde{๐ฉ}$ intersecting $๐ฉ_2$.
- Otherwise: The center of mass $๐ฑ_i$ is between the perpendicular lines.

We define the variables as follows: Skin-to-skin distance is

Normal vector is

Tangent vector $\hat{๐ญ}=R(270ยฐ)โ \hat{๐ง}$, and relative velocity $\tilde{๐ฏ}=๐ฏ_i-๐ฏ_w=๐ฏ_i.$ We do not need the relative position $\tilde{๐ฑ}$ between agent and line-obstacle.

## Total Force

The total force on an agent $iโA$ consists of the adjusting force, physical contact forces, social forces and fluctuation force. We define it as $$ ๐_i = ๐_{i}^{adj} + โ_{jโA,jโ i} (๐_{i,j}^{soc} + ๐_{i,j}^{con}) + โ_{wโW} ๐_{i,w}^{con} + \mathbf{ฮพ}_i. \label{eq:1} \tag{1} $$ This section discuss about computing each of these constituents.

### Adjusting Force

The adjusting force accounts for the agent’s desire to reach its destination. It works by adjusting the total force to align with the desired direction $\hat{๐}$ in the characteristic time $ฯ_{adj}$. We define it as $$ ๐_{adj} = \frac{m}{ฯ_{adj}} (v \hat{๐} - ๐ฏ), \label{eq:2} \tag{2} $$ where $m$ is the mass, $v$ is the desired velocity, $๐ฏ$ is the current velocity, and $ฯ_{adj}$ is the characteristic time.

### Contact Force

The contact force account for the physical contact with agents and other objects. We define the contact force as the sum of *tangential friction*, *counter compressive*, and *damping* forces
$$
๐_{con} =
\begin{cases}
h ฮบ (\tilde{๐ฏ} โ
\hat{๐ญ}) \hat{๐ญ} - h ฮผ \hat{๐ง} + ฮณ (\tilde{๐ฏ} โ
\hat{๐ง}) \hat{๐ง} & hโค0 \\

0, & h>0
\end{cases}
\label{eq:3}
\tag{3}
$$
where $ฮบ$ is tangential friction coefficient, $ฮผ$ is the counter compressive coefficient, $ฮณ$ is the damping coefficient, $\hat{๐ง}$ is the normal unit vector and $\hat{๐ญ}$ is the normal tangent vector. It follows the original formula used by ^{1} with the addition of the damping by ^{3}.

It is important to note that the friction model models the whole friction between two objects. Therefore, if an agent touches multiple line obstacles, we must compute the force as a weighted average of the contact forces between all of the line-obstacles. Otherwise, the model would give inconsistent results for friction between similar objects. For example, if we compute we naively computed the friction between the agent and line obstacle, and then split the line obstacle into two line-obstacles in the point of contact, the model would give twice as high friction between the split obstacle.

### Fluctuation Force

Fluctuation force is a stochastic force analogous to heat in particle systems. It accounts for the lateral movement of the agents and breaks symmetries in the simulation. $$ \mathbf{ฮพ}=ฮพ\hat{๐ฎ}(ฯ) \label{eq:4} \tag{4} $$

We draw the magnitude of the fluctuation force from a truncated normal distribution $ฮพโผN(0, ฯ_ฮพ^2)$ where $ฯ_ฮพ^2$ is the variance and the direction from a uniform distribution $ฯโผU(0, 2ฯ)$. The normal distribution $N$ is truncated to three standard deviations. We define $\hat{๐ฎ}(ฯ)$ to be a unit vector to direction $ฯ$.

### Social Force

The social force between agents accounts for collision avoidance. The original social force model used social force only dependent on the distance between agents. However, ^{7} ^{8} introduced a more realistic, statistical mechanics-based approach for modeling this interaction. The interaction depends on the linearly estimated future positions on the agents, which takes into account the positions and velocities of the agents. The authors derived it from a pair-wise interaction potential that follows the power law, and it is backed up by experimental data. The potential is defined
$$
E(ฯ) = m \frac{k}{ฯ^2} \exp \left( -\frac{ฯ}{ฯ_{soc}} \right),
\label{eq:5}
\tag{5}
$$
where $k$ is a social force scaling coefficient and $ฯ_{soc}$ is the interaction time horizon. The time-to-collision $ฯ$ is obtained by linear extrapolation of the current trajectories of the agents and solving the root for *skin-to-skin* distance $h$ as a function of $ฯ$ between these agents. The resulting $ฯ$ is a function of the relative position $\tilde{๐ฑ}$. We can derive the corresponding force by taking the spatial gradient of the potential
$$
๐_{soc} = -โ_{\tilde{๐ฑ}} E(ฯ) = - m k \left(\frac{1}{ฯ^2}\right) \left(\frac{2}{ฯ} + \frac{1}{ฯ_{soc}}\right) \exp\left (-\frac{ฯ}{ฯ_{soc}}\right ) โ_{\tilde{๐ฑ}} ฯ.
\label{eq:6}
\tag{6}
$$
If $ฯ$ is not defined, the force $๐_{soc}$ is equal to zero.

For two circular agents with radius $r_i$ and $r_j$ the skin-to-skin distance is the distance between their positions subtracted by the sum of the radii
$$
\begin{aligned}
h(ฯ) &= |(๐ฑ_i + ฯ ๐ฏ_i) - (๐ฑ_j + ฯ ๐ฏ_j)| - (r_i+r_j) \\

&= |\tilde{๐ฑ} + ฯ \tilde{๐ฏ}| - (r_i+r_j).
\end{aligned}
$$
We receive the quadratic equation in terms of $ฯ$ by settings $h(ฯ)=0$
$$
\begin{aligned}
|\tilde{๐ฑ} + ฯ \tilde{๐ฏ}| - (r_i+r_j) &= 0 \\

|\tilde{๐ฑ} + ฯ \tilde{๐ฏ}|^2 &= (r_i+r_j)^2 \\

|\tilde{๐ฏ}|^2 ฯ^2 + (\tilde{๐ฑ}โ
\tilde{๐ฏ}) ฯ + |\tilde{๐ฑ}|^2 - (r_i+r_j)^2 &= 0.
\end{aligned}
$$
We solve the quadratic equation
$$
ฯ = \frac{b-\sqrt{b^2-ac}}{a},
\label{eq:7}
\tag{7}
$$
where $a=|\tilde{๐ฏ}|^2$, $b=-\tilde{๐ฑ}โ
\tilde{๐ฏ}$, and $c=|\tilde{๐ฑ}|^2-(r_i+r_j)^2$. Now we can solve the gradient
$$
โ_{\tilde{๐ฑ}} ฯ = \frac{1}{a} \left(\tilde{๐ฏ} - \frac{a \tilde{๐ฑ}+b\tilde{๐ฏ}}{\sqrt{b^2-a c}}\right).
\label{eq:8}
\tag{8}
$$
By substituting the gradient to the social force, we can use it in the simulation.

## Total Torque

The total torque on agent $iโA$ consists of the adjusting, contact and fluctuation torque. We define it as $$ M_i = M_i^{adj} + โ_j M_{i,j}^{con} + โ_w M_{i,w}^{con} + ฮท. \label{eq:9} \tag{9} $$ This section discuss about computing each of these constituents.

### Adjusting Torque

The adjusting torque accounts for the agent’s desire to rotate its body to aligh with the desired orientation $ฯ_0$ in the characteristic time $ฯ_{rot}$. We define it as $$ M_{adj} = \frac{I}{ฯ_{rot}} \left(\frac{w_ฯ(ฯ_0-ฯ)}{ฯ} ฯ_0 - ฯ\right) \label{eq:10} \tag{10} $$

The function $w_ฯ$ wraps an arbitrary angle $ฯ$ in radians to the interval $(-ฯ,ฯ]$. The function $w_ฯ$ is defined as
$$
\begin{aligned}
w_ฯ(ฯ) &= w_ฯ'(ฯ \mod 2ฯ) \\

w_ฯ'(\hat{ฯ}) &=
\begin{cases}
\hat{ฯ} & \hat{ฯ}โคฯ \\

\hat{ฯ}-2ฯ & \hat{ฯ}>ฯ
\end{cases}.
\end{aligned}
\label{eq:11}
\tag{11}
$$

We can set the desired orientation $ฯ_0$ to walking direction. In dense crowds, the realistic behavior of the desired orientation might be different. For example, an agent might try to squeeze sideways through two agents. However, we have not implemented algorithms for modeling such behavior.

### Contact Torque

The contact torque account for the torque caused by the contact force $๐_{con}$. We define it as $$ M_{con} = ๐ซ ร ๐_{con}, \label{eq:12} \tag{12} $$ where $๐ซ$ is the contact position vector. The contact position vector is a vector from the agents center of mass $๐ฑ$ to the contact point.

The cross product between two $2$-dimensional vectors $๐ฎ$ and $๐ฏ$ is defined as $$ ๐ฎร๐ฏ = u_0 v_1 - u_1 v_0. \label{eq:13} \tag{13} $$ It is the $z$-component of $3$-dimensional cross product of vectors $๐ฎ$ and $๐ฏ$ with $z$-components that are zero.

### Fluctuation Torque

The fluctuation torque accounts for the stochastic fluctuation in the agent’s orientation. We define it as $$ ฮท = I ฮถ, \label{eq:14} \tag{14} $$ where the magnitude is drawn from a truncated normal distribution $ฮถ โผ N(0, ฯ_ฮถ^2)$ with variance $ฯ_ฮถ^2.$ The normal distribution is truncated to three standard deviations.

## Pathfinding and Obstacle Avoidance

In this simulation, we handle the psychological avoidance of obstacles by altering the agent’s desired direction away from obstacles if they get too close. We handle the pathfinding by assuming that agents follow the shortest path to their targets, taking into account the finite size of the agents. The approach is computationally efficient since we have to compute the shortest path only once. ^{9} ^{10}

### Continuous Shortest Path

The solution to *eikonal equation* gives us the minimal amount of time required to travel from $x$ to the boundary $โฮฉ$. The eikonal equation is defined as
$$
\begin{aligned}
|S(๐ฑ)| &= \frac{1}{f(๐ฑ)},\quad ๐ฑโฮฉ, \\

S(๐ฑ) &= q(๐ฑ),\quad ๐ฑโโฮฉ
\end{aligned}
\label{eq:15}
\tag{15}
$$
where $f:ฮฉโ(0,+โ)$ is the *speed of travel* and $q:โฮฉโ[0,+โ)$ is the *exit-time penalty*. The solution with a constant speed of travel $f$ is a continuous shortest path, refered as a *distance map*.

We will also define a unit-vector field, referred as a *direction map*, which points towards the boundary with the lowest exit-time penalty, as the normalized gradient
$$
D(S(๐ฑ))=-\frac{โS(๐ฑ)}{| โS(๐ฑ) |}.
\label{eq:16}
\tag{16}
$$

We can solve the Eikonal equation using existing algorithms, such as *Fast Marching Method*.

### Desired Direction

We define a *avoidance radius* as $r>0$ and *buffered obstacles* $๐'$ as the set of points at a distance less or equal to $r$ from obstacles $๐.$ We define the *desired direction* as a weighted average of two direction maps, one pointing to the target and one away from obstacles such that the weight depends on the proximity from the obstacles. When agents are closer to the obstacles, they desire more strongly away from them, and when they are further away from obstacles, they desire more strongly towards targets.

Let $T(๐ฑ)$ be a distance map to target $๐ฃ$ with domain $ฮฉ$ and exit-time penalty $T(๐ฑ) = 0,, ๐ฑ โ ๐ฃ$ and $T(๐ฑ) โ โ,, ๐ฑ โ ๐'$. Let $O(๐ฑ)$ be the distance map to obstacles $๐$ with domain $ฮฉ$ and exit-time penalty $O(๐ฑ) = 0,, ๐ฑ โ ๐$. The desired direction is then defined as $$ \hat{๐}(๐ฑ) = ฮป(O(๐ฑ)) D(T(๐ฑ)) + (1 - ฮป(O(๐ฑ))) D(O(๐ฑ)), \label{eq:17} \tag{17} $$ where $ฮป: โ^+ โ [0, 1]$ is decreasing function, $\frac{d}{dx} ฮป(x)โค0$, such that $ฮป(0) = 1$ and $\lim_{xโโ} ฮป(x) = 0$.

We will set the avoidance radius $r$ sufficiently larger than $\max_{iโA} r_i$, which ensures that agents tend away from the obstacles.

### Exponential ฮป-function

For example, we can use the exponential function, defined as $ฮป(x) = ab^{c x},$ with coefficients $a,b,c$ where $b>0$, which we can solve from the properties of function $ฮป$. Additionally, we require $ฮป(r)=ฯต>0$ such that $ฯต$ is close to zero.

- $ฮป(0)=a=1$
- $ฮป(r)=b^{c r}=ฯต$ gives us $c=\frac{\log_b ฯต}{r}$ therefore $ฮป(x)=ฯต^{x/r}$
- $\lim_{xโโ} ฮป(x) = 0$ implies $ฯต<1$
- Lastly, we verify that the function is decreasing, that is, $\frac{d}{dx} ฮป(x)=\frac{\log ฯต}{r} ฯต^{x/r}<0.$ Note that $\log ฯต<0$ because $0<ฯต<1.$

The resulting function $ฮป$ is defined as
$$
ฮป(x)=ฯต^{x/r},
\label{eq:18}
\tag{18}
$$
where $0<ฯต<1$ is referred as *strength* and $r>0$ is the avoidance radius.

### Piecewise Linear ฮป-function

We can also use a piecewise linear function such that $ฮป(r)=0$. Then we have
$$
ฮป(x)=
\begin{cases}
-x/r+1 & xโคr \\

0 & x>r
\end{cases}.
\label{eq:19}
\tag{19}
$$
This formulation is also a logical choice since the smoothness of the exponential formulation does not bring any additional value in the agent-based simulation.

## Computing Interactions

In the simulation, we have two types of *interactions*, agent-agent and agent-obstacle interactions. Interactions are computationally the most challenging part of the simulation, and therefore it is worth-while to optimize how we compute them.

The equations $\eqref{eq:1}$ and $\eqref{eq:9}$ express the interactions as the summation of the interaction potentials, that is, contact and social forces and torques. Both interactions depend on the distance between the objects.

Contact interactions only take place if the objects collide with each other.

Social force approaches zero at exponential speed as time-to-collision increases, that is, as the distance between agents increases. Therefore, we can define a cutoff distance $r_{c}โฅ0$ such that when $h>r_{c}$, the social force is close enough to zero that we can approximate it as zero.

We will refer to these interactions that have a finite range as *local interactions*. For local interactions, computing the total interaction potential requires summing the interaction potentials with the pairs of objects within the interaction radius $r_{c}โฅ0$. Fortunately, there exist efficient methods for finding such pairs of objects, which we will briefly discuss in this section.

### Agent-Agent Interactions

To find neighboring agent pairs within distance $r$, we can use *fixed-radius nearest neighbor search*. In particular, we use the *cell lists* algorithm, which improves the difficulty over brute force, from $O(n^2)$ to $O(n).$ Agents are required to have two properties:

- Agent-agent interactions are local.
- Agents have a finite size. Therefore, only a finite amount of agents can fit inside a finite area.

In practice, Cell Lists is better once the number of agent $n$ is large enough due to the size of the constant term in the $O$-notation and overhead of the cell lists algorithm. The exact size of $n$ when Cell lists is better must be determined experimentally. We have implemented the Cell lists algorithm in the link.

### Agent-Obstacle Interactions

Implementing efficient collision detection between agents and obstacles is not covered in this article. We recommend looking at spatial hashing techniques or how other physics engines implement collision detection.

## Updating the System

In order to update the system, we need to solve the second-order differential equations numerically. The translational motion is defined by the equation $$ m \frac{d^2}{dt^2} ๐ฑ(t) = ๐(t). \label{eq:20} \tag{20} $$

The rotational motion is defined by the equation $$ I \frac{d^2}{dt^2} ฯ(t) = M(t). \label{eq:21} \tag{21} $$

Because computing the total force on the agents is computationally expensive, we will use the first-order integrator. We chose to use *Verlet integration* with *adaptive timestep* $ฮt$.

### Adaptive Timestep

Adaptive timestep is selected from interval $[ฮt_{min}, ฮt_{max}]$ such that we limit the maximum distance that an agent can move in each iteration. First, we define
$$
ฮt_{mid} = ฮt_{max} \frac{\max_{iโA} v_i}{\max_{iโA} |๐ฏ_i|}
\label{eq:22}
\tag{22}
$$
where $|๐ฏ_i|$ is agent’s current velocity and $v_i$ is agent’s target velocity. Then, the adaptive timestep is defined
$$
ฮ t =
\begin{cases}
ฮt_{min} & ฮt_{mid} < ฮt_{min} \\

ฮt_{mid} & \mathrm{otherwise} \\

ฮt_{max} & ฮt_{mid} > ฮt_{max}
\end{cases}.
\label{eq:23}
\tag{23}
$$

The adaptive timestep helps to avoid numerical errors in the simulation.

### Verlet Integration

Verlet integration algorithm updates the new velocity and position at iteration $k=0,1,…$ as follows

$$
\begin{aligned}
๐ฏ'&=๐ฏ_k+๐_k ฮt/2 \\

๐_{k+1}&=๐_{k+1}/m \\

๐ฑ_{k+1}&=๐ฑ_k+๐ฏ' ฮt \\

๐ฏ_{k+1}&=๐ฏ'+๐_{k+1} ฮt/2
\end{aligned}
\label{eq:24}
\tag{24}
$$

The main difference to Euler integration is that Verlet integration updates the velocity using the mean of current and previous accelerations. Verlet integration for the rotational motion work similarly.
$$
\begin{aligned}
ฯ'&=ฯ_k+ฮฑ_k ฮt/2 \\

ฮฑ_{k+1}&=M_{k+1}/I \\

ฯ_{k+1}&=ฯ_k+ฯ' ฮt \\

ฯ_{k+1}&=ฯ'+ฮฑ_{k+1} ฮt/2
\end{aligned}
\label{eq:25}
\tag{25}
$$

Also, we wrap the new orientation to pi $$ ฯ_{k+1} = w_ฯ(ฯ_{k+1}) $$

## Parameter and Attribute Values

adult | male | female | child | eldery | |
---|---|---|---|---|---|

torso ratio $k_{t}$ [%] | 0.5882 | 0.5926 | 0.5833 | 0.5714 | 0.6 |

shoulder ratio $k_{s}$ [%] | 0.3725 | 0.3704 | 0.375 | 0.3333 | 0.36 |

torso-to-shoulder ratio $k_{ts}$ [%] | 0.6275 | 0.6296 | 0.625 | 0.6667 | 0.64 |

average radius $r$ [m] | 0.255ยฑ0.035 | 0.27ยฑ0.02 | 0.24ยฑ0.02 | 0.21ยฑ0.015 | 0.25ยฑ0.02 |

average velocity $v$ [m/s] | 1.25ยฑ0.3 | 1.35ยฑ0.2 | 1.15ยฑ0.2 | 0.9ยฑ0.3 | 0.8ยฑ0.3 |

average mass $m$ [kg] | 73.5 | 80 | 67 | 57 | 70 |

Table: Anthropological data of human dimensions. ^{11}

Agent $iโA$ initial values and attributes are:

- $๐ฑ_i$: The center of mass is drawn from random uniform distribution inside a designated source without overlapping with other agents or obstacles.
- $๐ฏ_i$: The velocity is set zero vector or to an instance dependent value.
- $๐_i$: The force is set to zero vector.
- $\hat{๐}_i$: The desired direction is set to zero vector or to an instance dependent value.
- $v_i$: The desired velocity is set according to the anthropological data, for example, as the average walking speed of a human.
- $r_i$: The radius is set according to the anthropological data.
- $m_i$: The mass is set according to the anthropological data.

Orientable agent model initial values and attributes:

- $ฯ_i$: The orientation is set to same angle an initial direction.
- $ฯ_i$: The angular velocity is set to zero.
- $M_i$: The torque is set to zero.
- $ฯ_i^0$: The desired orientation set to same angle an initial direction.
- $ฯ_i^0$: The desired angular velocity is set to $2ฯ/3.$
- $I_i$: The rotational moment is set to value $4ฯ (m_i/m_0) (r_i/r_0)^2$ where $m_0=80$ and $r_0=0.27$, which means that agent with mass $m_0$ and radius $r_0$ has moment $4ฯ$ and other values are scaled according to the equation for moment of inertia.

Three-circle agent model attributes:

- $r_i^t=k_{t} r_i$: Torso radius is set according to the anthropological data.
- $r_i^s=k_{s} r_i$: Shoulder radius is set according to the anthropological data.
- $r_i^{ts}=k_{ts} r_i$: Torso-to-shoulder radius is set according to the anthropological data.

Force parameters:

- $t_{adj}=0.5$
- $ฮบ=4โ 10^4$
- $ฮผ=1.2โ 10^5$
- $ฮณ=500$
- $ฯ_ฮพ^2=0.1^2$
- $k=1.5$
- $ฯ_{soc}=3.0$

Torque parameters:

- $ฯ_{rot}=0.2$
- $ฯ_ฮถ^2=0.1$

Integrator parameters:

- $ฮt_{min}=0.001$
- $ฮt_{max}=0.01$

## Implementation

We have experimented with the algorithms in the crowddynamics library and verified that they operate correctly. However, based on the experiences of implementing the library, we learned that creating large, complex simulations is implausible without a *graphical user interface* (GUI) and *extendable program architecture*. For example, if we wish to create and edit the field, logic, and set initial parameter and attribute values graphically, we require these features. For this reason, we advise building the simulation on top of existing software, such as a *physics engine* or *game engine*. We have the following suggestion for the software features required of this software.

*Optional GUI*. Ability to simulate with or without the graphical user interface. For example, we may wish to develop the simulation on our local machine with graphics and then run multiple scenarios in parallel in the cloud without the graphics only to obtain numerical results.*IO*. Support for saving and loading simulation data and configurations from and to file. Essential for sharing and reproducing simulations.*Interactive graphics*. Display the simulation progress in real-time. Buttons for starting, stopping, and resuming the simulation. Menu options for choosing the simulation data and configurations files.*Field, parameter, and attribute editor*. Ability to manipulate simulation fields graphically. Create, edit, and remove obstacles, agents, sources, and targets. Set attribute values for agents and select their source positions and targets.*Logic editor*. Ability to create individual simulation logic blocks (e.g., forces, torque, integrator) and then connect them into the full simulation logic graphically. Allows changing individual parts of simulation logic for developing and testing new algorithms without having to break simulations that depend on different logic configurations.

The programming language choice depends on the chosen engine. Engines typically use C++. Also, installation and getting started should be made as easy as possible to attract users outside the research community.

## Conclusion

In reality, human movement can be very complicated. For example, the motivation to move to a specific direction can be affected by many factors, such as sensory inputs and internal motivations. It would be difficult and even undesirable to attempt to model all the aspects of human movement. Therefore, the models we have presented in this article focus on capturing only the essential aspects of the human movement. It requires making simplifying assumptions in the models.

*Newtonian model*. Numerical simulation of Newtonian mechanics can create artifacts and non-realistic behavior, which we should account for in the implementation. For these reasons, we use adaptive timestep and Verlet integration.*Planar motion*. In reality, human movement is 3-dimensional, but creating 3-dimensional simulation is significantly more challenging.*Agent models*. The circular agent model is sufficient for modeling low-density crowds, but in high-density crowds, it is unrealistic. The orientable, three-circle model improves this shortcoming by more accurately representing the human torso and allowing body rotation. However, it requires algorithms for rotational motion and improved collision detection. Also, desired body rotation is not a completely solved problem.

The results from simulations are stochastic. Therefore, we need to run multiple simulations with different initial configurations and analyze the resulting data using statistical methods to obtain meaningful information.

## Contribute

If you enjoyed or found benefit from this article, it would help me share it with other people who might be interested. If you have feedback, questions, or ideas related to the article, you can write to my GitHub Discussions forum.

***

*For more content, you can follow my YouTube channel and join my newsletter. Since creating content and open-source libraries take time and effort, consider supporting the effort by subscribing or giving a one-time donation.*

## References

Helbing, D., Farkas, I., & Vicsek, T. (2000). Simulating dynamical features of escape panic. Nature, 407(6803), 487โ490. https://doi.org/10.1038/35035023 ↩︎

Gibelli, L., & Bellomo, N. (2019). Crowd Dynamics, Volume 1: Theory, Models, and Safety Problems. Springer International Publishing. Retrieved from https://books.google.com.br/books?id=Pd2EDwAAQBAJ ↩︎

Langston, P. A., Masling, R., & Asmar, B. N. (2006). Crowd dynamics discrete element multi-circle model. Safety Science. https://doi.org/10.1016/j.ssci.2005.11.007 ↩︎

Korhonen, T., Hostikka, S., Heliรถvaara, S., & Ehtamo, H. (2008). FDS+ Evac: An Agent Based Fire Evacuation Model. Pedestrian and Evacuation Dynamics 2008, 109โ120. https://doi.org/10.1007/978-3-642-04504-2 ↩︎

Stรผvel, S. A. (2016). Dense Crowds of Virtual Humans. ↩︎

Hidalgo, R. C., Parisi, D. R., & Zuriguel, I. (2017). Simulating competitive egress of noncircular pedestrians. PHYSICAL REVIEW E, 95. https://doi.org/10.1103/PhysRevE.95.042319 ↩︎

Karamouzas, I., Skinner, B., & Guy, S. J. (2014). Universal power law governing pedestrian interactions. Physical Review Letters. https://doi.org/10.1103/PhysRevLett.113.238701 ↩︎

Karamouzas, I., Skinner, B., & Guy, S. J. (2014). Supplemental material for: Universal Power Law Governing Pedestrian Interactions Ioannis. Physical Review Letters. https://doi.org/10.1103/PhysRevLett.113.238701 ↩︎

Kretz, T., Groรe, A., Hengst, S., Kautzsch, L., Pohlmann, A., & Vortisch, P. (2011). Quickest Paths in Simulations of Pedestrians. Advances in Complex Systems, 14(5), 733โ759. https://doi.org/10.1142/S0219525911003281 ↩︎

Cristiani, E., & Peri, D. (2015). Handling obstacles in pedestrian simulations: Models and optimization. Retrieved from http://arxiv.org/abs/1512.08528 ↩︎

Korhonen, T., Hostikka, S., Technical, V. T. T., Ehtamo, H., Analysis, S., Box, T. P. O., โฆ The, I. (2008). FDS+Evac: Modelling Social Interactions in Fire Evacuation. https://doi.org/10.1007/978-3-642-04504-2_8 ↩︎