How to Implement ContinuousTime MultiAgent Crowd Simulation
Table of Contents
Introduction
In this article, we discuss the mathematics on how to implement a continuoustime multiagent crowd simulation. The article covers two agent models, circular, unorientable model, and threecircle, 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 continuoustime, Newtonian mechanics with an agentbased model, where each agent is a rigid body, in a twodimensional lattice. In the model, the agents are intelligent selfpropelled particles, whose movement is produced by nonphysical and physical forces. The nonphysical 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 continuumbased models which model the crowd as a continuous mass. The agentbased model has the advantage of being the most realistic model, but it is also computationally the most challenging model.
There are reallife 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 indepth the mathematics behind the agentbased model and how to implement the simulation.
Geometric Primitives
We define twodimensional 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,…,m1$ 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 agentbased 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 threecircle, 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 mediumdensity crowds. However, it is not a very realistic model of the crosssection 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
Threecircle Model
The threecircle model consists of three circles representing the shoulders and midsection. The threecircle model is a more realistic representation of the human crosssection. 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 threecircle agent has attributes as follows.
 $r_i^{t}$ torso radius
 $r_i^{s}$ shoulder radius
 $r_i^{ts}$ torsotoshoulder 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 threecircle agent model by deriving the variables, social force, and contact force between threecircle agents and threecircle 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{๐ฏ}$, skintoskin 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 skintoskin distance $h=\tilde{๐ฑ}(r_j+r_i).$
Between Circular Agent and Lineobstacle
Between a circular agent $i$ and static lineobstacle $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: Skintoskin 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 lineobstacle.
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 lineobstacles. 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 lineobstacles 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 mechanicsbased 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 pairwise 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 timetocollision $ฯ$ is obtained by linear extrapolation of the current trajectories of the agents and solving the root for skintoskin 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 skintoskin 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^2ac}}{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^2a 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 exittime 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 unitvector field, referred as a direction map, which points towards the boundary with the lowest exittime 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 exittime penalty $T(๐ฑ) = 0,, ๐ฑ โ ๐ฃ$ and $T(๐ฑ) โ โ,, ๐ฑ โ ๐'$. Let $O(๐ฑ)$ be the distance map to obstacles $๐$ with domain $ฮฉ$ and exittime 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 agentbased simulation.
Computing Interactions
In the simulation, we have two types of interactions, agentagent and agentobstacle interactions. Interactions are computationally the most challenging part of the simulation, and therefore it is worthwhile 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 timetocollision 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.
AgentAgent Interactions
To find neighboring agent pairs within distance $r$, we can use fixedradius 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:
 Agentagent 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.
AgentObstacle 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 secondorder 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 firstorder 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 
torsotoshoulder 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.
Threecircle 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$: Torsotoshoulder 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 realtime. 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 nonrealistic 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 3dimensional, but creating 3dimensional simulation is significantly more challenging.

Agent models. The circular agent model is sufficient for modeling lowdensity crowds, but in highdensity crowds, it is unrealistic. The orientable, threecircle 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 opensource libraries take time and effort, consider supporting the effort by subscribing or giving a onetime 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 multicircle 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/9783642045042 ↩︎

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/9783642045042_8 ↩︎