Data Assimilation for Agent-Based Models Mathematical and Computational Challenges Jonathan Ward
[email protected] i-Sim 2017 Collaborators: Nick Malleson & Andy Evans (Leeds)
Goals Do DA Do ABMs
Don’t do ABMs
Don’t do DA
Goals Do DA Do ABMs
Don’t do ABMs
Let’s talk!
Don’t do DA
Goals Do DA Do ABMs
Don’t do ABMs
Don’t do DA Intro to DA: It’s useful!
Goals Do DA Do ABMs
Don’t do ABMs
ABM challenges. You could help!
Don’t do DA
Goals Do DA
Don’t do DA
Do ABMs
Don’t do ABMs
Enjoy!
Overview
Overview •
What are Agent-Based Models (ABMs)?
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
•
Some comments on DA for ABMs
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
•
Some comments on DA for ABMs
ABMs •
Agents as objects, an environment
•
Complexity, emergence
•
ABMs to explain; multi-agent simulation to solve practical problems; individual Based Models: less/no interaction
•
Informed by “data”
•
Heterogeneity, interaction networks, learning,… (Kitchen sink approach)
ABMs •
Epidemics
•
Voting behaviour
•
Innovation diffusion
•
Statistical physics models
•
Cellular Automata
I
S
I
µ
ABMs •
NetLogo: Ants, DaisyWorld, forest fires…
•
Altruism
•
Co-operation
•
Social influence
•
Segregation
ABMs •
Commercial planning models (multi-modal transport)
http://vision-traffic.ptvgroup.com
•
Others? (Uber, Dunnhumby…)
Simple ABM Example •
Simple example to illustrate how the Ensemble Kalman Filter (EnKF) works
•
Simple enough to know everything about the model dynamics and parameters
•
Use the model to indicate issues that could arise
Model Motivation Footfall counts in Leeds UK city centre:
https://datamillnorth.org
•
10000
-0.5 -1
8000
-1.5 -2
6000
-2.5
z¯w
-3
4000
-3.5 -4
2000
-4.5 -5
0
-5.5
0
12
24
36
48
60
72
84
t (hours)
96
108
120
132
144
156
Dynamic calibration of agent-based models using data assimilation
rsos.royalsocietypublishing.org
Model Motivation Jonathan A. Ward1 , Andrew J. Evans 2 and
Research
Nicolas S. Malleson2
Cite this article: Ward JA, Evans AJ, Malleson NS. 2016 Dynamic calibration of agent-based models using data assimilation. R. Soc. open sci. 3: 150703. http://dx.doi.org/10.1098/rsos.150703
1 School of Mathematics, and 2 School of Geography, University of Leeds,
Leeds LS2 9JT, UK
Footfall counts in Leeds UK city centre:
https://datamillnorth.org
•
A widespread approach to investigating the dynamical behaviour of complex social systems is via agent-based models (ABMs). In this paper, we describe how such models can be dynamically calibrated using the ensemble Kalman filter (EnKF), a standard method of data assimilation. Our goal is twofold. First, we want to present the EnKF in a simple setting for the benefit of ABM practitioners who are unfamiliar with it. Second, we want to illustrate to data assimilation experts the value of using such methods in the context of ABMs of complex social systems and the new challenges these types of model present. We work towards these goals within the context of a simple question of practical value: how many people are there in Leeds (or any other major city) right now? We build a hierarchy of exemplar models that we use to demonstrate how to apply the EnKF and calibrate these using open data of footfall counts in Leeds.
Received: 19 December 2015 Accepted: 11 March 2016
10000
Subject Category: Mathematics Subject Areas: complexity
8000
Keywords: agent-based models, data assimilation, complex systems
2000
12
24
36
48
-2.5
Agent-based models (ABMs) are characterized by a set of rules, typically compiled together in a computer program, which describe the evolution of the model in time from an initial condition. The models are distinguished from others by representing individual actors in the modelled system as individuals, with their own behaviours and histories, often embedded with their own positions within a data environment. Such ‘agents’ may enact behaviour based on the same, or different, parameters and the rules they work to may vary from mathematical equations to adaptable artificial intelligence routines. While agent-based systems can be used for a wide variety of purposes, including systems control, Web-search and robotics [1], ABMs are commonly implemented in science to understand social and environmental systems, either as simple in silico thought experiments or, increasingly, as detailed models of the real world [2]. In the latter case, the aim is not always
4000
0
-1.5
1. Introduction
Author for correspondence: Jonathan A. Ward e-mail:
[email protected] 0
-1
-2
6000
z¯w
-0.5
60
72
84
96
108
120
2016 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.
t (hours)
-3 -3.5 -4 -4.5 -5 -5.5
132
144
156
Simple ABM Example City
Simple ABM Example Agent pool
City
Simple ABM Example Agent pool
City
Simple ABM Example Agent pool
Arrive ↵
City
Simple ABM Example City
Agent pool Depart
Adding Complexity
Malleson, Heppenstall, See, Evans; Environment and Planning B: Planning and Design (2013)
Adding Complexity
Malleson, Heppenstall, See, Evans; Environment and Planning B: Planning and Design (2013)
Adding Complexity
Malleson, Heppenstall, See, Evans; Environment and Planning B: Planning and Design (2013)
Simple ABM Example •
Agent states
•
Population in city
xi (t) 2 {0, 1} X n(t) = xi (t) i
•
Simulate using Gillespie algorithm
Simple ABM Example Number 40 35 30
n
25 20 15 10 5 0 0
10
20
30
40
50
t Time
60
70
80
90
100
Mathematics •
Master equation p˙ n = ↵pn
•
1
Mean (t) =
•
(↵ + n)pn + (n + 1)pn+1
Distribution
↵
✓
↵
+ n0
n
(t) e pn (t) = n!
◆
e
(t)
t
↵ = 50,
= 0.25
250
250
200
nx
150
nx 200
100
50
0
0
5
10
150 50
15
55
t
60
65
70
15
20
t
0.035
1
0.03
0.8
0.025 0.6
p(x)
p
0.02
ah
0.015
0.4
0.01
0.2
0.005 0 140
0 160
180
200
x
220
240
260
n
0
5
10
h
↵ = 1,
= 0.2
70
16
60
14 12
50
10 40 x n
nx
8
30
6
20
4
10
2 0
0 0
10
20
30
40
50
100
110
120
t
130
140
150
t
0.2
1 0.8
0.15 0.6
p(x) p
ah
0.1
0.4 0.2
0.05
0 0
0
5
10
x
15
n
0
5
10
15
h
20
25
30
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
•
Some comments on DA for ABMs
Fitting to Data •
System state not observed directly
•
System state typically much larger than observations
•
Observations have error
•
Goals: •
Fit the model to observations
•
Use the model to make forecasts
Data Assimilation Forecast • Run model(s) forward from estimated state(s) • Compute forecast uncertainty
Analysis • Use observations to adjust model states • Compute Analysis uncertainty
Data Assimilation •
DA routinely used in weather forecasting
•
Well established area with lots of methods:
e.g. Kalman filters, variational methods, particle filters, hidden Markov models •
Typically used when forecasting
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
•
Some comments on DA for ABMs
Ensemble Kalman Filter % Given: initial xa, x->y map H, % obs times t_k, obs z, obs uncertainty R. for t_k in {t_0,…,t_T} % Forecast for i in {1,…,ensemble_size} xf(i)=run_model(xa(i),[t_k, t_k+1]) xf_mean=mean(xf)% Best guess Pf=cov(xf)% Uncertainty % Analysis K=kalman_gain(Pf,H,R) % Solve: (H*Pf*H’+R)K=Pf*H’ for i in {1,…,ensemble_size} xa(i)=xf(i)+K*(y(i)-H*xf(i)) xa_mean=mean(xa) % Best guess Pa=cov(xa) % Uncertainty
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
•
Some comments on DA for ABMs
EnKF - Model Observation times
tk 2 {t0 , t1 , . . . , tT }
•
State vector
x k 2 RN
•
Consider models of the form
•
xk+1 = F (xk ) + ⌫k
Noise
EnKF - Model Observation times
tk 2 {t0 , t1 , . . . , tT }
•
State vector
x k 2 RN
•
Consider models of the form
•
xk+1 = M (xk ; ⇠k )
Random numbers
EnKF - Analysis •
Virtual observations: zk (i) = yk + ⌫k (i)
•
Analysis step:
⌫k (i) ⇠ N (0, R)
Data Forecast
EnKF - Analysis •
Virtual observations: zk (i) = yk + ⌫k (i)
•
Analysis step: Analysis
Forecast
⌫k (i) ⇠ N (0, R)
Data
EnKF - Analysis •
Virtual observations: zk (i) = yk + ⌫k (i)
•
Analysis step: Analysis
Forecast
⌫k (i) ⇠ N (0, R)
Data
Bayes p(x|y) / p(y|x)p(x)
EnKF - Analysis •
Virtual observations: zk (i) = yk + ⌫k (i)
•
Analysis step: Analysis
Forecast
⌫k (i) ⇠ N (0, R)
Data
Bayes p(x|y) / p(y|x)p(x) model stochasticity
EnKF - Analysis •
Virtual observations: zk (i) = yk + ⌫k (i)
•
Analysis step: Analysis
Forecast
⌫k (i) ⇠ N (0, R)
Data
Bayes p(x|y) / p(y|x)p(x) observation error
model stochasticity
EnKF - Analysis •
Virtual observations: zk (i) = yk + ⌫k (i)
•
Analysis step: Analysis
Forecast
⌫k (i) ⇠ N (0, R)
Data
Bayes p(x|y) / p(y|x)p(x) observation error
model stochasticity
Normal distributions
State Estimation
EnKF Example N = 50,
n0 = 10,
↵ = 3,
= 0.1,
R = 25,
Nreal = 10
50 40
x
30 20 10 0 0
10
20
30
40
50
60
70
80
t
RMSEs Steady state
Observations
Forecast
Analysis
8.20
4.32
9.55
4.87
90
100
EnKF Example N = 50,
n0 = 10,
↵ = 3,
= 0.1,
R = 25,
Nreal = 100
50 40
x
30 20 10 0 0
10
20
30
40
50
60
70
80
t
RMSEs Steady state
Observations
Forecast
Analysis
8.20
5.01
8.13
3.39
90
100
EnKF Example N = 50,
n0 = 10,
↵ = 3,
= 0.1,
R = 25,
Nreal = 10
50 40
x
30 20 10 0 0
10
20
30
40
50
60
70
80
t
RMSEs Steady state
Observations
Forecast
Analysis
7.17
5.10
5.62
4.29
90
100
EnKF Example N = 50,
n0 = 10,
↵ = 3,
= 0.1,
R = 25,
Nreal = 100
50 40
x
30 20 10 0 0
10
20
30
40
50
60
70
80
t
RMSEs Steady state
Observations
Forecast
Analysis
7.17
5.11
6.90
4.19
90
100
EnKF Example N = 50,
n0 = 10,
↵ = 3,
= 0.1,
R = 25,
Nreal = 10
50 40
x
30 20 10 0
0
10
20
30
40
50
60
70
80
t
RMSEs Steady state
Observations
Forecast
Analysis
3.30
4.24
3.58
3.11
90
100
Comments •
If measurements accurate, analysis performance will be similar
•
Forecasts will perform better if less noise
•
Is the model a good forecast of the data?
•
(Taking more history doesn’t necessarily help since Markovian)
•
Forecasting involves luck! (Hence errors important)
Exact Analysis t = 50 0.2
y(tk ) = 22.88
Analysis 0.15
p
Data
0.1
Forecast
0.05
0 0
10
20
30
n
40
50
Exact Analysis t = 50 0.2
y(tk ) = 22.88
Analysis 0.15
p
Data
0.1
Forecast
0.05
0 0
10
20
30
n
40
50
Exact Analysis t = 50 0.2
y(tk ) = 22.88
Analysis 0.15
p
Data
0.1
Forecast
0.05
0 0
10
20
30
n
40
50
Exact Analysis 50 40
µn 30 20 10 0 0
10
20
30
40
50
60
70
80
90
100
60
70
80
90
100
t 6 5
σn
4 3 2 0
10
20
30
40
50
t
Observations •
GPS trackers on a sample of agents N X y(t) = xi (t) |S| i2S
40 35 30
n
25 20 15 10 5 0 0
10
20
30
40
50
t
60
70
80
90
100
Error Statistics
•
Variance quadratic in n, but n unknown
•
Resample approximation at each time step
Non-Normality N = 50,
n0 = 40,
↵ = 0.3,
= 0.1,
M = 35,
Nreal = 100
50 40
x
30 20 10 0 0
50
100
150
200
250
300
350
t
RMSEs Steady state
Observations
Forecast
Analysis
Posterior
1.97
1.13
1.94
1.18
0.83
400
Non-Normality N = 50,
n0 = 40,
↵ = 0.3,
= 0.1,
M = 35,
Nreal = 100
10 8
x
6 4 2 0 0
50
100
150
200
250
300
350
t
RMSEs Steady state
Observations
Forecast
Analysis
Posterior
1.97
1.13
1.94
1.18
0.83
400
Non-Normality 60 40
µn
20 0 -20 0
50
100
150
200
250
300
350
400
250
300
350
400
t 3
2
σn 1
0 0
50
100
150
200
t
Parameter Estimation
Parameter Estimation N = 50,
n0 = 10,
↵ = 3.65,
= 0.12,
M = 35,
Nreal = 100
50 40
x
30 20 10 0 0
50
100
150
200
250
300
t
Steady state
Observations
Forecast
Analysis
5.11
1.84
5.26
2.26
350
400
Parameter Estimation 1.5 1
efα
0.5 0 -0.5 -1 -1.5 0
50
100
150
200
250
300
350
400
250
300
350
400
t 1.5 1
eaα
0.5 0 -0.5 -1 -1.5 0
50
100
150
200
t
Parameter Estimation 0.04 0.02
efβ
0 -0.02 -0.04 0
50
100
150
200
250
300
350
400
250
300
350
400
t 0.04 0.02
eaβ
0 -0.02 -0.04 0
50
100
150
200
t
Parameter Estimation 38 36 34
α β 32 30 28 26 0
50
100
150
200
t
250
300
350
400
Parameter Estimation 0
5
-2000
4
-4000 -6000
α
3 -8000 -10000
2
-12000
1
-14000 -16000
0 0
0.1
0.2
0.3
β
0.4
0.5
Other Issues •
Agent heterogeneity
•
Exogenous events & dynamic adaptation
•
Single shot calibration vs dynamic DA
•
Fitting the wrong model
•
Non-Markovian models
•
(Variable rescaling, rounding errors, practical issues)
Overview •
What are Agent-Based Models (ABMs)?
•
What is Data Assimilation (DA)?
•
DA is easy to implement…
•
…but be careful
•
Some comments on DA for ABMs
Calibration Model validation Data assimilation
Data Analysis Statistical properties Benchmark forecasts
Model Mechanisms Avoid `black swans’ via multiple mechanisms How do you combine them?
Mapping Parameter Space Identify dynamical regimes via equation free methods Parameter sensitivity
Model hierarchy Analysis of simple models to inform higher complexity models Dimension reduction
Conclusions •
Understanding model dynamics is important
•
ABMs to forecast?!
•
Gaussian mixtures, mixed discrete/continuous?
•
Identifying when it’s not working
•
Translation of expertise from developed disciplines
•
Can ABMs compare with machine learning?
Thanks!