Parameter estimation for macroscopic pedestrian dynamics

models from microscopic data

∗ † ‡
Susana N. Gomes Andrew M. Stuart Marie-Therese Wolfram

March 26, 2019

arXiv:1809.08046v2 [math.AP] 22 Mar 2019

Abstract
In this paper we develop a framework for parameter estimation in macroscopic pedestrian
models using individual trajectories – microscopic data. We consider a unidirectional flow
of pedestrians in a corridor and assume that the velocity decreases with the average density
according to the fundamental diagram. Our model is formed from a coupling between a
density dependent stochastic differential equation and a nonlinear partial differential equa-
tion for the density, and is hence of McKean–Vlasov type. We discuss identifiability of the
parameters appearing in the fundamental diagram from trajectories of individuals, and we
introduce optimization and Bayesian methods to perform the identification. We analyze
the performance of the developed methodologies in various situations, such as for different
in- and outflow conditions, for varying numbers of individual trajectories and for differing
channel geometries.
Keywords: Macroscopic pedestrian models, generalized McKean–Vlasov equations, pa-
rameter estimation, optimization-based and Bayesian inversion.
AMS subject classifications:

1 Introduction
The complex dynamics of large pedestrian crowds as well as constant urbanization has initiated
many developments in the field of transportation research, physics, urban planning and more
recently applied mathematics. Originally, model development was mostly based on empirical
observations. More recently a plethora of data, from video surveillance cameras as well as
experiments, is starting to become available. This data, usually individual trajectories, is used to
calibrate individual based models or to identify characteristic relations, such as the fundamental
diagram, which have the potential to be observed in large pedestrian crowds. In this paper we
study such a calibration approach, setting up parameter estimation methodologies which allow
us to estimate parameters in macroscopic pedestrian models using individual trajectories.
There is a rich literature on mathematical models for pedestrian dynamics, see [12]. Indi-
vidual based, also known as microscopic models, are among the most popular in the engineering
and transportation literature. In these models individuals are characterized by their position
and sometimes velocity, see for example [23]. On the macroscopic level the crowd is described
by a density. These models are often derived from first principles, since the rigorous transition
from the micro- to the macroscopic model is still an open problem in certain scaling regimes.
∗
University of Warwick, Coventry CV4 7AL, UK (susana.gomes@warwick.ac.uk)
†
California Institute of Technology, 1200 E. California Blvd, Pasadena, CA 9112.
‡
University of Warwick, Coventry CV4 7AL, UK and RICAM, Austrian Academy of Sciences, Altenbergerstr.
66, 4040 Linz, AT

1
 (0, `) ΓN (L, `)

Γin Γout
x2

x1
(0, −`) ΓN (L, −`)

(a) Sketch of the domain (b) Sample of 5 generated trajectories.

Figure 1: Modeling setup considered

Then more general interactions, as proposed in [1, 8, 21, 25], can be considered. In this pa-
per we focus on mean-field interactions only. Here individual dynamics are influenced by the
averaged behavior of the crowd which can be characterized by the average velocity, flow or
density. Often, for example in high density regimes, individual dynamics have little influence
on the overall flow. For this reason averaged quantities are of specific interest to characterize
the crowd behavior. While individual based models depend on a large number of parameters,
which allow a detailed description of individual trajectories, their identification and estimation
is extremely challenging and computationally costly. Hence we focus on the identification of
macroscopic relations from individual trajectories in the following.
An important and commonly used macroscopic relation is the so-called fundamental dia-
gram. The fundamental diagram relates an averaged observed pedestrian density to either the
measured velocity or outflow. It is a well-established characteristic quantity in vehicle traffic
flow theory, see [30], but also plays an important role in pedestrian dynamics. Examples of its
use include the quantification of the capacity of pedestrian facilities or the evaluation of pedes-
trian models, see [39]. To understand the latter we may use trajectories which are collected
in controlled experiments. These experiments are conducted for various conditions – different
domains (corridors, junctions, ...), uni- and bidirectional flows, varying inflow and outflow con-
ditions, in addition to considering the effects of diverse social and cultural backgrounds; see
for example [5]. In experiments trajectory recordings usually start once the experiment has
equilibrated. The collected video or sensor data is used to extract individual trajectory data.
More recently also images from motion sensing cameras, such as Kinect have been used, see
[10], to collect individual trajectories in public spaces. The fundamental diagram is then cal-
culated by evaluating the individual trajectory data using suitable averaging techniques. In
unidirectional flows different regimes can be observed – at very low densities pedestrians walk
with a (maximum) velocity (denoted by vmax later on). Then the average speed decreases as
the density increases. At a certain point the velocity approaches zero due to overcrowding. We
refer to the density at which overcrowding occurs as ρmax later on. Different values for ρmax
and vmax can be found in the literature, ranging from 3.8 to 10 pedestrians per m2 for ρmax and
0.98 m/s to 1.5 m/s for vmax . These deviations can be explained by the experimental setup and
measurement techniques used, as well as psychological and cultural factors, cf. [7, 39, 43].
In this paper we consider a unidirectional flow in a corridor as illustrated in Figure 1a. We
assume that pedestrians enter the corridor on the left through Γin and exit on the right through
Γout , with no possible entrance or exit at the walls ΓN . The evolution of the overall density
ρ = ρ(x, t) can be modeled by a nonlinear Fokker–Planck equation, with (nonlinear) mixed

2
boundary conditions, describing the inflow and outflow at the entrance and exit. Here the drift
term accounts for the nonlinear velocity dependence on the density, as described by the funda-
mental diagram. The scaled density ρ allows us to describe individual trajectories X = X(t) as
realizations of a generalized McKean–Vlasov process. This introduces a coupling between the
(microscopic) trajectories X and the (macroscopic) density ρ of the system, with X governed
by a density dependent stochastic differential equation (SDE) (6) and ρ governed by a Fokker-
Planck type partial differential equation (PDE) (1). It is important to note that the density ρ
is not a probability density function, since the total mass changes in time due to the inflow and
outflow boundary conditions. Rather than directly employ the empirical density of the process
to perform parameter estimation we instead try to learn parameters from individual trajectories
whose guiding velocity depends on the density via the fundamental diagram; this avoids the
difficult problem of creating smooth densities from collections of individual trajectories. We
study the coupling between the PDE and the trajectories of each pedestrian with the goal of
learning two parameters from individual trajectories of the resulting nonlinear Markov process:
the maximum speed vmax and the maximum density ρmax . Whilst there is some work concern-
ing parameter estimation of a nonlinear Markov process this is primarily focused on settings in
which the data concerns macroscopic or mean properties [20]; see also the recent analyses of
nonlinear Markov processes in [17, 18]. In contrast we focus on learning from microscopic data;
closely related work concerns the learning of (macroscopic) ocean properties from Lagrangian
(microscopic) trajectories [2].
We will see that the maximum density ρmax is not identifiable in our setup. The parameter
vmax is learnable and will be determined using a maximum a posteriori (MAP) estimator as
well as a fully Bayesian methodology. Estimates such as MAP (as well as the related maximum
likelihood estimator, MLE) are asymptotically unbiased and asymptotically normal as long as
the path observations are sufficiently long (i.e., observe the evolution of the generalized McKean–
Vlasov process for T → ∞) [28, 36]. However, since individuals leave the domain after a finite
time, we use multiple trajectories and show that this leads to similar properties to those found
when observing one very long trajectory of an SDE.
The considered mathematical model is rather simplistic from a modeling perspective. At
the moment we do not account for many known features of individual behavior, such as collision
avoidance or small group dynamics. However, the proposed setup is an appropriate model for
the density of the process and still mathematically tractable. The developed framework will
allow us to extend and generalize our parameter learning approach to more complicated models,
geometries, and situations in the future. Furthermore we observe that the SDE model considered
describes individual motion only. The (Brownian) noise models the tendency of individuals to
not move along completely straight lines. However, it does not correspond to measurement
errors in experiments. Including such errors would involve an extra equation for the trajectory
observations and requires the use of filtering methods. This problem could be considered in the
future; however it is likely that the model error itself is larger than the observation error in such
situations, and investment in refining the model is arguably more important than accounting
for noise in the trajectory data.
This paper builds a systematic methodology to identify parameters of interest in systems
which are modeled by SDEs which depend on their own density. The methodology is robust to
improvements in the model, and can be extended to more complex sensing of individual trajec-
tories, and to the incorporation of noise into the trajectory data. This approach can also be used
to identify functions instead of parameters and quantify uncertainty in our estimates. There has
been an increasing interest in developing connections between trajectory data and mathematical
models; see for example [4]. However, in our approach the likelihood is constructed explicitly,
using properties of stochastic differential equations, allowing for fully Bayesian inference, and

3
avoiding the need for approximate Bayesian inference, as employed in the recent work by Bode
et al. [4].
Our main contributions are the following:
• set up of parameter estimation methodologies for SDEs which depend on their own density,
with trajectories confined to bounded domains via boundary conditions;

• discussion of identifiability in the considered model of pedestrian dynamics;

• development of numerical methods which allow for parameter estimation via optimization
and for uncertainty quantification via a Bayesian approach;

• analyzing the impact of flow conditions (inflow vs. outflow rate, transient vs. steady state
profiles) on the parameter estimation;

• comparison between transient and steady state density regimes, providing novel insights
about the applicability of the fundamental diagram in transient regimes (we recall that
the fundamental diagram is obtained from equilibrated experimental data).
The paper is organized as follows. We start by introducing the underlying mathematical
models in Section 2. Next we present analytic results for the PDE and SDE models in Section 3
and introduce the parameter identification problem in Section 4. The computational complexity
of the considered identification problem requires efficient and robust numerical solvers for the
PDE and SDE models. We discuss these solvers as well as the optimization methodologies used
in Section 5. We conclude by discussing the identifiability of vmax in various density regimes
and for different experimental setups with numerous computational experiments in Section 6.
We present our conclusions in Section 7. Further details on the analysis of the PDE will be
given in the Appendix.

2 The microscopic and the macroscopic models

We consider the uni-directional flow of a large number of pedestrians moving in a corridor
Ω = {(x1 , x2 ) ∈ [0, L] × [−`, `]} as illustrated in Figure 1a. On the macroscopic level the
evolution of the density ρ(x, t) can be described by a nonlinear Fokker–Planck (FP) equation

∂t ρ(x, t) = div(Σ∇ρ(x, t) − ρ(x, t)F (ρ)), (1)

where the diffusion matrix is of the form Σ = diag(σ12 , σ22 ). The diffusion accounts for the
tendency of pedestrians to not move along perfectly straight lines, and its structure allows for
different diffusivities in the vertical and horizontal direction. Since individuals move from the
left to the right, we choose a convective field F : R 7→ R2 of the form
 
+ + 1
F (ρ) = f (ρ)e1 , where f : R 7→ R and e1 = . (2)
0

This convective field depends on density only, and will be evaluated pointwise to drive individual
trajectories.1 We choose the velocity as
 
ρ
f (ρ) = vmax 1 − , (3)
ρmax
1
Later we will replace e1 , the direction of motion, by a gradient vector field to allow for more complex domains
such as bottlenecks; this formulation will be used in the existence and uniqueness proof of the PDE as well as
the bottleneck example presented in Section 6.

4
which corresponds to the density-flow relation observed in the fundamental diagram in traffic
flow and pedestrian dynamics. Thus the model states that all individuals want to move at a
maximum speed vmax in the absence of other pedestrians, and the actual velocity decreases
linearly from this value with the density ρ evaluated locally. We assume that the corridor is
initially empty, that is ρ(x, 0) = 0 for all x ∈ Ω. Individuals enter at a certain rate a on the
left through the entrance Γin and leave on the right through the exit Γout with a rate b. On the
rest of the boundary we impose no flux boundary conditions. Let j = −Σ∇ρ + F (ρ)ρ denote
the flux. Then the corresponding boundary conditions are given by


j · n = −a ρmax − ρ , for all (x1 , x2 ) ∈ Γin , (4a)
j · n = bρ, for all (x1 , x2 ) ∈ Γout , (4b)
j · n = 0, for all (x1 , x2 ) ∈ ΓN , (4c)

where

Γin = {(x1 , x2 ) ∈ Ω : x1 = 0} , Γout = {(x1 , x2 ) ∈ Ω : x1 = L}

ΓN = {(x1 , x2 ) ∈ Ω : x2 = ±`}

and n denotes the unit outer normal vector. Note that the inflow condition (4a) includes the
additional factor (ρmax − ρ) due to volume exclusion. This prefactor arises in the formal limit,
when particles are only allowed to enter or move to a certain position if enough physical space is
available, see [42]. We note that a and b are rates of entrance and exit in the domain. Balancing
the left- and righthand sides of equations (4a) and (4b) gives us their units – both a and b are
in m/s. A small modification of the maximum principle calculations in [6] shows that it is
sufficient for the problem to be well posed that 0 ≤ a, b ≤ vmax .
In the following we will see that trajectory data generated by the model (1) is independent
ρ
of the parameter ρmax . Let ρ̃ = ρmax , then equation (1) can be rescaled as
 
∂t ρ̃ = div Σ∇ρ̃ − ρ̃F̃ (ρ̃) (5a)

where

F̃ (ρ̃) = f˜(ρ̃)e1 , f˜(ρ̃) = vmax (1 − ρ̃) . (5b)

The system is supplemented with the boundary conditions

j̃ · n = −a 1 − ρ̃ , for all (x1 , x2 ) ∈ Γin , (5c)
j̃ · n = bρ̃, for all (x1 , x2 ) ∈ Γout , (5d)
j̃ · n = 0, for all (x1 , x2 ) ∈ ΓN . (5e)

Here the scaled flux is given by j̃ = −Σ∇ρ̃ + vmax (1 − ρ̃)ρ̃e1 . We see that the maximum density
ρmax is not present in the scaled formulation. Furthermore the convective field F , which governs
individual trajectories, is given by (2), (3) and can be expressed entirely in terms of ρ̃, with no
reference to ρmax .
Our stated aim is to identify vmax and ρmax in f (ρ) (3) using individual trajectories. How-
ever the preceding argument shows that the parameter ρmax does not influence the SDE for
trajectories. Since our parameter inference is based only on these trajectories this means that
ρmax cannot be learned from the data available to us. Hence we make the following remark:

5
Remark 1. The parameter ρmax can not be identified within our adopted microscopic macro-
scopic data-model framework. This limitation is not caused by the identification methodologies
proposed, but rather by the invariance of the PDE-SDE model to scaling in ρ and the fact that
the measured trajectory data is independent of the value of ρmax . Indeed we initially conducted
numerical experiments using the unscaled model, which led to the understanding that ρmax is not
identifiabile. Similar identifiability issues are well-known in the literature relating to inference
for diffusion processes; see for example the paper [37] in which it is shown that the quadratic
variation of sample paths cannot contain information about time-rescaling.
For this reason the rest of this paper is based only on the scaled Fokker–Planck equation
for ρ̃. However we wish to drop the ∼ notation for ease of presentation. This corresponds to
using equations (1)–(4) with ρmax = 1. Since the preceding arguments show that ρmax is not
identifiable from trajectory data, making any specific choice of ρmax , including the value 1, will
not affect the inference for vmax . The parameter ρmax thus plays no further part in the paper.
Our focus, then, is on identification of vmax . We will use the scaled Fokker–Planck equation
(5) in the following, which coincide with the system (1)–(4) when setting ρmax = 1.
The existence of steady states as well as the different stationary regimes – so-called influx lim-
ited, outflux limited and maximum current regime – are discussed by Burger and Pietschmann
in [6]. We present corresponding existence results for the time dependent problem in Section
3. Note that the steady state as well as the time dependent solutions of (5) satisfy 0 ≤ ρ ≤ 1,
which ensures that f stays non-negative for all x ∈ Ω and times t > 0. This also allows
us to define individual trajectories as realizations of the following generalized McKean–Vlasov
equation
√
dX(t) = F ρ(X(t), t) dt + 2ΣdW (t), (6)
where W (·) is a unit Brownian motion and ρ solves equation (5). These realizations correspond
to individual trajectories as we see in Figure 1b.
In the presented PDE model the boundary conditions drive the dynamics of the process.
Hence they need to be included at the SDE level as well if we wish for self-consistency. Boundary
conditions for SDEs are a delicate issue and a general existence theory is not available in closed
domains. It is out of the scope of this paper to state the SDE problem in a rigorous manner
(e.g., using local time [41]). However, it is important to implement the boundary conditions
in a consistent way with the PDE, and we discuss their implementation in more detail in
Section 5.2. Once this is done we note that the scaled FP equation is self-consistent with the
forward Kolmogorov equation for this SDE. Indeed any sufficiently regular solution ρ of the
Fokker–Planck equation (5) can be used as a drift in (6). The corresponding Fokker–Planck
equation is then a linear PDE with the same solution as the nonlinear PDE (5). See [32, Thm.
1],[16, Thm. 7.3.1],[13, 19] for theory relating to such PDEs in related settings.

3 Analysis of the model

In this section we discuss the analysis of the generalized McKean–Vlasov SDE (6) as well as
the rescaled Fokker–Planck equation (5). We start with the existence and regularity results for
the Fokker–Planck equation, as they ensure the well-posedness of the corresponding generalized
McKean–Vlasov process. We re-emphasize the important fact, discussed prior to Remark 1,
that (5) is identical to equation (1)–(4), with ρmax = 1.

3.1 Analysis of the Fokker–Planck equation

The nonlinear boundary conditions (5c) and (5d) strongly influence the time dependent and
steady state solutions. The relation of the maximum velocity vmax to the inflow and outflow

6
Figure 2: Examples of the steady states of the Fokker–Planck equation in one dimension for
σ = 0.5, 0.1, 0.05 and (left) the influx limited phase – a = 0.2, b = 0.4, (center) the outflux
limited phase – a = 0.4, b = 0.2, and (right) the maximal current phase – a = 0.9, b = 0.975.

parameters a and b defines different stationary regimes, in which boundary layers arise at the
entrance and/or the exit.
The following regimes were introduced in [42] and were characterized and analyzed by Burger
and Pietschmann in [6]:

(1) Influx limited phase in the case a < b and min (a, b) < vmax
2 . We observe an asymptotically
1
low density (ρ < 2 ) and a boundary layer at the exit.
(2) Outflux limited phase for a > b and min (a, b) < vmax
2 . This creates an asymptotically high
density (ρ > 12 ) and a boundary layer at the entrance.
(3) Maximal current phase if a, b ≥ vmax
2 . Here we observe an asymptotic density with ρ ≈ 2
1

in most of the domain when σ is small and boundary layers at the entrance and exit. In
the case a = b = vmax 1
2 , the unique steady state is given by ρ(x) = 2 .

Figure 2 illustrates the different regimes in 1D in the case of a corridor of length L = 3,

vmax = 1.5 and different values of σ = σ1 . We observe the predicted boundary layers at the
entrance and/or exit, whose width is determined by the diffusivity σ. The existence of at
least one steady state solution under suitable assumptions on the in- and outflow boundary
conditions as well as a more general regular potential was proven in [6]. Note that these steady
state profiles are unique in 1D.
In the following we state and discuss the respective existence result for the time dependent
problem. The proof uses similar techniques as in [6] and can be found in the Appendix. Without
loss of generality we assume σ1 = σ2 = 1 and vmax = 1. Furthermore, we consider the more
general case in which e1 is replaced by ∇V in (2). We make the following assumptions:

(A1) Ω ⊂ Rn bounded, n = 1, 2, 3, with boundary ∂Ω ∈ C 2 .

(A2) 0 ≤ a, b ≤ 1, and min{a, b} < 1.

(A3) V ∈ W 1,∞ (Ω) and ∆V = 0.

Assumption (A1) actually excludes domains like the corridor shown in Figure 1a. We still
employ the boundary conditions (5c)–(5e), simply assuming that the union of (the disjoint)
Γin , Γout and ΓN comprise the entire boundary ∂Ω. One could consider a rectangular domain
with rounded off corners instead, but since we did not observe any issues in the computational
experiments with corners, we retained them for a more realistic simulation setting.

7
 We seek weak solutions in the set

S = {ρ ∈ R : 0 ≤ ρ ≤ 1}, (7)

and refer to its interior as S 0 . Then the weak solution ρ : Ω × [0, T ] → S satisfies
Z T  Z Z Z 
h∂t ρ, ϕiH −1 ,H 1 − j · ∇ϕ dx − a (1 − ρ)ϕ ds + b ρϕ ds dt = 0, (8)
0 Ω Γin Γout

for all test functions ϕ ∈ L2 (0, T ; H 1 (Ω)).

Theorem 1. Let assumptions (A1)-(A3) be satisfied. Let ρ0 : Ω → S 0 be a measurable function

such that the initial entropy E defined by (20) satisfies E(ρ0 ) < ∞. Then there exists a weak
solution to system (5) in the sense of (8) satisfying

∂t ρ ∈ L2 (0, T ; H 1 (Ω)∗ )
ρ ∈ L2 (0, T ; H 1 (Ω)).

In the preceding theorem, the star denotes the dual operation. The existence proof as well
as improved regularity results can be found in the Appendix.

3.2 Analytic considerations on the generalized McKean–Vlasov equation

The analytic results for the Fokker–Planck equation (5) would allow us to prove existence
and uniqueness of strong solutions for (6) in R2 . However, since we consider the setup of a
corridor with a mix of boundary conditions (reflected in the walls, and partially reflected in
the entrance and exit), standard results are not applicable. The only available results for SDEs
with boundaries are slight generalizations of Skorokhod’s problem [41] and can only be used
for SDEs in the half plane. We are not aware of any results for more complicated geometries.
However, our numerical experiments show that a careful discretization of the process, and
in particular treatment of trajectories impinging on the boundary of Ω, yields bounded and
continuous solutions of equation (6).

4 Parameter estimation
In this section we introduce the parameter estimation framework for SDEs which depend on
their density. We will use a Bayesian approach based on sampling from a posterior distribution,
and confirm our results with the computation of a maximum a posteriori (MAP) estimator.
We recall that our initial goal was to estimate both ρmax and vmax . However, as we have seen,
determination of ρmax is not possible using the current setting and only determination of vmax
is feasible.
We estimate v := vmax from a collection of sample paths {Xi (t)}i=1,...,J
t∈[0,T ] which are realizations
of the McKean–Vlasov equation (6). In what follows we parameterize f (·) by v and write f (·; v).
For ease of presentation, we discuss the estimation from a single trajectory before generalizing
it for multiple trajectories.
Let X = X(t) be a realization of (6). Throughout the remainder of the paper we use the
notation
1
| · |A = |A− 2 · |
where | · | is the Euclidean norm and A any positive-definite symmetric matrix; a corresponding
inner-product may be defined by polarization. Since Ẇ is white, it is intuitive that finding the

8
best value of the parameter v, given an observation of a trajectory X, corresponds to minimizing
the function
1 T
Z
|Ẋ − F ρ(X(t), t); v |2Σ ,


Φ(v; X) = (9)
4 0
over all possible values of v. However, the function Φ(v; ·) is almost surely infinite. In order to
avoid this problem, we note that we can write
1 T
Z
Φ(v; X) = Ψ(v; X) + |Ẋ|2 dt, (10)
4 0
where the function Ψ is defined as
1 T
Z 
|F ρ(X(t), t); v |2Σ dt − 2hF ρ(X(t), t); v , dX(t)iΣ .



Ψ(v; X) = (11)
4 0
Note that the last term in Ψ is an Itô stochastic integral. It is preferable to perform the
parameter estimation using Ψ instead ofRΦ, since the latter is infinite almost surely, whilst the
T
former is finite almost surely, and since 0 |Ẋ|2 dt does not depend explicitly on v both cases
yield the same results when applied in discretized form.
We will also add prior information and perform the parametric estimation of vmax by mini-
mizing
1
J (v; X) := Ψ(v; X) + |v − m|2 , (12)
2c
subject to the constraint that v is positive. The parameters m and c correspond, respectively, to
a prior estimate of the velocity and the variance associated with this estimate. This optimization
problem will now be derived through a detailed description of the Bayesian formulation of
inversion. Specifically we define precisely the Bayesian problem for the posterior distribution
of v given a trajectory {X(t)}t∈[0,T ] ; this posterior distribution is maximized at the constrained
minimization problem just defined.
In this Bayesian context the functional J can be interpreted as follows. The first term
measures the misfit between the observed trajectory and the predicted density regime obtained
from the FP equation. The second term is a regularization term, which is weighted by prior
knowledge on v. In particular, we assume that v is normally distributed with mean m and
variance c, conditioned to be positive. The probabilistic interpretation of J is as follows. The
function exp(−J (v; X))1(v > 0), appropriately normalized, is a probability distribution and in
fact the posterior distribution, P(v|X), of v given a realization X. This conditional distribution
can be obtained, up to a normalization constant which is independent of v, by computing the
joint probability of the process X and the random variable v, which is given by
P(X, v) = P(v|X)P(X);
thus
P(v|X) ∝ P(X, v), (13)
where the constant of proportionality is independent of v. Consider the probability measure
λ = Q ⊗ Leb(R), where Q is the law of the Brownian motion driving the process and Leb(R) is
the Lebesgue measure in R. Lemma 5.3 in [22] states that
dP dP
(X, v) = (X|v)π0 (v), (14)
dλ dQ
where π0 is the prior Lebesgue density for v. Using Girsanov’s Theorem [14, 33], we see that
dP
(X|v) = exp(−Ψ(v; X)), (15)
dQ

9
where Ψ is defined in (11). If we choose π0 (dv) = 1(v > 0)N (m, c)(dv), which ensures that the
velocity is almost surely positive both a priori and hence a posteriori, equation (13) gives
 1 
P(v|X) ∝ exp −Ψ(v; X) − |v − m|2 1(v > 0), (16)
2c
which is exactly exp(−J (v; X))1(v > 0).
We can either sample from the probability distribution given by (16) (the Bayesian ap-
proach) or we can simply minimize the negative log likelihood of the data J (v; X), subject to v
being positive. While the fully Bayesian framework allows us to quantify the uncertainty on our
estimate, the minimization approach only gives a single point estimator. It is however, much
cheaper to simply minimize. In the following sections we will compare the results of the two ap-
proaches by sampling from the posterior distribution using the pCN algorithm, and minimizing
the function J using the Nelder–Mead algorithm. Both methods are derivative free, so they
do not involve any extra function evaluation due to computing derivatives. We note that other
derivative free samplers or optimizers could be used.
The accuracy of the aforementioned estimators can be expected to improve with the lifetime
of the observations [36]. Since in our setting trajectories terminate when they exit the domain,
we are only able to observe them for a finite time. We will see that using multiple trajectories
to estimate parameters has similar desirable properties to observing a single trajectory of a
standard (boundaryless) SDE over a long time horizon.
We denote a family of such trajectories by Xj (t), j = 1, . . . , J. We assume that each trajec-
tory Xj starts at tj0 with Xj (tj0 ) = Xj0 and exits at time t = tjf . We set F (ρ(Xj (t), t); vmax ) = 0
/ [tj0 , tjf ], which allows us to define
for t ∈
Z T
1 
|F ρ(Xj (t), t); v |2Σ dt − 2hF ρ(Xj (t), t); v , dXj (t)iΣ .



ψj (v; Xj (t)) =
4 0

Then we can sum the misfit of all trajectories and obtain

J
X
ΨJ (v; X1 , . . . , XJ ) = ψj (v; Xj ). (17)
j=1

Since all ψj are independent, we can apply the previous argument, replacing Ψ(v, X) by (17).
The averaging over all the trajectories has the same effect as considering large values of T for
a single trajectory. This can be shown formally, using techniques similar to those used in [28]
for a single trajectory.
It could be of interest to also estimate the diffusion coefficient matrix Σ. Estimating con-
stant diffusion coefficients for sufficiently frequent observations is, in theory, a well-understood
problem, see [26, 36]. However our model is likely to be inconsistent with the data at small
scales and so it is important to appreciate that estimation of Σ might well be a non-trivial
problem, requiring techniques such as those introduced in [34].

5 Computational methods
In this section we discuss the discretization of the nonlinear PDE (5) and SDE (6) before con-
tinuing with parameter estimation in the following section. Note that the parameter estimation
is computationally costly in the sense that it involves the solution of a PDE in every iteration of
the sampling and optimization algorithms. Motivated by relevant experimental settings (noting
that data collection starts after the experiment is equilibrated), we explore the effect of using

10
both stationary and transient density regimes on the quality of our estimates. Using steady
state density profiles also has the advantage of a lower computational cost. However we will
show that the steady state regimes do not produce consistent estimates across all parameter
regimes, an effect which is mitigated by using (the more expensive) time-dependent density
based estimation.

5.1 The time dependent and steady state Fokker–Planck solvers

The time dependent solver is based on the decomposition of (5a) into a diffusive and a convective
part. For the nonlinear convective operator an upwind discontinuous Galerkin (DG) method
and an explicit time integration scheme is used. For the discretization of the linear diffusive part
we use a hybrid discontinuous Galerkin (HDG) method and an implicit time discretization. This
additive time splitting, also known as Implicit-Explicit (IMEX) method, allows to treat the stiff
diffusive term implicitly while the nonlinear non-stiff hyperbolic problem is solved explicitly,
see [3]. The IMEX scheme allows for larger time steps, while the (H)DG discretization ensures
stability and parallelizablilty.
In some cases we employ 1D steady state profiles ρs = ρs (x) (extended as constants in the
orthogonal direction) and these are calculated using an H 1 conforming finite element discretiza-
tion. We use a damped Newton method to solve the resulting nonlinear system, with an initial
guess depending on the relation of a and b. Both solvers use the finite element software package
Netgen/Ngsolve, see [38].

5.2 SDE solver and trajectory generation

We approximate all trajectories of (6) by using the explicit Euler-Maruyama scheme [24], which
is a forward Euler time-stepping method. It consists of defining the approximation Xk ≈
X(tk ) := X(k∆t) and updating the individual position as
√
Xk+1 = Xk + F (ρ(Xk , tk ))∆t + 2Σ∆tξk , (18)

where ξk ∼ N (0, 1) and ρ is the solution to the Fokker–Planck equation (which we advance
in time simultaneously). To generate trajectories using the steady state profiles, we replace
ρ(Xk , tk ) by the steady state solution ρs (Xk ).
Since the diffusion coefficient is constant, the Euler-Maruyama scheme has strong order of
convergence one [24], and we do not need to use more complicated time-stepping methods. At
the boundaries we implement individual rules consistent with (5). Singer et al [40] and Erban
and Chapman [15] derived appropriate scaled rates of entrance and exit for partially reflected
boundaries in the half plane. These ensure that the Robin boundary conditions are satisfied
by the corresponding probability densities. Even though it is not established that their results
are directly applicable in the case of the more complex geometries we use here, we nonetheless
follow their approach, and implement the entrance and exit behavior as described in Table 1.

5.3 Derivative-free optimization – the Nelder–Mead algorithm

As detailed in Section 4, we will use derivative-free methods for the minimization of (12), corre-
sponding to finding the MAP estimator of the posterior probability distribution: the maximizer
of (16). We choose the Nelder–Mead algorithm [31], which is a widely used nonlinear uncon-
strained optimization algorithm. It is a direct search method which is based on the function
evaluation of n + 1 points in an n−dimensional space. A rough outline of the method (see [9] for
more details) is given in Algorithm 1. The number in boldface is the total number of function
evaluations necessary if that particular step is performed. If the objective function is strictly

11
 Boundary PDE SDE Implementation

ΓN Neumann Reflecting
p=1

Γin Robin Wait to get in p = Pin

(from outside) p = 1 − Pin

p = 1 − Pin

Γin Partially p = Pin

Robin
(from inside) reflecting back to Γin

p = 1 − Pout
Partially p = Pout
Γout Robin
reflecting leave domain

Table 1: Implementation of the boundary conditions at the SDE level. The filled dot represents
the current state, while the empty circles represent the possible position in p
the following time
step, which is realized
p with probability p. Here we use the rescalings P in = π∆t/(2σ12 )a(1 −
2
ρ(0, x2 )) and Pout = π∆t/σ1 bρ(L, x2 ) according to [15].

Algorithm 1 Nelder–Mead algorithm

1: Choose an initial simplex {v 0 , . . . , v n } and compute J i = J (v i ; X).
2: while Diameter of simplex > tolerance, do
3: Order the n + 1 vertices so that J 0 ≤ J 1 ≤ · · · ≤ J n .
4: Reflect: reflect the worst vertex v n over the centroid of the simplex, obtaining v r . If
J 0 ≤ J r < J n−1 , replace v n by v r . [1]
5: If J r < J 0 , Expand the reflected point, obtaining v e .
If J e ≤ J r ,replace v n by v e . Otherwise replace v n by v r . [2]
6: If J r ≥ J n−1 , Contract: the reflected point, obtaining v c .
If v c improves on J n then replace v n by v c .[2]
7: Otherwise, perform a Shrink. [n+1]
8: end while

12
convex then the algorithm only requires 1 or 2 function evaluations per iteration. Furthermore,
it converges in 1D for strictly convex functions [29]. In 2D only limited convergence results
are known. We use Matlab’s implementation of the Nelder–Mead algorithm (via the function
fminsearch) with a modification of J that penalizes negative values of v. We observe a large
decrease of the objective function in the first few iterations, as reported in the literature for
many other applications.

5.4 MCMC method – the pCN algorithm

The MAP estimator obtained with the Nelder–Mead algorithm does not provide us with any
quantification of uncertainty on our estimates. In order to obtain such uncertainty estimates
we sample from the posterior distribution given by (16). For this purpose we use the pCN
(pre-conditioned Crank–Nicholson) method described in [11]; it is a general purpose method for
sampling from any distribution which is formed as the product of a (not necessarily smooth)
function f and a Gaussian density and has the advantage that derivatives of f are not required. In
our problem this means that we avoid computing derivatives of Ψ. Looking forward, the method
also has the potential to scaling up to cases where the objective is to estimate the shape of f (ρ)
in (3) nonparametrically. Furthermore, the algorithm has a single tuneable parameter β which
can be explored to maximize efficiency. Algorithm 2 gives the details.

Algorithm 2 The pCN algorithm

1: Set k = 0 and pick v (0) .
2: for k = 1, ..., N , wherepN is the number of iterations, do
3: (k)
Propose y = m + (1 − β 2 )(v (k) − m) + βξ (k) , ξ (k) ∼ N (0, c).
4: Set v (k+1) = y (k) with probability αk := α(v (k) , y (k) )
5: Set v (k+1) = v (k) otherwise.
6: k → k + 1.
7: end for

The acceptance probability is given by

α(v, y) = min{1, exp Ψ(v; X) − Ψ(y; X) }1(y > 0). (19)

Notice, in particular, that if Ψ(y (k) ; X) ≤ Ψ(v (k) ; X) and y (k) > 0 then the proposed move is
accepted with probability one.

6 Numerical results
In this section we present and discuss numerical results based on the solvers and methodologies
presented in Section 5. In particular we want to estimate vmax using multiple trajectories of
the coupled SDE-FP system. These trajectories are generated using both the time-dependent
and the steady state solution of the FPE, and for different inflow and outflow conditions. The
corresponding parameter estimation problem is then solved using the Nelder–Mead optimizer
or the pCN Bayesian sampler. Our numerical simulations lead to the following conclusions:

1. the maximum speed vmax is learnable; see Figure 3

2. results are more accurate and less uncertain if the influx and outflux rate a and b differ
considerably; see Figure 11

13
 3. we do not obtain reliable estimates for all regimes using stationary density profiles; see
Figure 7

4. estimates obtained from trajectories experiencing steady state densities are less uncertain
than from time dependent ones (in the regimes where steady state estimation does give
consistent results); see Figures 11 and 12.

In order to test the parameter estimation methodology, we generate a collection of trajecto-

ries for a range of parameters which represent the different density regimes. To create this data
we choose the following parameters for all trajectories: corridor length L = 3, corridor width
2` = 0.5, noise strength σ1 = σ2 = σ0 = 0.05 and a final time T = 2. The time steps in the
SDE solver are set to ∆tSDE = 10−3 . We assume that the true value for the maximum velocity
is vmax = 1.5 m/s, which is chosen close to values obtained from experiments. The final time
is chosen sufficiently large to ensure that some individuals have the time to leave the domain.
The size of the time steps corresponds to the setting in a high frequency regime. We generate
trajectories using a time dependent density ρ = ρ(x, t) satisfying (5a) as well as the steady state
ρs = ρs (x). Note that the latter choice is computationally less expensive and corresponds more
closely to widely used experimental conditions, since data is typically collected once the pedes-
trian flow has equilibrated. We consider five sets of parameters which span the three steady
state density regimes:

- Outflux limited: a = 0.4, b = 0.2 and a = 0.45, b = 0.4,

- Influx limited: a = 0.2, b = 0.4 and a = 0.1, b = 0.15,
- Maximal current: a = 0.9, b = 0.975.

Since the variance coefficients might not be known in practice, we will estimate the pa-
rameters using values for the diffusion coefficient which differ significantly from those present
in the data; we take σ1 = σ2 = σ = 1 in our algorithm. This corresponds to a form of
model-mispecification and avoids committing an inverse crime [27].

6.1 Benchmarking
In the following, we present a first set of numerical results which confirm that the optimization
and MCMC methodologies perform as expected. In particular, we will see that the value of
the estimator for vmax is consistent across both methodologies, and is also independent of the
various parameters used, such as time steps or spatial mesh, and the parameter β of the pCN
algorithm. Throughout this subsection we use J = 20 trajectories.

Influence of the discretization This test case is of particular interest with respect to the
computational efficiency, since each iteration of either minimization algorithm requires one (or
more) PDE solves. Therefore we wish to use as large time steps as possible as well as a mesh
which is as coarse as possible. We first present numerical results related to the influence of the
PDE time step on the value of the estimator for vmax . We run the pCN (and Nelder–Mead)
algorithm using three different time steps for the PDE solve: ∆t = 0.001, ∆t = 0.005 and
∆t = 0.01 and interpolate the obtained values to fit the time step of the SDE trajectories. We
tested the influence of the time step in the maximal current regime with a = 0.9, b = 0.975.
The posterior distributions for vmax as well as the corresponding prior distribution and the
true value are depicted in Figure 3. These results are consistent with those obtained using the
Nelder–Mead algorithm to compute the MAP estimator (successive iterations of the Nelder–
Mead algorithm are depicted in the right panel with orange circles, and we point out that the

14
Figure 3: Influence of the PDE time step on the estimator for a maximal current regime with
a = 0.9, b = 0.975. Left: true value (gray dotted line), prior distribution (full gray line), and
posterior distributions for different PDE time steps. Right: prior mean (full gray line), true
value (gray dotted line), and posterior mean (orange dotted line) and MAP estimator iterations
(orange circles) for ∆t = 0.005.

number of iterations for the MAP estimator was scaled by a factor of approximately 10−3 for
comparison). We observe that all the test cases produce similar results, with a slightly better
agreement in the case of the finer time steps, where the results are also less uncertain, since the
posterior distribution has a smaller variance. For this reason, all the pCN and Nelder–Mead
simulations presented below will be performed with a PDE time step of ∆t = 0.005.
We observe a similar agreement between the mean of the posterior distribution and the MAP
estimator for all the results presented below, henceforth we present the posterior distribution
only. We also point out that the MAP estimator (mode of the posterior distribution) is approx-
imately the same as the running average (mean of the posterior distribution). This suggests
that the posterior distribution is approximately Gaussian as is to be expected when the data
is highly informative. Finally, we also tested the influence of the SDE time step (frequency of
observations) and spatial mesh on the estimates, obtaining similarly good results.

Influence of the parameter β The value of the estimator obtained using the pCN algorithm
should be independent of the parameter β since β is a tuneable parameter in the algorithm
and not a parameter of the posterior distribution. In Figures 4 and 5 we plot the prior and
posterior distributions for vmax in two influx limited regimes, a = 0.2, b = 0.4 (left panels) and
a = 0.1, b = 0.15 (right panels). Figure 4 was obtained in a time dependent regime, while
Figure 5 uses trajectories experiencing a steady state density. We observe in both figures that
the posterior distribution is independent of the parameter β, as expected.

6.2 Influence of problem-specific parameters

In the following, we present our numerical results which confirm the conclusions presented at the
beginning of this section. In particular, we will show that both approaches give consistent results
with respect to the initial guess, prior mean and variance, and number of used trajectories. All
the examples presented use β = 0.1 and J = 20 trajectories, except the case in which we vary
the number of trajectories.

15
Figure 4: Influence of the parameter β on the estimator for two influx limited regimes and
using the solution of the time dependent Fokker–Planck equation. Left: a = 0.2, b = 0.4. Right:
a = 0.1, b = 0.15.

Figure 5: Influence of the parameter β in the posterior distribution of vmax guesses for influx
limited regimes and using the solution of the steady state Fokker–Planck equation. Left: a =
0.2, b = 0.4. Right: a = 0.1, b = 0.15.

Influence of the initial guess and prior mean We assume that vmax has prior distributions
N (m, c) with m ∈ {1, 2} and c = 0.25 and test using two initial guesses v = 1 and v = 2. Figure 6
depicts the prior distribution (gray full line), true value (gray dotted line) and the posterior
distributions for each case tested. The MAP estimator agrees with the means of the posterior
distributions for all cases. The left panel depicts the results for m = 1 while the right panel
has m = 2. In each case we observe that the estimates are independent of the initial guess, and
both prior distributions produce similar results. Next we use trajectories generated with steady
state density profiles. Figure 7 shows the influence of initial guess and prior mean for three
different regimes: an outflux regime with a = 0.4, b = 0.2 (left panel), an influx limited regime
with a = 0.2, b = 0.4 (middle panel) and a maximal current regime with a = 0.9, b = 0.975
(right panel). We use the same prior distributions and initial guesses as in Figure 6. We find
perfect agreement between the Nelder–Mead and pCN estimators, and the Bayesian estimator
has converged to the posterior distribution. The maximal current (right panel) and influx
limited (left panel) regimes produce similar results as the corresponding time dependent case;
however, in the outflux limited regime we can observe that the posterior distribution stays

16
Figure 6: Influence of prior distribution and initial guess for an outflux limited regime with
a = 0.45, b = 0.4 using solutions of the time dependent Fokker–Planck equation. Left: Prior
mean m = 1. Right: Prior mean m = 2. Both figures have c = 0.25 and depict the prior
distribution (full gray line), true value (dashed gray line) and posterior distributions for initial
guesses vmax = 1, 2.

Figure 7: Influence of prior distribution and initial guess for three different regimes using
solutions of the steady state equation. Left: outflux limited regime with a = 0.4, b = 0.2.
Middle: influx limited regime with a = 0.2, b = 0.4. Right: maximal current regime with
a = 0.9, b = 0.975.

close to its prior distribution. The reason for this may be found in the trajectories used – in
the outflux limited regimes, individuals experience very high densities and get ‘stuck’. Since
they can not move, they do not experience the maximum velocity and carry little information
about it. In the influx limited regimes, the overall density is smaller and does not influence the
velocity that much. This can be seen from the depicted trajectories in Figure 8. Note that we
do not observe similar problems using time dependent density profiles, since most trajectories
experience significant velocities before the flow has equilibrated.

Influence of the prior variance The influence of the prior variance is depicted in Figure 9
for the time dependent case and in Figure 10 for steady states. Again, we test this for an
outflux (left panels) and influx (right panels) limited regimes. As before, we observe that the
estimate for vmax is accurate for both regimes in the time dependent case, except when the prior
variance is too small, while for the steady state case we observe again that the outflux limited
case produces a posterior distribution which mimics the prior.

17
Figure 8: Examples of five trajectories for steady state densities. Left: outflux limited regime.
Right: influx limited regime.

Figure 9: Influence of prior variance for two different regimes using solutions of the time depen-
dent Fokker–Planck equation. Left: outflux limited regime with a = 0.4, b = 0.2. Right: influx
limited regime with a = 0.2, b = 0.4.

Figure 10: Influence of prior variance for two different regimes using solutions of the steady
state equation. Left: outflux limited regime with a = 0.4, b = 0.2. Right: influx limited regime
with a = 0.2, b = 0.4.

Influence of information and the proximity of the parameters a and b Next we

discuss how the amount of information, that is the number of trajectories, and the ratio of the

18
parameters a and b influence the quality of the estimates. The ratio of a and b determines
the location of the boundary layer and the density range experienced by the trajectories. For
example, the pair a = 0.2, b = 0.4 gives us a steady state profile with values for ρs varying from
[0.2, 0.5], while a = 0.1, b = 0.15 gives values in [0.075, 0.65].
We will show that the number of trajectories J plays a large influence when a ≈ b, while this
is not the case when these parameters differ by a “large” amount. In this set of experiments only
we use the same value of Σ to both generate the data and in the inference method employed:
σ1 = σ2 = σ0 = 0.05. Figures 11 and 12 illustrate the role of J and of the relative size of a and
b. In each of the panels, we vary the total number of trajectories, that is J = 5, 10, 15, 20, and
study the effect in two influx limited and regimes.

Figure 11: Influence of the number of trajectories on the estimate for vmax for influx limited
regimes and in the time dependent case. Left: a = 0.2, b = 0.4. Right: a = 0.1, b = 0.15.

We observe that, especially in the regime when the ratio between a and b is further from 1,
the estimates become better as the number of trajectories increases, both in the estimate itself
(mean of the posterior distribution) and the variance of the posterior distribution. A similar
behavior is observed for the steady state case, which is presented in Figure 12 below.

Figure 12: Influence of the number of trajectories on the estimate for vmax for influx limited
regimes, in the steady state regime. Left: a = 0.2, b = 0.4. Right: a = 0.1, b = 0.15.

19
 Again, the estimate becomes closer to the true value, and with a smaller variance, as the
number of trajectories increases. Furthermore, the regime where a/b is further from one has
better estimates. This suggests that when planning experiments, these are preferable regimes
to consider. We point out that Figures 11 and 12 are zoomed in so that the effect of the number
of trajectories is clearly seen and, in particular, the prior is supported on a much larger length
scale than that displayed, demonstrating that it has been forgotten.

Time dependent vs steady state regimes We have already pointed out that in a steady
state situation, outflux limited regimes do not carry information about vmax due to the trajec-
tories getting stuck. However, an important thing to notice is that in the influx limited and
maximal current regimes the estimates for vmax are always accurate, and more importantly, less
uncertain. This can be observed in, e.g., Figures 6 and 7.

6.3 Corridor with a bottleneck

We conclude with a more realistic example, in which pedestrians move through a corridor
with a bottleneck. The corridor has length L = 3 and a maximum width of 2` = 0.5, which
reduces to 2w = 0.1 inside the bottleneck. The exit boundary is now split into a rigid wall
(Neumann/reflecting boundary conditions) and a door (Robin/partially reflecting boundary
conditions) which has width 0.3. We assume that diffusion in the vertical direction is smaller
due to limited space inside the bottleneck. In particular we set σ1 = 0.05 and σ2 = 0.03. In
this example the vector e1 in (3) is replaced by the negative gradient of a potential φ. This
potential is computed from the eikonal equation with suitable boundary conditions, see [35].
The function φ corresponds to the shortest distance to the exit, and its negative gradient to the
optimal trajectory to navigate towards a door and/or through a bottleneck, see Figure 13.

Figure 13: Solution of the eikonal equation in the corridor with a bottleneck described in this
section. The arrows depict the direction of the negative gradient, which is the direction with
which individuals move.

We generate trajectories using a time step ∆tSDE = 10−3 and set the final time T = 1. In
the PDE solver the time step was set to ∆tP DE = 4×10−3 . The true value of vmax was, as before
vmax = 1.5, and we used two of the previous regimes: a = 0.2, b = 0.4 and a = 0.4, b = 0.2.
In Figure 14 we present our numerical results. The figure depicts the true value of vmax
(dashed gray line), prior distribution (full gray line), posterior distribution (full blue line) and
the Nelder–Mead algorithm iterates (orange circles). We observe that both regimes produce
very good estimates for vmax . This demonstrates that the framework we have developed can be
applied to more realistic setups; we intend to study the impact of parameters and geometry in
more detail in future work.

20
Figure 14: Results of the parameter estimation methodology in a corridor with a bottleneck.
Left: outflux limited regime with a = 0.4, b = 0.2. Right: influx limited regime with a =
0.2, b = 0.4.

7 Conclusions/Discussion
We have studied a macroscopic model for a unidirectional flow of pedestrians in a corridor.
The evolution of the pedestrian density is given by a nonlinear Fokker–Planck equation, whose
coefficients depend on the so-called fundamental diagram. We formulated and analyzed the
problem of estimating two parameters of interest in the fundamental diagram using individual
trajectories. We assume that these trajectories are realizations of a generalized McKean–Vlasov
equation. This identification problem was solved using derivative-free methodologies. We have
shown that the first parameter – the maximum pedestrian density – cannot be estimated from
the model considered. The second characteristic quantity – the maximum speed – can be
accurately learned from a variety of inflow and outflow conditions, both in time dependent and
steady state settings. We have also seen that boundary conditions play an important role. We
believe that the proposed framework may help to understand their impact in estimation, as well
as experimental design.
We also believe that the developed framework provides the basis for future developments in
applied mathematics as well as transportation research. In particular, the next steps include the
identification of parameters in different forms of the fundamental diagram, or the application
of nonparametric estimation techniques to learn its functional form. Furthermore, we want
to use pedestrian trajectory data, which requires the framework to be generalized to noisy
observations.
These future developments will contribute to the validation of the fundamental diagram
adopted in many models in the transportation literature. It will also justify its use in certain
parameter regimes for microscopic pedestrian models.
Acknolwedgements The work of S.G. and A.S. was supported by the EPSRC Programme
Grant EQUIP. The work of M.-T.W. and AS was supported by a Royal Society international
collaboration grant. M.-T. W. acknowledges partial support from the Austrian Academy of
Sciences via the New Frontier’s grant NST-001. SG is grateful to Imperial College London for
use of computer facilities. The authors are grateful to Grigorios Pavliotis for helpful discussions.

21
References
[1] A. Aggarwal, R.M. Colombo, and P. Goatin. Nonlocal systems of conservation laws in
several space dimensions. SIAM Journal on Numerical Analysis, 53(2):963–983, 2015.

[2] A. Apte, C.K.R.T. Jones, and A.M. Stuart. A Bayesian approach to Lagrangian data
assimilation. Tellus A: Dynamic Meteorology and Oceanography, 60(2):336–347, 2008.

[3] U. M. Ascher, S. J. Ruuth, and B.T.R. Wetton. Implicit-explicit methods for time-
dependent partial differential equations. SIAM Journal on Numerical Analysis, 32(3):797–
823, 1995.

[4] N.W.F. Bode, M. Chraibi, and S. Holl. The emergence of macroscopic interactions between
intersecting pedestrian streams. Transportation Research Part B: Methodological, 119:197–
210, 2019.

[5] M. Boltes, A. Seyfried, B. Steffen, and A. Schadschneider. Automatic extraction of pedes-

trian trajectories from video recordings. In Pedestrian and Evacuation Dynamics 2008,
pages 43–54. Springer, 2010.

[6] M. Burger and J.-F. Pietschmann. Flow characteristics in a crowded transport model.
Nonlinearity, 29(11):3528, 2016.

[7] U. Chattaraj, A. Seyfried, and P. Chakroborty. Comparison of pedestrian fundamental

diagram across cultures. Advances in complex systems, 12(03):393–405, 2009.

[8] R.M. Colombo and M. Lécureux-Mercier. Nonlocal crowd dynamics models for several
populations. Acta Mathematica Scientia, 32(1):177–196, 2012.

[9] A. Conn, K. Scheinberg, and L. N. Vicente. Introduction to derivative-free optimization.

Society for Industrial and Applied Mathematics, 2009.

[10] A. Corbetta, C.-M. Lee, R. Benzi, A. Muntean, and F. Toschi. Fluctuations around mean
walking behaviors in diluted pedestrian flows. Phys Rev E, 95:032316, 2017.

[11] S. L. Cotter, G. O. Roberts, A. M. Stuart, and D. White. MCMC methods for functions:
modifying old algorithms to make them faster. Statistical Science, 28(3):424–446, 2013.

[12] E. Cristiani, B. Piccoli, and A. Tosin. Multiscale modeling of pedestrian dynamics, vol-
ume 12. Springer, 2014.

[13] D. A Dawson. Critical dynamics and fluctuations for a mean-field model of cooperative
behavior. Journal of Statistical Physics, 31(1):29–85, 1983.

[14] K.D. Elworthy. Stochastic differential equations on manifolds, volume 70. Cambridge
University Press, 1982.

[15] R. Erban and J. Chapman. Reactive boundary conditions for stochastic simulations of
reaction-diffusion processes. Physical Biology, 4(1):16, 2007.

[16] S. Ermakov, V. V. Nekrutkin, and A. S. Sipin. Random processes for classical equations of
mathematical physics, volume 34. Springer Science & Business Media, 1989.

[17] J. Garnier, G. Papanicolaou, and T.-W. Yang. Mean field model for collective motion
bistability. arXiv preprint arXiv:1611.02194, 2016.

22
[18] J. Garnier, G. Papanicolaou, and T.-W. Yang. Consensus convergence with stochastic
effects. Vietnam Journal of Mathematics, 45(1-2):51–75, 2017.

[19] J. Gärtner. On the McKean-Vlasov limit for interacting diffusions. Mathematische

Nachrichten, 137(1):197–248, 1988.

[20] K. Giesecke, G. Schwenkler, and J. Sirignano. Inference for large financial systems. 2017.
Boston University Questrom School of Business Research.

[21] P. Goatin and M. Mimault. A mixed system modeling two-directional pedestrian flows.
Mathematical biosciences and engineering, 12(2):375–392, 2015.

[22] M. Hairer, A. M. Stuart, J. Voss, and P. Wiberg. Analysis of SPDEs arising in path
sampling. part 2: The Nonlinear Case. Ann. Appl. Prob, 17(5):1657–1706, 2007.

[23] D. Helbing and P. Molnar. Social force model for pedestrian dynamics. Phys Rev E,
51(5):4282–4286, 1995.

[24] D.J. Higham. An algorithmic introduction to numerical simulation of stochastic differential

equations. SIAM Review, 43(3):525–546, 2001.

[25] R.L. Hughes. The flow of human crowds. Annual review of fluid mechanics, 35(1):169–182,
2003.

[26] S.M. Iacus. Simulation and inference for stochastic differential equations: with R examples.
Springer Science & Business Media, 2009.

[27] J. P. Kaipio and E. Somersalo. Statistical and Computational Inverse Problems, volume
160. Springer, 2005.

[28] Y.A. Kutoyants. Statistical Inference for Ergodic Diffusion Processes. Spinger, 2004.

[29] J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright. Convergence properties

of the Nelder–Mead simplex method in low dimensions. SIAM Journal on optimization,
9(1):112–147, 1998.

[30] M. J. Lighthill and G. B. Whitham. On kinematic waves. ii. a theory of traffic flow on long
crowded roads. Proceedings of the Royal Society of London. Series A, Mathematical and
Physical Sciences, pages 317–345, 1955.

[31] J. A. Nelder and R. Mead. A simplex method for function minimization. The Computer
Journal, 7(4):308–313, 1965.

[32] K. Oelschlager. A martingale approach to the law of large numbers for weakly interacting
stochastic processes. The Annals of Probability, pages 458–479, 1984.

[33] B. Oksendal. Stochastic Differential Equations: An Introduction with Applications.

Springer, 2000.

[34] G.A. Pavliotis and A.M. Stuart. Parameter estimation for multiscale diffusions, 2007.

[35] J. Qian, Y. Zhang, and H. Zhao. Fast sweeping methods for eikonal equations on triangular
meshes. SIAM Journal on Numerical Analysis, 45(1):83–107, 2007.

[36] B.L.S. Prakasa Rao. Statistical inference for diffusion type processes. Kendall’s Lib. Statist.
8, 1999.

23
[37] G.O. Roberts and O. Stramer. On inference for partially observed nonlinear diffusion
models using the Metropolis–Hastings algorithm. Biometrika, 88(3):603–621, 2001.

[38] J. Schöberl. Netgen/Ngsolve. https://ngsolve.org/, 2010.

[39] A. Seyfried, B. Steffen, W. Klingsch, and M. Boltes. The fundamental diagram of pedes-
trian movement revisited. Journal of Statistical Mechanics: Theory and Experiment,
2005(10):P10002, 2005.

[40] A. Singer, Z. Schuss, A. Osipov, and D. Holcman. Partially reflected diffusion. SIAM J.
Appl. Math., 68(3):844–868, 2008.

[41] A.V. Skorokhod. Stochastic equations for diffusion processes in a bounded region. Theory
Probab. Appl., 6(3):264–274, 1961.

[42] A.J. Wood. A totally asymmetric exclusion process with stochastically mediated entrance
and exit. J. Phys. A: Math. Theor., 42:445002, 2009.

[43] J. Zhang, W. Klingsch, A. Schadschneider, and A. Seyfried. Ordering in bidirectional pedes-

trian flows and its influence on the fundamental diagram. Journal of Statistical Mechanics:
Theory and Experiment, 2012(02):P02002, 2012.

Appendix
7.1 Proof of Theorem 1 and Corollary 1
Equation (5a) is a gradient flow with respect to the entropy functional
Z
E(ρ) = (ρ log ρ + (1 − ρ) log(1 − ρ) − ρV )dx (20)
Ω

δE
with V (x) = x1 . The corresponding entropy variable u = δρ = log ρ − log(1 − ρ) − V allows us
to write (5a) as

∂t ρ = div(m(ρ)∇u)

with a nonlinear mobility m(ρ) = ρ(1 − ρ). This gradient flow structure provides the necessary
a-priori estimates to prove global in time existence of solutions.

Proof of Theorem 1. Let N ∈ N and consider a discretization of (0, T ] into subintervals (0, T ] =
T
∪Nk=1 ((k − 1)τ, kτ ] with time steps τ = N . We look for a weak solution ρ : Ω × [0, T ] → S of
(5a) in the sense of (8). The proof is based on the implicit Euler discretization, which gives us
a recursive sequence of elliptic problems. We consider its regularized version
ρk − ρk−1
= div(m(ρk )∇uk ) + τ ∆uk . (21)
τ
Existence of at least one weak solution ρ ∈ H 1 (Ω) ∩ L∞ (Ω) with 0 ≤ ρ ≤ 1 to (21) follows from
Theorem 3.5 in [6], if assumptions (A1)-(A3) are satisfied. Note that the transformation from
ρ to the entropy variables u is one to one and given by

eu+V
ρ= . (22)
1 + eu+V

24
Therefore solutions ρ lie automatically in the set S. Furthermore the entropy density
h(ρ) = ρ log ρ − (1 − ρ) log(1 − ρ) − ρV
ρ
is strictly convex for ρ ∈ S 0 since h0 (ρ) = log 1−ρ − V and h00 (ρ) = ρ1 + 1−ρ
1
. Since h is convex
0
we have that h(q1 ) − h(q2 ) ≤ h (q1 )(q1 − q2 ) and therefore for q1 = ρk and q2 = ρk−1 that
h(ρk ) − h(ρk−1 ) ≤ h0 (ρk )(ρk − ρk−1 ). (23)
Entropy dissipation We consider the weak formulation of (21) to obtain
Z Z Z
1 T
(ρ − ρk−1 )ϕ dx + ∇ϕ m(ρk )∇uk dx−a (1 − ρk )ϕ ds+
τ Ω Ω
Z Γin Z
b ρk ϕ ds + τ ∇ϕ∇uk = 0,
Γout Ω

for ϕ ∈ H 1 (Ω) and use that ρk = h0−1 (uk ). Since uk = h0 (ρk ) we can rewrite (23) as
Z Z
1 1
(ρk − ρk−1 )uk dx ≥ [h(ρk ) − h(ρk−1 )]dx. (24)
τ Ω τ Ω
Using (24) and choosing test functions ϕ = uk gives
Z Z Z
T
h(ρk )dx + τ ∇uk m(ρk )∇uk dx − aτ (1 − ρk )uk ds+
Ω Ω Γin
Z Z Z
+bτ ρk uk ds + τ 2 ∇uk ∇uk ≤ h(ρk−1 )dx.
Γout Ω Ω

Next we prove the following a-priori estimate:

Lemma 1. Let ρ ∈ L2 (Ω) and ρ ∈ S 0 a.e. be such that u = h0 (ρ) ∈ H 1 (Ω). Then there exist
constants such that
Z Z Z Z
∇uT m(ρ)∇u − a (1 − ρ)u ds + b ρu ds ≥ |∇ρ|2 dx − C
Ω Γin Γout Ω

Proof. For the first term we obtain

Z Z
∇uT m(ρ)∇u dx = ρ(1 − ρ)|∇u|2 dx.
Ω Ω
∇ρ
Since ∇u = − ∇V we have
ρ(1−ρ)

|∇ρ|2
Z Z Z Z
∇uT m(ρ)∇udx = dx − 2 ∇V · ∇ρ + ρ(1 − ρ)|∇V |2 dx =
Ω Ω ρ(1 − ρ) Ω Ω
|∇ρ|2
Z
≥ − ρ(1 − ρ)|∇V |2 dx
Ω 2ρ(1 − ρ)
Z Z
1
≥ 2 |∇ρ|2 dx − |∇V |2 dx.
Ω 4 Ω

The estimates follow from Cauchy’s inequality and the fact that ρ(1 − ρ) ≤ 14 . The inflow
boundary term can be rewritten as
Z Z
− (1 − ρ)uds = − (1 − ρ)(log ρ − log(1 − ρ) − V )ds
Γin Γin
Z  
1−ρ
Z
= (1 − ρ) log + 2ρ − 1 ds + [(1 − ρ)V − 2ρ + 1] ds.
Γin ρ Γin

25
The first term on the right hand side is a Kullback–Leibler distance and therefore non-negative.
Since V ∈ H 1 (Ω) we know that its trace satisfies V |∂Ω ∈ L2 (∂Ω) and since ρ ∈ S the second
term is bounded. We can use a similar argument for the outflow term, which can be written as
Z Z
− ρu ds = − ρ(log ρ − log(1 − ρ) − V ) ds =
Γout Γout
 
1−ρ
Z Z
= ρ log + 2ρ − 1 ds + [ρV − 2ρ + 1] ds.
Γout ρ Γout

Again we have a non-negative and a bounded term on the right hand side.

The dissipation inequality and the recursion yields

Z k Z
X n
X Z
h(ρk )dx + τ |∇ρk |2 dx + τ 2 R(uk , uk ) ≤ h(ρ0 )dx + T C. (25)
Ω j=1 Ω k=1 Ω

R
where R(u, ϕ) = Ω ∇u∇ϕdx. This discrete entropy dissipation relation allows us to pass to
the limit τ → 0.

The limit τ → 0: Let ρk denote a sequence of solutions to (21). We define ρτ (x, t) = ρk (x) for
x ∈ Ω and t ∈ ((k − 1)τ, kτ ] Then ρτ solves the following problem where στ denotes the shift
operator, that is (στ ρτ )(x, t) = ρτ (x, t − τ ) for τ ≤ t ≤ T , and
Z TZ  
1
(ρτ − στ ρτ )ϕ + ∇ρτ ∇Φ − ρτ (1 − ρτ )∇ϕ dx
0 Ω τ
Z TZ Z (26)
−a (1 − ρτ )ϕds + b ρτ ϕds = 0.
0 Γin Γout

for test functions ϕ ∈ L2 (0, T ; H 1 (Ω)). Then the entropy dissipation inequality becomes
Z Z T Z Z T Z
2
h(ρτ (T ))dx + |∇ρτ | dx + τ R(uk , uk ) ≤ h(ρ0 )dx + T C. (27)
Ω 0 Ω 0 Ω

which gives the a-priori estimate kρτ kL2 (0,T ;H 1 (Ω)) ≤ K. Using the a-priori estimates we use
(26) to obtain
Z T Z
hρτ − στ ρ, ϕidt ≤ KkϕkL2 (0,T ;H 1 (Ω)) .
0 Ω

We know that ρτ is in L2 (0, T ; H 1 (Ω)) and τ1 hρτ − στ ρτ i in L2 (0, T ; H 1 (Ω)∗ ). Therefore we can
use Aubin’s lemma to conclude the existence of a subsequence, also denoted by ρτ such that for
τ →0

ρτ → ρ strongly in L2 (0, T ; L2 (Ω)). (28)

Finally we check that all terms in (26) converge to the right limit as τ → 0. Because of (28) we
know that

1 − ρτ → 1 − ρ strongly in L2 (0, T, L2 (Ω)).

Since V ∈ H 1 (Ω) we can pass to the limit in the boundary terms as well. This concludes the
existence proof.

26
Corollary 1. Let all assumptions of Theorem 1 be satisfied. Then every weak solution ρ ∈
L2 (0, T ; H 1 (Ω)) is also in C(Ω̄ × (0, T )).

Proof of Corollary (1). The assertion follows from a bootstrap argument. Since ρ ∈ L2 (0, T ; H 1 (Ω))
we can rewrite equation (5a) as

∂t ρ(x, t) − div(Σ∇ρ) = − div(h) in Ω × (0, T )

ρ(x, t) = g(x) in Γ × (0, T )

where h = ρF (ρ) ∈ L2 (0, T ; L2 (Ω)) and g ∈ L2 (0, T ; L2 (Γ)) then

kρkL∞ (0,T ;L2 (Ω)) ≤ khkL2 (0,T ;L2 (Ω)) + kgkL2 (0,T ;L2 (Γ)) .

Since ∆V = 0 we can rewrite (5a) as

∂t ρ − div(Σ∇ρ) + v · ∇ρ = 0 in Ω × (0, T )
(∇ρ + vρ) · n = g in Γ × (0, T )
ρ(x, 0) = 0 in Ω,

with g ∈ L2 (0, T ; Lq (Γ)). This equation has a unique solution ρ ∈ Lq (0, T ; W 2,q (Ω))∩W 1,q (0, T ; Lq (Ω)).
From Sobolev embeddings we deduce that ρ ∈ C 0 (Ω̄ × (0, T )).

27