How to Implement Continuous-Time Multi-Agent Crowd Simulation

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 ¹.

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 ² 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 $𝐩\in {ℝ}^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,\dots ,{𝐩}_m),m\ge 2$ where each pair $({𝐩}_i,{𝐩}_{i+1})$ for $i=1,\dots ,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 $\Omega \subset {ℝ}^2$, a subset of the lattice which contains all the other objects.

• The set of agents ${\mathcal{A}}_i\subset \Omega$ for $i\in A$ where $A=1,\dots ,n$ are the indices. An agent is an impassable, dynamic object. We will denote two distinct agents using the indices as $i,j\in A$ where $i\ne j$.

• The set of obstacles $\mathcal{O}\subset \Omega$. An obstacle is an impassable, static object. We denote obstacles using the subscript $w\in W$, that is, ${\mathcal{O}}_w\subset \mathcal{O}$.

• The set of targets $\mathcal{T}\subset \Omega$. A target is a passable, static object.

• The set of sources $\mathcal{S}\subset \Omega$. 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\in A$ has the following attributes.

• ${𝐱}_i$ center of mass
• ${𝐯}_i$ velocity
• ${𝐟}_i$ force
• ${\widehat{𝐞}}_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:

• ${\phi }_i$ orientation
• ${\omega }_i$ angular velocity
• $M_i$ torque
• ${\phi }_i^0$ desired orientation
• ${\omega }_i^0$ desired angular velocity
• $I_i$ rotational moment

We constrained the values of the orientations $\phi$ and ${\phi }^0$ into interval $[-\pi ,\pi ].$

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 $\widetilde{𝐱}$, relative velocity $\widetilde{𝐯}$, skin-to-skin distance $h$ (can be negative), normal (unit) vector $\widehat{𝐧}$ and tangent (unit) vector $\widehat{𝐭}$. We use these variables for computing the interaction potentials. We use $R(\varphi )$ to define the rotation matrix in two dimensions where $\varphi \in [0°,360°]$ is the angle in degrees.

Between Two Distinct Circular Agents

Between two distinct circular agents $i,j\in A,i\ne j$ we have the variables as follows: Relative position $\widetilde{𝐱}={𝐱}_j-{𝐱}_i$, relative velocity $\widetilde{𝐯}={𝐯}_j-{𝐯}_i$, normal vector $\widehat{𝐧}=\widetilde{𝐱}/|\widetilde{𝐱}|$, tangent vector $\widehat{𝐭}=R(270°)\cdot \widehat{𝐧}$ and skin-to-skin distance $h=|\widetilde{𝐱}|-(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 $\widetilde{𝐩}={𝐩}_2-{𝐩}_1$, line length $l=|\widetilde{𝐩}|$, line tangent ${𝐭}_w=\widetilde{𝐩}/l$, line normal ${𝐧}_w=R(90°)\cdot {𝐭}_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\cdot ({𝐪}_1+{𝐪}_2).$

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

1. $l_t>l$: Center of mass ${𝐱}_i$ is on the left of the line perpendicular to $\widetilde{𝐩}$ intersecting ${𝐩}_1$.
2. $l_t<-l$: Center of mass ${𝐱}_i$ is on the right of the line perpendicular to $\widetilde{𝐩}$ intersecting ${𝐩}_2$.
3. 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 $\widehat{𝐭}=R(270°)\cdot \widehat{𝐧}$, and relative velocity $\widetilde{𝐯}={𝐯}_i-{𝐯}_w={𝐯}_i.$ We do not need the relative position $\widetilde{𝐱}$ between agent and line-obstacle.

Total Force

The total force on an agent $i\in 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 $\widehat{𝐞}$ in the characteristic time ${\tau }_{adj}$. We define it as $$\begin{array}{}\text{(2)}& {𝐟}_{adj}=\frac{m}{{\tau }_{adj}}(v\widehat{𝐞}-𝐯),\end{array}$$ where $m$ is the mass, $v$ is the desired velocity, $𝐯$ is the current velocity, and ${\tau }_{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 $$\begin{array}{}\text{(3)}& {𝐟}_{con}=\{\begin{array}{ll}h\kappa (\widetilde{𝐯}\cdot \widehat{𝐭})\widehat{𝐭}-h\mu \widehat{𝐧}+\gamma (\widetilde{𝐯}\cdot \widehat{𝐧})\widehat{𝐧}& h\le 0\\ 0,& h>0\end{array}\end{array}$$ where $\kappa$ is tangential friction coefficient, $\mu$ is the counter compressive coefficient, $\gamma$ is the damping coefficient, $\widehat{𝐧}$ is the normal unit vector and $\widehat{𝐭}$ is the normal tangent vector. It follows the original formula used by ¹ with the addition of the damping by ³.

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. $$\begin{array}{}\text{(4)}& \mathbf{\xi }=\xi \widehat{𝐮}(\phi )\end{array}$$

We draw the magnitude of the fluctuation force from a truncated normal distribution $\xi \sim N(0,{\sigma }_{\xi }^2)$ where ${\sigma }_{\xi }^2$ is the variance and the direction from a uniform distribution $\phi \sim U(0,2\pi )$. The normal distribution $N$ is truncated to three standard deviations. We define $\widehat{𝐮}(\phi )$ to be a unit vector to direction $\phi$.

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 $$\begin{array}{}\text{(5)}& E(\tau )=m\frac{k}{{\tau }^2}\exp (-\frac{\tau }{{\tau }_{soc}}),\end{array}$$ where $k$ is a social force scaling coefficient and ${\tau }_{soc}$ is the interaction time horizon. The time-to-collision $\tau$ 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 $\tau$ between these agents. The resulting $\tau$ is a function of the relative position $\widetilde{𝐱}$. 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 $\tau$ 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{array}{rl}h(\tau )& =|({𝐱}_i+\tau {𝐯}_i)-({𝐱}_j+\tau {𝐯}_j)|-(r_i+r_j)\\ & =|\widetilde{𝐱}+\tau \widetilde{𝐯}|-(r_i+r_j).\end{array}$$ We receive the quadratic equation in terms of $\tau$ by settings $h(\tau )=0$ $$\begin{array}{rl}|\widetilde{𝐱}+\tau \widetilde{𝐯}|-(r_i+r_j)& =0\\ |\widetilde{𝐱}+\tau \widetilde{𝐯}{|}^2& =(r_i+r_j{)}^2\\ |\widetilde{𝐯}{|}^2{\tau }^2+(\widetilde{𝐱}\cdot \widetilde{𝐯})\tau +|\widetilde{𝐱}{|}^2-(r_i+r_j{)}^2& =0.\end{array}$$ We solve the quadratic equation $$\begin{array}{}\text{(7)}& \tau =\frac{b-\sqrt{b^2-ac}}{a},\end{array}$$ where $a=|\widetilde{𝐯}{|}^2$, $b=-\widetilde{𝐱}\cdot \widetilde{𝐯}$, and $c=|\widetilde{𝐱}{|}^2-(r_i+r_j{)}^2$. Now we can solve the gradient $$\begin{array}{}\text{(8)}& {\nabla }_{\widetilde{𝐱}}\tau =\frac{1}{a}(\widetilde{𝐯}-\frac{a\widetilde{𝐱}+b\widetilde{𝐯}}{\sqrt{b^2-ac}}).\end{array}$$ By substituting the gradient to the social force, we can use it in the simulation.

Total Torque

The total torque on agent $i\in 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 ${\phi }_0$ in the characteristic time ${\tau }_{rot}$. We define it as $$\begin{array}{}\text{(10)}& M_{adj}=\frac{I}{{\tau }_{rot}}(\frac{w_{\pi }({\phi }_0-\phi )}{\pi }{\omega }_0-\omega )\end{array}$$

The function $w_{\pi }$ wraps an arbitrary angle $\phi$ in radians to the interval $(-\pi ,\pi ]$. The function $w_{\pi }$ is defined as $$\begin{array}{}\text{(11)}& \begin{array}{rl}{w}_{\pi }(\phi )& =w_{\pi }^{\prime }(\phi \mathrm{mod}2\pi )\\ w_{\pi }^{\prime }(\hat{\phi })& =\{\begin{array}{ll}\hat{\phi }& \hat{\phi }\le \pi \\ \hat{\phi }-2\pi & \hat{\phi }>\pi \end{array}.\end{array}\end{array}$$

We can set the desired orientation ${\phi }_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 $$\begin{array}{}\text{(12)}& M_{con}=𝐫\times {𝐟}_{con},\end{array}$$ 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 $$\begin{array}{}\text{(13)}& 𝐮\times 𝐯=u_0{v}_1-u_1{v}_0.\end{array}$$ 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 $$\begin{array}{}\text{(14)}& \eta =I\zeta ,\end{array}$$ where the magnitude is drawn from a truncated normal distribution $\zeta \sim N(0,{\sigma }_{\zeta }^2)$ with variance ${\sigma }_{\zeta }^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 $\partial \Omega$. The eikonal equation is defined as $$\begin{array}{}\text{(15)}& \begin{array}{rl}|S(𝐱)|& =\frac{1}{f(𝐱)},𝐱\in \Omega ,\\ S(𝐱)& =q(𝐱),𝐱\in \partial \Omega \end{array}\end{array}$$ where $f:\Omega \to (0,+\infty )$ is the speed of travel and $q:\partial \Omega \to [0,+\infty )$ 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 $$\begin{array}{}\text{(16)}& D(S(𝐱))=-\frac{\nabla S(𝐱)}{|\nabla S(𝐱)|}.\end{array}$$

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 ${𝓞}^{\prime }$ 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 $\Omega$ and exit-time penalty $T(𝐱)=0,,𝐱\in 𝓣$ and $T(𝐱)\to \infty ,,𝐱\in {𝓞}^{\prime }$. Let $O(𝐱)$ be the distance map to obstacles $𝓞$ with domain $\Omega$ and exit-time penalty $O(𝐱)=0,,𝐱\in 𝓞$. The desired direction is then defined as $$\begin{array}{}\text{(17)}& \widehat{𝐞}(𝐱)=\lambda (O(𝐱))D(T(𝐱))+(1-\lambda (O(𝐱)))D(O(𝐱)),\end{array}$$ where $\lambda :{ℝ}^{+}\to [0,1]$ is decreasing function, $\frac{d}{dx}\lambda (x)\le 0$, such that $\lambda (0)=1$ and $\underset{x\to \infty }{lim}\lambda (x)=0$.

We will set the avoidance radius $r$ sufficiently larger than $\underset{i\in A}{max}{r}_i$, which ensures that agents tend away from the obstacles.

Exponential λ-function

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

1. $\lambda (0)=a=1$
2. $\lambda (r)=b^{cr}=ϵ$ gives us $c=\frac{{\log }_bϵ}{r}$ therefore $\lambda (x)={ϵ}^{x/r}$
3. $\underset{x\to \infty }{lim}\lambda (x)=0$ implies $ϵ<1$
4. Lastly, we verify that the function is decreasing, that is, $\frac{d}{dx}\lambda (x)=\frac{\log ϵ}{r}{ϵ}^{x/r}<0.$ Note that $\log ϵ<0$ because $0<ϵ<1.$

The resulting function $\lambda$ is defined as $$\begin{array}{}\text{(18)}& \lambda (x)={ϵ}^{x/r},\end{array}$$ 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 $\lambda (r)=0$. Then we have $$\begin{array}{}\text{(19)}& \lambda (x)=\{\begin{array}{ll}-x/r+1& x\le r\\ 0& x>r\end{array}.\end{array}$$ 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 $\text{(???)}$ and $\text{(???)}$ 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\ge 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\ge 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:

1. Agent-agent interactions are local.
2. 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 $$\begin{array}{}\text{(20)}& m\frac{d^2}{d{t}^2}𝐱(t)=𝐟(t).\end{array}$$

The rotational motion is defined by the equation $$\begin{array}{}\text{(21)}& I\frac{d^2}{d{t}^2}\phi (t)=M(t).\end{array}$$

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 $\Delta t$.

Adaptive Timestep

Adaptive timestep is selected from interval $[\Delta t_{min},\Delta t_{max}]$ such that we limit the maximum distance that an agent can move in each iteration. First, we define $$\begin{array}{}\text{(22)}& \Delta t_{mid}=\Delta t_{max}\frac{\underset{i\in A}{max}{v}_i}{\underset{i\in A}{max}|{𝐯}_i|}\end{array}$$ where $|{𝐯}_i|$ is agent’s current velocity and $v_i$ is agent’s target velocity. Then, the adaptive timestep is defined $$\begin{array}{}\text{(23)}& \Delta t=\{\begin{array}{ll}\Delta t_{min}& \Delta t_{mid}<\Delta t_{min}\\ \Delta t_{mid}& \mathrm{otherwise}\\ \Delta t_{max}& \Delta t_{mid}>\Delta t_{max}\end{array}.\end{array}$$

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,\dots$ as follows

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{𝐯}^{\prime }& ={𝐯}_k+{𝐚}_k\Delta t/2\\ {𝐚}_{k+1}& ={𝐟}_{k+1}/m\\ {𝐱}_{k+1}& ={𝐱}_k+{𝐯}^{\prime }\Delta t\\ {𝐯}_{k+1}& ={𝐯}^{\prime }+{𝐚}_{k+1}\Delta t/2\end{array}\end{array}$$

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{array}{}\text{(25)}& \begin{array}{rl}{\omega }^{\prime }& ={\omega }_k+{\alpha }_k\Delta t/2\\ {\alpha }_{k+1}& =M_{k+1}/I\\ {\phi }_{k+1}& ={\phi }_k+{\omega }^{\prime }\Delta t\\ {\omega }_{k+1}& ={\omega }^{\prime }+{\alpha }_{k+1}\Delta t/2\end{array}\end{array}$$

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

Parameter and Attribute Values

Table: Anthropological data of human dimensions. ¹¹

Agent $i\in 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.
• ${\widehat{𝐞}}_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:

• ${\phi }_i$: The orientation is set to same angle an initial direction.
• ${\omega }_i$: The angular velocity is set to zero.
• $M_i$: The torque is set to zero.
• ${\phi }_i^0$: The desired orientation set to same angle an initial direction.
• ${\omega }_i^0$: The desired angular velocity is set to $2\pi /3.$
• $I_i$: The rotational moment is set to value $4\pi (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\pi$ 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$
• $\kappa =4\cdot {10}^4$
• $\mu =1.2\cdot {10}^5$
• $\gamma =500$
• ${\sigma }_{\xi }^2={0.1}^2$
• $k=1.5$
• ${\tau }_{soc}=3.0$

Torque parameters:

• ${\tau }_{rot}=0.2$
• ${\sigma }_{\zeta }^2=0.1$

Integrator parameters:

• $\Delta t_{min}=0.001$
• $\Delta 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 contact me via email. For more content, you can follow me on YouTube or join my newsletter. Creating content takes time and effort, so consider supporting me with a one-time donation.

References