{'fraction_non_alphanumeric': 0.031016882606988615, 'fraction_numerical': 0.009567499431735994, 'mean_word_length': 4.643615160349854, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Parameter estimation for macroscopic pedestrian dynamics models from microscopic data

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


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]. Individual 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.  Then more general interactions, as proposed in [1,8,21,25], can be considered. In this paper 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 diagram. 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 pedestrian 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 conditions, 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 calculated 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 v max 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 v max can be found in the literature, ranging from 3.8 to 10 pedestrians per m 2 for ρ max and 0.98 m/s to 1.5 m/s for v max . 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 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 fundamental 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 v max and the maximum density ρ max . Whilst there is some work concerning 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 v max 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 trajectories, 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 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 v max 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.  Figure 1a. On the macroscopic level the evolution of the density ρ(x, t) can be described by a nonlinear Fokker-Planck (FP) equation

The microscopic and the macroscopic models
where the diffusion matrix is of the form Σ = diag(σ 2 1 , σ 2 2 ). 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 → R 2 of the form This convective field depends on density only, and will be evaluated pointwise to drive individual trajectories. 1 We choose the velocity as 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 v max 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 where 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 ≤ v max . 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 The system is supplemented with the boundary conditions Here the scaled flux is given byj = −Σ∇ρ + v max (1 −ρ)ρe 1 . 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 v max and ρ max in f (ρ) (3) using individual trajectories. However 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: Remark 1. The parameter ρ max can not be identified within our adopted microscopic macroscopic 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 v max . The parameter ρ max thus plays no further part in the paper.
Our focus, then, is on identification of v max . 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 limited, 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 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 [13,19] for theory relating to such PDEs in related settings.

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.

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 v max to the inflow and outflow (1) Influx limited phase in the case a < b and min (a, b) < vmax 2 . We observe an asymptotically low density (ρ < 1 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 (ρ > 1 2 ) and a boundary layer at the entrance.
Here we observe an asymptotic density with ρ ≈ 1 2 in most of the domain when σ is small and boundary layers at the entrance and exit. In the case a = b = vmax 2 , the unique steady state is given by ρ(x) = 1 2 . Figure 2 illustrates the different regimes in 1D in the case of a corridor of length L = 3, v max = 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 v max = 1. Furthermore, we consider the more general case in which e 1 is replaced by ∇V in (2). We make the following assumptions: (A2) 0 ≤ a, b ≤ 1, and min{a, b} < 1.
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.
We seek weak solutions in the set and refer to its interior as S 0 . Then the weak solution ρ : for all test functions ϕ ∈ L 2 (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 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.

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 R 2 . 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).

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 v max . However, as we have seen, determination of ρ max is not possible using the current setting and only determination of v max is feasible.
We ] 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 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 best value of the parameter v, given an observation of a trajectory X, corresponds to minimizing the function 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 where the function Ψ is defined as Note that the last term in Ψ is an Itô stochastic integral. It is preferable to perform the parameter estimation using Ψ instead of Φ, since the latter is infinite almost surely, whilst the former is finite almost surely, and since T 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 v max by minimizing 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 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 where π 0 is the prior Lebesgue density for v. Using Girsanov's Theorem [14,33], we see that 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 which is exactly exp(−J (v; X))1(v > 0). We can either sample from the probability distribution given by (16) (the Bayesian approach) 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 approaches 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 X j (t), j = 1, . . . , J. We assume that each trajectory X j starts at t j 0 with X j (t j 0 ) = X 0 j and exits at time Then we can sum the misfit of all trajectories and obtain 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 constant 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].

Computational methods
In this section we discuss the discretization of the nonlinear PDE (5) and SDE (6) before continuing 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 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.

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 discretization. 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].

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 X k ≈ X(t k ) := X(k∆t) and updating the individual position as 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 ρ(X k , t k ) by the steady state solution ρ s (X k ).
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.

Derivative-free optimization -the Nelder-Mead algorithm
As detailed in Section 4, we will use derivative-free methods for the minimization of (12), corresponding 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 unconstrained 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 Wait to get in

Robin
Partially reflecting 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 the following time step, which is realized with probability p. Here we use the rescalings P in = π∆t/(2σ 2 1 )a(1 − ρ(0, x 2 )) and P out = π∆t/σ 2 1 bρ(L, x 2 ) according to [15]. 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

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

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. Propose N (0, c).

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 v max 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 v max 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 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 trajectories 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 ∆t SDE = 10 −3 . We assume that the true value for the maximum velocity is v max = 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 pedestrian flow has equilibrated. We consider five sets of parameters which span the three steady state density regimes: Since the variance coefficients might not be known in practice, we will estimate the parameters 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].

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 v max 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 v max . 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 v max 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 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 approximately 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 v max 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.

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.  Influence of the initial guess and prior mean We assume that v max 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  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 v max 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.   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. 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. 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 v max due to the trajectories getting stuck. However, an important thing to notice is that in the influx limited and maximal current regimes the estimates for v max are always accurate, and more importantly, less uncertain. This can be observed in, e.g., Figures 6 and 7.

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 e 1 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 ∆t SDE = 10 −3 and set the final time T = 1. In the PDE solver the time step was set to ∆t P DE = 4×10 −3 . The true value of v max was, as before v max = 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 v max (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 v max . 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.

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.