Parameter estimation for pedestrian dynamics models Susana Gomes - - PowerPoint PPT Presentation
Parameter estimation for pedestrian dynamics models Susana Gomes Mathematics Institute, University of Warwick Joint work with A.M. Stuart (Caltech) and M.-T. Wolfram (Warwick) Crowds: Models and Control CIRM, Marseille, 5 th of June 2019
Mathematics Institute, University of Warwick Joint work with A.M. Stuart (Caltech) and M.-T. Wolfram (Warwick)
Crowds: Models and Control CIRM, Marseille, 5th of June 2019
hello world We use a simple model exploring a coupling between
pedestrians (data from experiments1), and
pedestrians to identify and estimate quantities of interest for the dynamics of the crowd.
1BaSiGo experiments (Forschungszentrum Jülich) 2Figure from Chattaraj, Seyfried, Chakroborty, Advances in Complex Systems, 2009. 1 / 17
hello world We use a simple model exploring a coupling between
pedestrians (data from experiments1), and
pedestrians to identify and estimate quantities of interest for the dynamics of the crowd. In particular, we would like to
diagram2, or
diagram, and
diagram for the time dependent evolution of crowds. hello world hello world hello world
1BaSiGo experiments (Forschungszentrum Jülich) 2Figure from Chattaraj, Seyfried, Chakroborty, Advances in Complex Systems, 2009. 1 / 17
(0, −ℓ) (L, −ℓ) (0, ℓ) (L, ℓ) ΓN ΓN Γin Γout F(ρ(x, t))
Individual trajectories are realizations of a generalized McKean–Vlasov equation dXt = F(ρ(Xt, t)) dt + √ 2Σ dWt.
ρmax
fundamental diagram). The microscopic equation depends on the density of the process.
2 / 17
The corresponding Fokker–Planck equation is ρt = ∇ ·
ρ ρmax
(0, −ℓ) (L, −ℓ) (0, ℓ) (L, ℓ) ΓN ΓN Γin Γout F(ρ(x, t))
with boundary conditions: (−Σ∇ρ + F(ρ)ρ) · n = −a (ρmax − ρ) , for x ∈ Γin, (−Σ∇ρ + F(ρ)ρ) · n = bρ, for x ∈ Γout, (−Σ∇ρ + F(ρ)ρ) · n = 0, for x ∈ ΓN. where n is the unit normal, a is the inflow rate and b the outflow rate.
3 / 17
By defining ˜ ρ =
ρ ρmax , we obtain a rescaled SDE-PDE system:
dXt = vmax (1 − ˜ ρ) dt + √ 2Σ dWt, ˜ ρt = ∇ · (Σ∇˜ ρ − vmax˜ ρ(1 − ˜ ρ)) , with the boundary conditions ˜ j · n = −a (1 − ˜ ρ) , for x ∈ Γin, ˜ j · n = b˜ ρ, for x ∈ Γout, ˜ j · n = 0, for x ∈ ΓN. and where ˜ j = (−Σ∇˜ ρ + vmax(1 − ˜ ρ)˜ ρ).
4 / 17
By defining ˜ ρ =
ρ ρmax , we obtain a rescaled SDE-PDE system:
dXt = vmax (1 − ˜ ρ) dt + √ 2Σ dWt, ˜ ρt = ∇ · (Σ∇˜ ρ − vmax˜ ρ(1 − ˜ ρ)) , with the boundary conditions ˜ j · n = −a (1 − ˜ ρ) , for x ∈ Γin, ˜ j · n = b˜ ρ, for x ∈ Γout, ˜ j · n = 0, for x ∈ ΓN. and where ˜ j = (−Σ∇˜ ρ + vmax(1 − ˜ ρ)˜ ρ). ρmax does not influence the dynamics in the rescaled SDE-PDE system! ⇒ we cannot estimate ρmax using this model.
4 / 17
Steady solutions of the FP equation are influenced by a and b.3
3Burger and Pietschmann, Nonlinearity, 2016. 5 / 17
Steady solutions of the FP equation are influenced by a and b.3
Properties of solution Solution in 1D
⋆ Maximal current if min(a, b) ≥ vmax
2 ;
asymptotic density ˜ ρ = 1
2 and boundary
layers at both entrance and exit.
3Burger and Pietschmann, Nonlinearity, 2016. 5 / 17
Steady solutions of the FP equation are influenced by a and b.3
Properties of solution Solution in 1D
⋆ Maximal current if min(a, b) ≥ vmax
2 ;
asymptotic density ˜ ρ = 1
2 and boundary
layers at both entrance and exit.
⋆ Influx limited if a < b and a < vmax
2 ;
asymptotic low ˜ ρ < 1
2 and boundary layer at
the exit.
3Burger and Pietschmann, Nonlinearity, 2016. 5 / 17
Steady solutions of the FP equation are influenced by a and b.3
Properties of solution Solution in 1D
⋆ Maximal current if min(a, b) ≥ vmax
2 ;
asymptotic density ˜ ρ = 1
2 and boundary
layers at both entrance and exit.
⋆ Influx limited if a < b and a < vmax
2 ;
asymptotic low ˜ ρ < 1
2 and boundary layer at
the exit.
⋆ Outflux limited if b < a and b < vmax
2 .
asymptotic high ˜ ρ > 1
2 and boundary layer at
the entrance.
3Burger and Pietschmann, Nonlinearity, 2016. 5 / 17
The Fokker–Planck equation has a gradient flow structure with respect to the entropy E(ρ) =
(ρ log ρ + (1 − ρ) log(1 − ρ) − ρV dx, withV = x. ❘
6 / 17
The Fokker–Planck equation has a gradient flow structure with respect to the entropy E(ρ) =
(ρ log ρ + (1 − ρ) log(1 − ρ) − ρV dx, withV = x. Using this gradient structure we can
existence of weak solutions.
solutions ρ ∈ C(¯ Ω × [0, T]). ❘
6 / 17
The Fokker–Planck equation has a gradient flow structure with respect to the entropy E(ρ) =
(ρ log ρ + (1 − ρ) log(1 − ρ) − ρV dx, withV = x. Using this gradient structure we can
existence of weak solutions.
solutions ρ ∈ C(¯ Ω × [0, T]). Further extensions on these results would allow us to prove existence and uniqueness of strong solutions for the SDE if it was posed in ❘2.
6 / 17
Goal: Estimate vmax given a set of trajectories.
7 / 17
Goal: Estimate vmax given a set of trajectories. Finding the best value of v given an observation of a trajectory X corresponds to minimizing the function Φ(v; X) = 1 4 T | ˙ X − F (ρ(X(t), t); v) |2
Σ,
7 / 17
Goal: Estimate vmax given a set of trajectories. Finding the best value of v given an observation of a trajectory X corresponds to minimizing the function Φ(v; X) = 1 4 T | ˙ X − F (ρ(X(t), t); v) |2
Σ,
7 / 17
Goal: Estimate vmax given a set of trajectories. Finding the best value of v given an observation of a trajectory X corresponds to minimizing the function Φ(v; X) = 1 4 T | ˙ X − F (ρ(X(t), t); v) |2
Σ,
Φ(v; X) = Ψ(v; X) + 1 4 T | ˙ X|2 dt, where the function Ψ is defined as Ψ(v; X) = 1 4 T
Σdt − 2F (ρ(X(t), t); v) , dX(t)Σ
Ψ(v, X) is the log-likelihood function of this process... ... which we can obtain it due to the PDE-SDE coupling in our model.
7 / 17
We will include prior information and consider instead the functional J (v; Xt) = Ψ(v; Xt) + 1 2|v − m|2
c. ✶
P P
8 / 17
We will include prior information and consider instead the functional J (v; Xt) = Ψ(v; Xt) + 1 2|v − m|2
c.
Using Bayes’ and Girsanov’s theorems, we can show that the function e−J (v;Xt )✶(v>0), appropriately normalised, is the posterior distribution P(v|X) of v given the
P
8 / 17
We will include prior information and consider instead the functional J (v; Xt) = Ψ(v; Xt) + 1 2|v − m|2
c.
Using Bayes’ and Girsanov’s theorems, we can show that the function e−J (v;Xt )✶(v>0), appropriately normalised, is the posterior distribution P(v|X) of v given the
8 / 17
We cannot take T → ∞ in our observations. Therefore we use N = 20 trajectories
t
i=1,...,20
t∈[0,T]
to compute our estimate of vmax.
⋆ Minimise J directly; ⋆ Provides the maximum a posteriori estimator; ⋆ Needs extra computations to quantify uncertainty; ⋆ Use derivative free algorithm: Nelder–Mead algorithm4.
⋆ Sample from the probability distribution P(v|Xi
t);
⋆ Provides posterior distribution; ⋆ Enables uncertainty quantification; ⋆ Use derivative free MCMC methodology: pCN algorithm5.
4Nelder and Mead, The Computer Journal, 1965. 5Cotter et al., Statistical Science, 2013. 9 / 17
We solve the steady and time-dependent Fokker–Planck equation efficiently using the Netgen/NgSolve interface.6
⋆ 4th order IMEX scheme ⋆ Explicit: nonlinear convective oper- ator – upwind discontinuous Galerkin discretization. ⋆ Implicit: linear diffusion – hybrid dis- continuous Galerkin discretization.
⋆ Damped Newton solver ⋆ H1 conforming finite elements of or- der 1 ⋆ Initial guess depends on the magni- tude and relation between the parame- ters a and b.
6https://ngsolve.org/ 10 / 17
To advance the microscopic model in time, we use the Euler–Maruyama scheme and use PDE solution in the SDE drift: Xi
t+∆t = Xi t + F(ρ(Xi t, t))∆t +
√ 2∆tΣξ, ξ ∼ N(0, 1),
7Erban and Chapman, Physical Biology, 2007. 11 / 17
To advance the microscopic model in time, we use the Euler–Maruyama scheme and use PDE solution in the SDE drift: Xi
t+∆t = Xi t + F(ρ(Xi t, t))∆t +
√ 2∆tΣξ, ξ ∼ N(0, 1), together with a careful application of the boundary conditions.7
Boundary PDE SDE Implementation ΓN Neumann Reflecting
p = 1
Γin (from outside) Robin Wait to get in
p = Pin p = 1 − Pin
Γin (from inside) Robin Partially reflecting
p = 1 − Pin p = Pin back to Γin
Γout Robin Partially reflecting
p = 1 − Pout p = Pout leave domain 7Erban and Chapman, Physical Biology, 2007. 11 / 17
(maximal current) vmax = 1.5, a = 0.9, b = 0.975 and σ1 = σ2 = 0.05. Left: Influence of the PDE time step on the posterior distribution. Right: Comparison between Bayesian and optimization approach for ∆t = 0.005. MAP estimator obtained by derivative-free optimization coincides with posterior mean for all test cases – suggests posterior distribution is Gaussian.
12 / 17
Influence of prior mean and initial guess for outflux limited regime (a = 0.45, b = 0.4). Estimates are independent of
13 / 17
Influence of prior mean and initial guess for (left) outflux limited regime, (middle) influx limited regime and (right) maximal current regime. Outflux limited regimes (left) produce estimates which mimic the prior distribution! Influx limited and maximal current regimes produce reliable and less uncertain esti- mates.
14 / 17
Influence of number of trajectories for two influx limited regimes Using time dependent solution Using steady state solution
15 / 17
We presented a parameter estimation framework for SDEs which depend on their density.
(experiments) cases.
regimes (trajectories get stuck and there is no information).
16 / 17
We presented a parameter estimation framework for SDEs which depend on their density.
(experiments) cases.
regimes (trajectories get stuck and there is no information).
the corridor, doors in the entrance and exit boundaries, and corridors with bottlenecks.
16 / 17
We presented a parameter estimation framework for SDEs which depend on their density.
(experiments) cases.
regimes (trajectories get stuck and there is no information).
the corridor, doors in the entrance and exit boundaries, and corridors with bottlenecks.
experimental data.
16 / 17
For more details, arXiv:1809.08046v2 susana.gomes@warwick.ac.uk
17 / 17