A Kinetic Theory Approach to Model Pedestrian Dynamics in Bounded Domains with Obstacles

We consider a kinetic theory approach to model the evacuation of a crowd from bounded domains. The interactions of a person with other pedestrians and the environment, which includes walls, exits, and obstacles, are modeled by using tools of game theory and are transferred to the crowd dynamics. The model allows to weight between two competing behaviors: the search for less congested areas and the tendency to follow the stream unconsciously in a panic situation. For the numerical approximation of the solution to our model, we apply an operator splitting scheme which breaks the problem into two pure advection problems and a problem involving the interactions. We compare our numerical results against the data reported in a recent empirical study on evacuation from a room with two exits. For medium and medium-to-large groups of people we achieve good agreement between the computed average people density and flow rate and the respective measured quantities. Through a series of numerical tests we also show that our approach is capable of handling evacuation from a room with one or more exits with variable size, with and without obstacles, and can reproduce lane formation in bidirectional flow in a corridor.


Introduction
The complex dynamical behavior of pedestrian crowds has fascinated researchers from various scientific fields since the early 1950's. Academic studies started with empirical observations and continued with the development of models in the field of applied physics and mathematics. The simulation of pedestrian flow has attracted increasing research attention in recent years since a reliable simulation model for pedestrian flow may greatly benefit engineers in mass transportation management, and designers in urban planning and architecture.
A very large variety of models have been developed over the years. The different mathematical models can be divided into three main categories depending on the scale of observation [13]. A first approach corresponds to the macroscopic description: evolution equations are derived for mass density and linear momentum, which are regarded as macroscopic observables of pedestrian flow. See, e.g., [27,36]. Such an approach is suitable for high density, large-scale systems, which are not the focus of our work.
A second approach looks at the problem at the microscopic level. Microscopic models can be further divided into models which are grid-based or grid-free. Cellular Automata [16,15,14,20,30] models belong to the first category. They describe pedestrian flow in space-time by assigning discrete states to a grid of space-cells. These cells can be occupied by a pedestrian or be empty. Thus, the movement of pedestrians in space is done by passing them from cell to cell (discrete space) in discrete time. Grid-free methods use Newtonian mechanics to interpret pedestrian movement as the physical interaction between the people and the environment, i.e. the action of other people and the environment on a given pedestrian is modeled with forces. These microscopic models, also called force-based models, are one of the most popular modeling paradigms of continuous models because they describe the movement of pedestrians qualitatively well. See, e.g., [24,25,26,38,32,18,42,31] and references therein. Collective phenomena, like unidirectional or bidirectional flow in a corridor, lane formation, oscillations at bottlenecks, the faster-is-slower effect, and emergency evacuation from buildings, are well reproduced. Force-based models can be extended into agent-based models by incorporating individual features. See, e.g., [19,17,3,2,4] and references therein. Agent-based models allow for flexibility, extensibility, and capability to realize heterogeneity in crowd dynamics. Both force-based and agent-based models may introduce artifacts due to the force representation of human behavior, leading to unrealistic backward movement or oscillating trajectories. These artifacts can be reduced by incorporating extra rules and/or elaborate calibrations, at the price of an increased computational cost.
The scale of observation for the third approach is between the previous two. Introduced in [12] and further developed in [10,6,1,7,8,5,9], this approach derives a Boltzmann-type evolution equation for the statistical distribution function of the position and velocity of the pedestrians, in a framework close to that of the kinetic theory of gases. See also [11] for a literature review on this approach. The model proposed in [12,10,6] is valid in unbounded domains and with a homogeneous distribution of walking ability for the pedestrians, while the extension to bounded domains is presented in [1] and further explored in [7,8,9]. In [7], more general features of behavioralsocial dynamics are taken into account. In [8], Monte Carlo simulations are introduced to study pedestrians behavior in complex scenarios. The methodology in [8] is validated by comparing the computed fundamental density-velocity diagrams with empirically observed ones and by checking that well known emerging properties are reproduced. A kinetic theory approach for modeling pedestrian dynamics in presence of social phenomena, such as the propagation of stress conditions, is presented in [9]. Finally, we refer to [5] for a thorough description of how kinetic theory and evolutionary game theory can be used to understand the dynamics of living systems.
In this work, we consider the model proposed in [1]. We first validate it against experimental data and then extend it to bounded domains with obstacles. It is worth noticing that most of the models and methods in the references cited so far have been shown to reproduce phenomena of pedestrian movement qualitatively through analysis and/or numerical simulations. However, before using a model to predict quantitative results like, e.g., the total evacuation time, the mathematical models have to be validated and the numerical methods have to be verified [21]. In the context of pedestrian dynamics, this is still difficult due to a lack of reliable experimental data. In addition, the few available datasets show large differences [35,33,41]. In this paper, we compare our numerical results against the data reported in a recent empirical study [39]. We have selected this study because it deals with egressing from a facility and thus it is the most directly related to our focus. With the model under consideration, for medium and medium-to-large groups of people we manage to achieve good agreement between the computed average people density and flow rate and the respective measured quantities. Finally, we mention that in order to make the models more reliable, evolutionary adjustment of the parameters and data assimilations have also been proposed in [28,40], respectively.
The strategy we propose to handle obstacles within the computational domain makes use of an effective obstacle area, which is an enlarged area that encloses the real obstacle, and a model parameter used to describe the quality of the environment. Thanks to those two ingredients, we can successfully exclude from the walkable area square and rectangular obstacles. We test our strategy in a square room that contains one or two obstacles and has one exit. In addition, through a series of numerical tests we show that our approach is capable of handling evacuation from a room with one or more exits with variable size, with and without obstacles, and can reproduce lane formation in bidirectional flow in a corridor.
The paper is organized as follows. Sec. 2 introduces the representation of the system and the modeling of interactions with pedestrians and with the envornment. In Sec. 3, we apply the Lie splitting algorithm to the model described in Sec. 2. Numerical results are presented in In Sec. 4 and conclusions are drawn in Sec. 5.

Mathematical model
The model we consider is based on the model proposed in [1]. Let Ω ⊂ R 2 denote a bounded domain. We assume that the boundary ∂Ω includes an exit E which could be the finite union of disjoint sets, and walls W . Here, E ∪ W = ∂Ω and E ∩ W = ∅. Let x = (x, y) ∈ Ω denote position and v = v(cos θ, sin θ) ∈ Ω v denote velocity, where v is the velocity modulus, θ is the velocity direction, and Ω v ⊂ R 2 is the velocity domain. For a system composed by a large number of pedestrians distributed inside Ω, the distribution function is given by Under suitable integrability conditions, f (t, x, v)dxdv represents the number of individuals who, at time t, are located in the infinitesimal rectangle [x, x + dx] × [y, y + dy] and have a velocity belonging to [v, v + dv] × [θ, θ + dθ]. Since we use polar coordinates for the velocity, we can write the distribution function as f = f (t, x, v, θ).
Following [1], we assume that variable θ is discrete. This assumption is motivated by the granularity of pedestrian dynamics when the crowd size is not enough to justify the continuity of the distribution function over the variable θ. For simplicity, we assume θ can take values in the set: where N d is the maxim number of possible directions. As for the velocity magnitude v, we model it as a continuous deterministic variable which evolves in time and space according to macroscopic effects determined by the overall dynamics. In fact, experimental studies show that in practical situations the speed of pedestrians depends mainly on the level of congestion around them. Due to the deterministic nature of the variable v, the kinetic type representation is given by the reduced distribution function where f i (t, x) = f (t, x, θ i ) represents the active particles that, at time t and position x, move with direction θ i . In eq. (1), δ denotes the Dirac delta function. Let us introduce some reference quantities that will be use to make the variable dimensionless. We define: -D: the largest distance a pedestrian can cover in domain Ω; -V M : the highest velocity modulus a pedestrian can reach in low density and optimal environmental conditions; -T : a reference time given by D/V M , -ρ M : the maximal admissible number of pedestrians per unit area. The dimensionless variables are then: positionx = x/D, timet = t/T , velocity modulusv = v/V M and distribution functionf = f /ρ M . From now on, all the variables will be dimensionless and hats will be omitted to simplify notation.
Due to the normalization of f , and of each f i , the dimensionless local density is obtained by summing the distribution functions over the set of directions: As mentioned above, we assume that pedestrians adjust their speed depending on the the level of congestion around them. This means that the velocity modulus depends formally on the local density, i.e. v = v[ρ](t, x), where square brackets are used to denote that v depends on ρ in a functional way. For instance, v can depend on ρ and on its gradient. A parameter α ∈ [0, 1] is introduced to represent the quality of the domain where pedestrians move: α = 0 corresponds to the worst quality which forces pedestrians to slow down or stop, while α = 1 corresponds to the best quality, which allows pedestrians to walk at the desired speed. We assume that the maximum dimensionless speed v M a pedestrian can reach depends linearly on the quality of the environment. For simplicity, we take v M = α. Let ρ c be a critical density value such that for ρ < ρ c we have free flow regime (i.e., low density condition), while for ρ > ρ c we have a slowdown zone (i.e., high density condition). We set ρ c = α/5. Note that this choice is compatible with the experimentally measured values of ρ c reported in [34]. Next, we set the velocity magnitude v equal to v M in the free flow regime and equal to a heuristic third-order polynomial in the slowdown zone: where a i is constant for i = 0, 1, 2, 3. To set the value of these constants, we impose the following conditions: v(ρ c ) = v M , ∂ ρ v(ρ c ) = 0, v(1) = 0 and ∂ ρ v(1) = 0. This leads to: (4) Fig. 1 (a) reports v as a function of ρ for α = 0.4, 0.7, 1.

Modeling interactions
Each pedestrian is modeled as a particle. Interactions involve three types of particles: test particles with state (x, θ i ): they are representative of the whole system; candidate particles with state (x, θ h ): they can reach in probability the state of the test particles after individual-based interactions with the environment or with field particles; field particles with state (x, θ k ): their presence triggers the interactions of the candidate particles. The process through which a pedestrian decides the direction to take is the results of several factors. We take into account four factors: (F1) The goal to reach the exit.
Given a candidate particle at the point x, we define its distance to the exit as and we consider the unit vector u E (x), pointing from x to the exit. See Fig. 1 (b). (F2) The desire to avoid the collision with walls.
Given a candidate particle at the point x moving with direction θ h , we define the distance  d W (x, θ h ) from the particle to a wall at a point x W (x, θ h ) where the particle is expected to collide with the wall. The unit tangent vector u W (x, θ h ) to ∂Ω at x W points to the direction of the the exit. Vector u W is used to avoid a collision with the walls. See Fig. 1 (b). (F3) The tendency to look for less congested area.
A candidate particle (x, θ h ) may decide to change direction in order to avoid congested areas. This is achieved with the direction that gives the minimal directional derivative of the density at the point x. We denote such direction by unit vector u C (θ h , ρ). (F4) The tendency to follow the stream.
A candidate particle modifies its state, in probability, into that of the test particle due to interactions with field particles, while the test particle loses its state as a result of these interactions. A candidate particle h interacting with a field particle k may decide to follow it and thus adopt its direction, denoted with unit vector u F = (cos θ k , sin θ k ).
Factors (F1) and (F2) are related to geometric aspects of the domain, while factors (F3) and (F4) consider that people's behavior is strongly affected by surrounding crowd. Note that the effects related to factors (F3) and (F4) compete with each other: (F4) is dominant in a panic situation, while (F3) characterizes rational behavior. As a weight between (F3) and (F4), we introduce parameter ε ∈ [0, 1]: ε = 0 corresponds to the situation in which only the research of less congested areas is considered (rational behavior), while ε = 1 corresponds to the situation in which only the tendency to follow the stream is taken into account (panic behavior).

Interaction with the bounding walls
The interaction with the bounding walls is modeled with two terms: -µ[ρ]: the interaction rate models the frequency of interactions between candidate particles and the boundary of the domain. If the local density is getting lower, it is easier for pedestrians to see the walls and doors. Thus, we set the transition probability gives the probability that a candidate particle h, i.e. with direction θ h , adjusts its direction into that of the test particle θ i due to the presence of the walls and/or an exit. The following constraint for A h (i) has to be satisfied: We assume that particles change direction, in probability, only to an adjacent clockwise or counterclockwise direction in the discrete set I θ . This means a candidate particle h may end up into the states h − 1, h + 1 or remain in the state h. In the case h = 1, we set ..,N d forms the so-called table of games that models the game played by active particles interacting with the bounding walls.
To take into account factors (F1) and (F2), we define the vector Here θ G is the geometrical preferred direction, which is the ideal direction that a pedestrian should take in order to reach the exit and avoid the walls in an optimal way. Notice that the closer a pedestrian is to an exit (resp., a wall), the more direction u E (resp., u W ) weights. A candidate particle h will update its direction by choosing the angle closest to θ G among the three allowed angles θ h−1 , θ h and θ h+1 . The transition probability is given by: where s = arg min In (6), δ denotes the Kronecker delta function. Coefficient β h , proportional to parameter α, is defined by: where ∆θ = 2π/N d . The role of β h is to allow for a transition to θ h−1 or θ h+1 even in the case that the geometrical preferred direction θ G is closer to θ h . Such a transition is more likely to occur the more distant θ h and θ G are. Notice that if θ G = θ h , then β h = 0 and A h (h) = 1, meaning that a pedestrian keeps the same direction (in the absence of interactions other than with the environment) with probability 1.

Interaction with obstacles
The strategy reported in the previous section to avoid collisions with the walls works well when the pedestrian is sufficiently far from the walls. If pedestrians get too close to the bounding walls, and in particular if they are close to an exit, the definition of u G in (5) does not prevent collisions with the walls. Thus, obstacles within the domain Ω cannot be avoided just by adjusting u W . In this section, we report an effective strategy to handle obstacles. Three ingredients are needed to exclude the real obstacle area from the walkable domain: 1. An effective area: an enlarged area that encloses the real obstacle. 2. A definition of u W to account for the effective area. 3. A setting of the parameter α in the effective area depending on the shape of the obstacle. The effective area is necessary especially if the obstacle is close to an exit: it allows to define u W with respect to a larger area than the one occupied by the obstacle to achieve the goal of having no pedestrian walking on the real obstacle area. In the numerical results reported in Sec. 4.3, we used an effective area that is four times bigger than the real obstacle area.
Since some pedestrians will walk on part of the effective area, one needs to set parameter α. By setting α = 1 (i.e., best quality of the environment) in the effective area, pedestrians can move with the maximal velocity modulus as they approach the obstacle and thus they quickly adapt to the effective area through u W . However, some pedestrians will walk close to the top, bottom, and rear (with respect to the pedestrian motion) boundary of the effective area. Thus, the real obstacle is located at the front of the effective area. From the numerical results reported in Sec. 4.3, we also see that the shape of the obstacle is square. By setting α = 0 (i.e., worst quality of the environment) in the effective area, pedestrians are forced to slow down at the front part of the effective area. The slow down leads to higher densities in the front part of the effective area, therefore direction u W competes with direction u C . As a result some pedestrians walk on the front part of the effective area. However, as the congestion decreases pedestrains avoid the rear part of the effective area. From the numerical results shown in Sec. 4.3, we see that the shape of the obstacle for α = 0 in the effective area is slender.

Interactions between pedestrians
The interaction with other pedestrians is modeled with two terms: -η[ρ]: the interaction rate defines the number of binary encounters per unit time. If the local density increases, then the interaction rate also increases. For simplicity, we take η[ρ] = ρ.
Notice that unlike the case of classical gas dynamics, this rate is not related to the relative particle velocity.
: the transition probability gives the probability that a candidate particle h modifies its direction θ h into that of the test particle i, i.e. θ i , due to the research of less congested areas and the interaction with a field particle k that moves with direction θ k . The following constrain for B hk (i) has to be satisfied: where again the square brackets denote the dependence on the density ρ. The game consists in choosing the less congested direction among the three admissible ones. This direction can be computed for a candidate pedestrian h situated at x, by taking where ∂ j ρ denotes the directional derivative of ρ in the direction given by angle θ j . We have u C (θ h , ρ) = (cos θ C , sin θ C ). As for the tendency to follow the crowd, we set u F = (cos θ k , sin θ k ). This means that a candidate particle follows the direction of a field particle.
To take into account (F3) and (F4), we define the vector where the subscript P stands for pedestrians. Direction θ P is the interaction-based perferred direction, obtained as a weighted combination between the trendency to follow the stream and the tendency to avoid crowded zones. The transition probability is given by: where r and β hk are defined by: We recall that d(·, ·) is defined in (7).

Equation of balance
The derivation of the mathematical model can be obtained by a suitable balance of particles in an elementary volume of the space of microscopic states, considering the net flow into such volume due to transport and interactions [1]. Taking into account factors (F1)-(F4), we obtain: represents the net balance of particles that move with direction θ i due to interactions. As explained in the previous subsection, we consider both the interaction with the environment and with the surrounding people. Thus, we can write J i as J i = J i G + J i P , where J i G is an interaction between candidate particles and the environment and J i P is an interaction between candidate and field particles. Eq. (8) is completed with eq. (2) for the density and eq. (3),(4) for the velocity. In the next section, we will discuss a numerical method for the solution of problem (2),(3),(4),(8).

Numerical method
The approach we consider is based on a splitting method that decouples the treatment the transport term and the interaction term in eq. (8). As usual with splitting methods, the idea is to split the model into a set of subproblems that are easier to solve and for which practical algorithms are readily available. The numerical method is then completed by picking an appropriate numerical scheme for each subproblem. Among the available operator-splitting methods, we chose the Lie splitting scheme because it provides a good compromise between accuracy and robustness, as shown in [22].

The Lie operator-splitting scheme
Although the Lie splitting scheme is quite well-known, it may be useful to present briefly this scheme before applying it to the solution of problem (2),(3),(4), (8).
Let us consider a first-order system in time: where A is an operator from a Hilbert space into itself. Operator A is then split, in a non-trivial decomposition, as The Lie scheme consists of the following. Let ∆t > 0 be a time discretization step for the time interval [0, T ]. Denote t k = k∆t, with k = 0, . . . , N t and let φ k be an approximation of φ(t k ). Set φ 0 = φ 0 . For n ≥ 0, compute φ k+1 by solving and then set φ k+i/I = φ i (t k+1 ), for i = 1, . . . .I. This method is first-order accurate in time. More precisely, if (9) is defined on a finitedimensional space, and if the operators A i are smooth enough, then φ(t k ) − φ k = O(∆t) [22].
In the next section, we will apply Lie splitting to problem (8). The whole problem will be split into three subproblems: 1. A pure advection problem in the x direction. 2. A pure advection problem in the y direction. 3. A problem involving the interaction with the environment and other pedestrians.

Lie scheme applied to problem (8)
Let us apply the Lie operator-splitting scheme described in the previous section to problem (8).
In order to simplify notation, let us set Let M be a positive integer (≥ 3, in practice). We associate with M a time discretization step τ = (t f − t 0 )/M and set t m = t 0 + mτ . Next, we proceed with the space and time discretization of each subproblem in Sec. 3.2.
Step 1 Let φ 0 = f i,k . Problem (11) can be rewritten as The finite difference method we use produces an approximation Φ m p,q ∈ R of the cell average Φ m p,q ≈ 1 ∆x ∆y where m = 1, . . . , M , 1 ≤ p ≤ N x − 1 and 1 ≤ q ≤ N y − 1. Given an initial condition φ 0 , function φ m will be approximated by Φ m with Φ m The Lax-Friedrichs method for problem (14) can be written in conservative form as follows: Step 2 Let φ 0 = f i,k+ 1 3 . Problem (12) can be rewritten as Similarly to step 1, we use the conservative Lax-Friedrichs scheme: Step 3 Let J = J i and φ 0 = f i,k+ 2 3 . Problem (13) can be rewritten as For the approximation of the above problem, we use the Forward Euler scheme: where F m is the approximation of the reduced distribution function (1) at time t m .

Numerical results
We present a series of numerical results to showcase the features of our model. We start from the simulation of evacuation from a room with one exit in Sec. 4.1. We use this first test to validate our implementation of the model presented in Sec. 2 against the results in [1]. In order to further validate our software, in Sec. 4.2 we compare our numerical results with the experimental data reported in [39], where the authors study the evacuation from a room with two exits. Our successfully validated code is then used to study evacuation from a room with obstacles in Sec. 4.3 and lane formation in a corridor in Sec. 4.4. For all the simulations, we consider eight different velocity directions N d = 8 in the discrete set:

Evacuation from a room with one exit
This first test case is taken from [1]. The computational domain encloses a square room with side 10 m with an exit door located in the middle of the right side. The exit size is 2.6 m. The computational domain is larger than the room itself to follow the motion of the pedestrian also once they have left the room. We aim at simulating the evacuation of 46 people located inside the room and initially distributed into two equal-area circular clusters. See Fig. 2, top left panel. The two groups are initially moving against the each other with opposite initial directions θ 3 and θ 7 .
In order to work with dimensionless quantities as described in Sec. 2, we define the following reference quantities: D = 10  Figure 2: Evacuation process of 46 pedestrians grouped into two clusters with opposite initial directions θ 3 and θ 7 using the medium mesh and ∆t medium for time t = 0, 1.5, 3, 6, 10.5, 13.5 s. The color refers to density. Fig. 3 (a) reports the number of pedestrians left in the room computed with the meshes and time steps mentioned above. In all the cases the total evacuation time is around 18 s, which agrees well with the results reported in [1]. However, we see that the evacuation dynamics varies when the mesh and time step change. In fact, from Fig. 3 (a) one can observe that as the time step gets smaller with a given mesh people walk slightly faster, while as the mesh gets finer with a given time step pedestrians walk slightly slower. When using operator splitting methods, it is advised to reduce the time step as the mesh is refined (see, e.g. [22]). So, for ease of comparison in Fig. 3 (b) we report only the results computed with coarse mesh and ∆t large , medium mesh and ∆t medium , fine mesh and ∆t small . We can see very good agreement for those three curves. Next, we consider the same room as in the previous test but we vary the exit size. We locate the exit symmetrically with respect to the room centerline and let the exit size vary from 1.5 m to 4 m. We consider the medium meshes and ∆t medium mentioned before since Fig. 3 show that is an appropriate choice for the problem under consideration. All the other model and discretization parameters are set like in the previous test. Fig. 4 shows the total evacuation time as a function of the exit size. First, we notice that our results are in very good agreement with the results reported in [1]. As expected, the total evacuation time decreases with the exit size, but that once the exit is large enough for the crowd contained in the room the evacuation time does not change significantly if the exit is further enlarged.
As last test for the room with one exit, we computed the total evacuation time from a room with variable exit size for different values of ε. The results did not show any visible difference from the ones displayed in Fig. 4, which corresponds to ε = 0.4, and thus are not reported here. This test case is probably too simple to see any difference in the crowd dynamics when (F3), i.e. the tendency avoid the crowd, is dominant over (F4), i.e. the tendency to follow the crowd, and viceversa.

Evacuation from a room with two exits
After validating our software against the numerical results in [1], we proceed with the validation against the experimental data reported in [39]. The computational domain corresponds to the setting studied in [39]: a room of side 10 m with different sized exits placed on on the right side: exit 1 with size 0.7 m and exit 2 with size 1.1 m. The distance between the two exits is 3 m. See  [39] and initial density and direction (i.e., θ 1 ) for the experiment with 138 pedestrians.
The data in [39] refer to ten trials: 2 trials with 18 pedestrians, 6 trials with 40 pedestrians, and 2 trials with 138 pedestrians. In the experiments with 18 and 40 pedestrians (and so in the simulation as well), the group is initially positioned in a square located symmetrically with respect to the horizontal axis of the room and towards the back of the room with linearly increasing people density from the front to the back. In the experiments with 138 pedestrians, pedestrians are initially distributed as follows: a first group of 90 people is positioned in a rectangle in the middle of the room, while the remaining 48 people are positioned in a square behind the first group. Thus, we adopt the same configuration in the simulations. We impose constant density in the rectangle and prescribe a linearly increasing density from the front to the back of the square. See Fig. 5. In all the cases, the pedestrians are given initial direction θ 1 . For all the simulations, we use the medium mesh and ∆t small considered in Sec. 4.1.
To compute the mean density D V and mean flow rate F V we use the Voronoi method [37]: where V V is the mean velocity modulus over the area ω and E is the exit width. We choose two 4 m 2 areas in front of the exits as ω. The mean density and mean flow rate computed from our simulations are shown in Fig. 6 (a) and (b). Fig. 6 (c) and (d) report the measured mean density and mean flow rate form the experiments in [39]. We see very good agreement between computed and measured quantities for the 40 and 138 people cases, while there is no good agreement for the 18 people case. This is to be expected since the kinetic approach is not meant to simulate the dynamics of a small number of pedestrians.
See Fig. 7 (a). We repeat the test with 138 pedestrians for all the above velocity moduli. The corresponding number of pedestrians in the room versus times is shown in Fig. 7 (b). For a given density ρ in the slowdown zone, we have v(ρ) > v purple (ρ) > v orange (ρ) > v blue (ρ). Thus, it is not surprising the the total evacuation time increases as we pass from the original definition of velocity modulus in (3) to v blue through v purple and v orange . To further illustrate the difference in the evacuation process when the above velocity moduli are used, in Fig. 8 we display the density and velocity modulus (with selected velocity vectors) at time t = 15 when using the three velocity moduli. Next, we investigate the relationship between evacuation time and the width ratio of the two exits. We fix the size and position of exit 1, and the position of the center of exit 2, while the size of exit 2 varies. See Fig. 9 (a) for the widths of exit 2 under consideration and corresponding width ratios. We consider two scenarios: the group of 40 people with velocity modulus (3) and the group of 138 people with the blue velocity modulus (16). In both cases, all other model and discretization parameters are set like for the results reported in Fig. 6. Fig. 9 (b) shows the total evacuation time versus the exit width ratio for both scenarios. As expected, when the ratio of exit 2 width/exit 1 width increases the total evacuation time decreases.

Evacuation from a room with obstacles
In the model introduced in Sec. 2.1, parameter α represents the quality of the environment, which influences also the maximal dimensionless velocity modulus (3) a pedestrian can reach. In theory, parameter α = 0 forces pedestrians to stop, while the value α = 1 contributes to keep the maximal velocity modulus. However, in practice this parameter alone is not sufficient to model obstacles within the domain. To effectively model the obstacles, we use the strategy reported in Sec. 2.1.2.
We consider a square room of side 10 m with a 2.6 m wide exit located on the right wall as in Sec. 4.1, with the following obstacle configuration: 1. One obstacle close to the exit, i.e. in the middle of the right wall, with effective area depicted in Fig. 10 (a). 2. Two obstacles close to the right wall, place symmetrically with respect to the exit. See Fig. 10 (b) for the effective area of both obstacles. Pedestrians are initially distributed in a rectangular region with constant density ρ = 0.80, as shown in Fig. 10. The total number of of pedestrians is 44. The initial direction is θ 1 . Fig. 11 and 12 show the evacuation process for configuration 1 when α = 1 and α = 0 in the effective area, respectively. As explained in Sec. 2.1.2, when α = 1 (resp., α = 0) pedestrians avoid the front (resp., back) part of the effective area. Moreover, the shape of the real obstacle is different: it is square for α = 1, while for α = 0 it is slender. These findings are confirmed by Fig. 13 and 14, which show the evacuation process for configuration 2 when α = 1 and α = 0 in the effective area, respectively. Finally, Fig. 15 compares the evacuation times for the room with no obstacles (α = 1 everywhere in the domain), and for configurations 1 and 2 for α = 1 and α = 0 in the effective area. Obviously, the shortest evacuation time is for the room with no obstacles and overall good quality of the  Figure 11: Configuration 1 with α = 1 in the effective area: computed density for t = 0, 3, 6, 7.5, 9, 13.5 s. The small square within the effective area represents the real obstacle. environment. The evacuation time is slightly larger when there is one or two obstacles and α is equal to 1 in the effective area. When we compare Fig. 11 and Fig. 12 at t = 13.5 s, it seems that roughly the same number of people is left in the room. However, the small group of people in the effective area in Fig. 12 has velocity modulus close 0 since α is equal to 0 there and thus the overall evacuation process takes longer. The same happens in Fig. 14.

Lane formation
In this section, we assess the ability of our model to reproduce an empirically observed phenomenon: formation of lanes in a corridor when two or more groups of pedestrians have opposite walking directions [23].  Figure 12: Configuration 1 with α = 0 in the effective area: computed density for t = 0, 3, 6, 7.5, 9, 13.5 s. The small rectangle within the effective area represents the real obstacle.
directions θ 1 and θ 5 . See the top left panel in Fig. 16 and 17. The rest of Fig. 16 and 17 shows the evolution of the pedestrian dynamics. From both pictures, we see that pedestrians try to avoid contact by changing the direction, which leads to sorting, separation, and lane formation.

Conclusion
We considered a kinetic theory approach to model pedestrian dynamics in bounded domain and adapted it to handle obstacles. For the numerical approximation of the solution to our model, we applied the Lie splitting scheme which breaks the problem into two pure advection problems and a problem involving the interaction with the environment and other pedestrians.
Several test cases have been considered in order to show the ability of the model to reproduce qualitatively: -evacuation from a room with one exit, without and with obstacles; -evacuation from a room with two exits and no obstacles; -lane formation. In the case of the room with two exits and no obstacles, we also presented a quantitative comparison with experimental data. Numerical results and experimental data are in very good agreement for medium and medium-to-large groups of people. With the confidence in the model given by the experimental validation, we performed numerical tests to study evacuation for different scenarios