Crowd Escape Panic Model

Adam Rumpf, 5/1/2018

Introduction

Based on the following paper:

D. Helbing, I. Farkas, and T. Vicsek. Simulating dynamical features of escape panic. Nature, 407:487—490, 2000.

The paper describes a dynamical systems model (see below for a full definition) for crowds of people running to an exit. The model consists of a system of ODEs based on the forces of people acting on each other. This program solves the ODEs using forward Euler’s method.

This demonstration consists mostly of a single Manipulate environment to produce an animation of a crowd attempting to escape a building. People are represented as randomly-sized dots. The color of the dot represents the amount of crush force experienced, with redder color indicating more force, which a planner would ideally like to avoid since this can lead to injury. A time series is also displayed in the lower right corner, with the green series indicating number of successful escapes (which occurs when a dot crosses the threshold of the exit) and the red series indicating total force.

There are also controls to change the layout of the room and the size and behavior of the people. Due to the computational intensity of this program changes are not applied in real time, and are instead applied by pressing the “apply changes” button, at which point the simulation will restart with the new parameters. Before applying the changes, the prospective layout is drawn over the current layout in a gray outline. This button can also be pressed even without applying changes simply to randomly restart the simulation with different starting positions.

ODE Model

Notation

t = time
N = number of pedestrians
(Subscript[m, i]) = mass of pedestrian i
((Subscript[v, i])^0)(t) = desired speed of pedestrian i at time t
((Subscript[(OverVector[e]), i])^0)(t) = desired direction of pedestrian i at time t
(Subscript[(OverVector[v]), i])(t) = actual velocity of pedestrian i at time t
(Subscript[τ, i]) = characteristic time for pedestrian i (time taken for actual velocity to match desired velocity)
W = set of walls
(Subscript[(OverVector[f]), ij]) = interaction force between pedestrians i and j
(Subscript[(OverVector[f]), iW]) = interaction foroce between pedestrian i and wall W
(Subscript[(OverVector[r]), i])(t) = position of pedestrian i at time t
(Subscript[A, i]),(Subscript[B, i]) = constants describing the repulsive interaction force acting on pedestrian i
(Subscript[d, ij]) = distance between pedestrian i and j’s centers of mass
(Subscript[d, iW]) = distance between pedestrain i and wall W
(Subscript[(OverVector[n]), ij])=(((Subscript[n, ij])^1),((Subscript[n, ij])^2)) = normalized vector pointing from pedestrian j to i
(Subscript[(OverVector[n]), iW]) = direction perpendicular to wall W
(Subscript[r, i]) = radius of pedestrian i
(Subscript[r, ij]) = sum of radii of pedestrians i and j
k,κ = large constants that control body force and sliding friction, respectively
(Subscript[(OverVector[t]), ij]) = tangential direction between pedestrians i and j
(Subscript[(OverVector[t]), iW]) = direction tangential to wall W
((Subscript[Δv, ij])^t) = difference in tangential velocities of pedestrians i and j
g(x)={({{0, if x<0}, {x, if x≥0}})=max{0,x} = function to zero out negative arguments

Estimated Parameters

(Subscript[m, i])=80kg
((Subscript[v, i])^0)=1(m/s) , which corresponds to the paper’s “normal” conditions
(Subscript[τ, i])=0.5s
(Subscript[A, i])=2×(10^3)N
(Subscript[B, i])=0.08m
k=1.2×(10^5)(kg/(s^2))
κ=2.4×(10^5)(kg/(m s))
2(Subscript[r, i])~U[0.5m,0.7m]

Model Definition

The change in pedestrian i’s velocity at time t is given by the acceleration equation

(Subscript[m, i])((d(Subscript[(OverVector[v]), i]))/(dt))=(Subscript[m, i])((((Subscript[v, i])^0)(t) ((Subscript[(OverVector[e]), i])^0)(t)-(Subscript[(OverVector[v]), i])(t))/(Subscript[τ, i]))+(∑_(j(≠i)))(Subscript[(OverVector[f]), ij])+(∑_W)(Subscript[(OverVector[f]), iW])

Velocity is change in position, so by definition ((d(Subscript[(OverVector[r]), i]))/(dt))=(Subscript[(OverVector[v]), i])(t).

The repulsive interaction force between i and j is

(Subscript[A, i])exp((((Subscript[r, ij])-(Subscript[d, ij]))/(Subscript[B, i])))(Subscript[(OverVector[n]), ij])

Some other simple definitional constraints include

(Subscript[d, ij])=‖(Subscript[(OverVector[r]), i])-(Subscript[(OverVector[r]), j])‖
(Subscript[(OverVector[n]), ij])=(((Subscript[(OverVector[r]), i])-(Subscript[(OverVector[r]), j]))/(Subscript[d, ij]))
(Subscript[r, ij])=(Subscript[r, i])+(Subscript[r, j])

Pedestrians touch each other whenever their distance (Subscript[d, ij]) is less than their combined radii (Subscript[r, ij]). In this case, two additional forces take effect. Body force counteracts body compression, and sliding friction impedes relative tangental motion. They are given, respectively, by

k((Subscript[r, ij])-(Subscript[d, ij]))(Subscript[(OverVector[n]), ij])
κ((Subscript[r, ij])-(Subscript[d, ij]))((Subscript[Δv, ji])^t)(Subscript[(OverVector[t]), ij])

The tangential direction and tangnetial velocity differences are given, respectively, by

(Subscript[(OverVector[t]), ij])=(-((Subscript[n, ij])^2),((Subscript[n, ij])^1))
((Subscript[Δv, ij])^t)=((Subscript[(OverVector[v]), j])-(Subscript[(OverVector[v]), i]))·(Subscript[(OverVector[t]), ij])

Combining all of the forces, we can define the interaction force between pedestrians i and j as

(Subscript[(OverVector[f]), ij])=((Subscript[A, i])exp((((Subscript[r, ij])-(Subscript[d, ij]))/(Subscript[B, i])))+k g((Subscript[r, ij])-(Subscript[d, ij])))(Subscript[(OverVector[n]), ij])+κ g((Subscript[r, ij])-(Subscript[d, ij]))((Subscript[Δv, ji])^t)(Subscript[(OverVector[t]), ij])

and between pedestrian i and a wall it is

(Subscript[(OverVector[f]), iW])=((Subscript[A, i])exp((((Subscript[r, i])-(Subscript[d, iW]))/(Subscript[B, i])))+k g((Subscript[r, i])-(Subscript[d, iW])))(Subscript[(OverVector[n]), iW])-κ g((Subscript[r, i])-(Subscript[d, iW]))((Subscript[(OverVector[v]), i])·(Subscript[(OverVector[t]), iW]))(Subscript[(OverVector[t]), iW])

Code

Initialization

(* constants *)
Dynamic[v0];
Dynamic[Np];
Dynamic[bb];
Dynamic[par1];
Dynamic[par2];
Dynamic[W];
Dynamic[Wg];
Dynamic[lo];
Dynamic[Ap];
Dynamic[Aw];
Dynamic[door];
Dynamic[κ];
m=80.0;(* pedestrian mass *)
v0=1.0; (* desired pedestrian velocity *)
τ=0.5;(* time to acquire desired velocity *)
Ap=2.0*(10^3);
Aw=1.0*(10^3);(* repulsive interaction force constant (people, wall) *)
B=0.08;(* repulsive interaction force constant *)
k=1.2*(10^5);(* body force constant *)
κ=2.4*(10^5);(* sliding friction force constant *)
Np=20;(* number of pedestrians *)
radmin=(0.5/2);(* minimum pedestrian radius *)
radmax=(0.7/2);(* maximum pedestrian radius *)
W={{{10.0,20.0},{10.0,5.75}},{{10.0,5.75},{11.0,5.75}},{{11.0,5.75},{11.0,20.0}},{{10.0,-10.0},{10.0,4.25}},{{10.0,4.25},{11.0,4.25}},{{11.0,4.25},{11.0,-10.0}}};(* list of wall segments, expressed as the coordinates of its endpoints *)
Wg={Polygon[{W[[1,1]],W[[2,1]],W[[3,1]],W[[3,2]]}],Polygon[{W[[4,1]],W[[5,1]],W[[6,1]],W[[6,2]]}]};(* graphics objects for the walls *)
door={{10.0,5.25},{10.0,4.75}};(* line segment that people will try to move towards, and directly right after that *)
dt=0.1;(* time increment for finite difference solution *)
cutoff=45.0;(* simulation time cutoff (s) *)
bb={{1.0,9.0},{1.0,9.0}};(* bounds of initial position box *)
par1=1.5;(* parameters of layout *)
par2=2.0;
colrad=0.6;(* radius of column in the column layout *)
Nw=3;(* maximum number of walls in any layout *)
lo=0;(* layout index *)
(* variables *)
Dynamic[rad];
Dynamic[radsum];
Dynamic[r];
Dynamic[v];
Dynamic[dvi];
Dynamic[dij];
Dynamic[nij];
Dynamic[tij];
Dynamic[dvijt];
Dynamic[e0];
Dynamic[fij];
Dynamic[pos];
Dynamic[pressure];
Dynamic[niw];
Dynamic[tiw];
Dynamic[diw];
Dynamic[fiw];
Dynamic[evac];
Dynamic[tpressure];
(* subroutines *)
g[x_]:=Max[0,x];(* function to cut off negative arguments *)
dijcalc[i_,j_]:=EuclideanDistance[r[[i]],r[[j]]];
dniwcalc[i_,w_]:=wallnorm[r[[i]],w];(* calculates distances and normal vectors for walls *)
fijcalc[i_,j_]:=If[i==j,{0,0},(Ap Exp[((radsum[[i,j]]-dij[[i,j]])/B)]+k g[radsum[[i,j]]-dij[[i,j]]])nij[[i,j]]+κ g[radsum[[i,j]]-dij[[i,j]]]dvijt[[j,i]]tij[[i,j]]];
fiwcalc[i_,w_]:=(Aw Exp[((rad[[i]]-diw[[i,w]])/B)]+k g[rad[[i]]-diw[[i,w]]])niw[[i,w]]-κ g[rad[[i]]-diw[[i,w]]](v[[i]].tiw[[i,w]])tiw[[i,w]];
nijcalc[i_,j_]:=If[i==j,{0,0},((r[[i]]-r[[j]])/(dij[[i,j]]))];
tijcalc[i_,j_]:={-nij[[i,j,2]],nij[[i,j,1]]};
tiwcalc[i_,w_]:={-niw[[i,w,2]],niw[[i,w,1]]};
dvijtcalc[i_,j_]:=(v[[i]]-v[[j]]).tij[[j,i]];
e0calc[i_]:=If[r[[i,1]]≥door[[1,1]],
(* past the door just go right *)
{1,0},
(* otherwise find the closest point on the door *)
Piecewise[{{door[[1]]-r[[i]],r[[i,2]]≥door[[1,2]]},{door[[2]]-r[[i]],r[[i,2]]≤door[[2,2]]}},{1,0}]
];
dvicalc[i_]:=(((v0 e0[[i]]-v[[i]])/τ)+(1/m)(Sum[If[j==i,0,fij[[i,j]]],{j,1,Np}]+Sum[fiw[[i,w]],{w,1,Nw}]))dt;(* finite difference calculation of (OverVector[(Subscript[Δv, i])]) *)
vcalc[i_]:=v[[i]]+dvi[[i]];
rcalc[i_]:=r[[i]]+v[[i]]dt;
vcap[i_]:=If[Norm[v[[i]]]≤v0,v[[i]],((v0 v[[i]])/(Norm[v[[i]]]))];(* caps speed to prevent people from moving unrealistically quickly *)
(* update process *)
update[dummy_:0]:=Module[{},
(* calculate desired directions *)
e0=Table[e0calc[i],{i,1,Np}];
(* calculate pairwise distances *)
Do[
Do[
dij[[i,j]]=dij[[j,i]]=dijcalc[i,j],
{j,1,i-1}],
{i,1,Np}];
(* calculate symmetric and antisymmetric pairwise values *)
Do[
Do[
dij[[i,j]]=dij[[j,i]]=dijcalc[i,j];(* distance *)
nij[[i,j]]=nijcalc[i,j];(* pointing vector *)
nij[[j,i]]=-nij[[i,j]];
tij[[i,j]]=tijcalc[i,j];(* tangential vector *)
tij[[j,i]]=-tij[[i,j]];
dvijt[[i,j]]=dvijt[[j,i]]=dvijtcalc[i,j];(* relative tangential velocity *)
(* forces *)
fij[[i,j]]=fijcalc[i,j];
fij[[j,i]]=-fij[[i,j]],
{j,1,i-1}],
{i,1,Np}];
(* calculate wall values *)
Do[
Do[
(*{niw[[i,w]],diw[[i,w]]}=dniwcalc[i,w];*)
niw[[i,w]]=dniwcalc[i,w][[1]];
diw[[i,w]]=dniwcalc[i,w][[2]];
tiw[[i,w]]=tiwcalc[i,w];
fiw[[i,w]]=fiwcalc[i,w],
{w,1,Nw}],
{i,1,Np}];
(* calculate changes in velocity and position *)
Do[
dvi[[i]]=dvicalc[i];(* change in velocity *)
v[[i]]=vcalc[i];(* velocity update *)
v[[i]]=vcap[i];(* cap speed (just in case) *)
r[[i]]=rcalc[i],(* position update *)
{i,1,Np}];
]
(* generate layout {wall lines, wall graphics, bounding box, door} *)
wallgen[layout_,p1_,p2_]:=Switch[layout,
0,(* single door *)
{(* W *){{{10.0,20.0},{10.0,5.0+(p1/2)}},{{10.0,5.0+(p1/2)},{11.0,5.0+(p1/2)}},{{11.0,5.0+(p1/2)},{11.0,20.0}},{{10.0,-10.0},{10.0,5.0-(p1/2)}},{{10.0,5.0-(p1/2)},{11.0,5.0-(p1/2)}},{{11.0,5.0-(p1/2)},{11.0,-10.0}}},
(* Wg *){Polygon[{{10.0,20.0},{10.0,5.0+(p1/2)},{11.0,5.0+(p1/2)},{11.0,20.0}}],Polygon[{{10.0,-10.0},{10.0,5.0-(p1/2)},{11.0,5.0-(p1/2)},{11.0,-10.0}}]},
(* bb *){{1.0,9.0},{1.0,9.0}},
(* door *){{10.0,5.0+(p1/2)-0.5},{10.0,5.0-(p1/2)+0.5}}},
1,(* column and door *)
{(* W *){{{10.0,20.0},{10.0,5.0+(p1/2)}},{{10.0,5.0+(p1/2)},{11.0,5.0+(p1/2)}},{{11.0,5.0+(p1/2)},{11.0,20.0}},{{10.0,-10.0},{10.0,5.0-(p1/2)}},{{10.0,5.0-(p1/2)},{11.0,5.0-(p1/2)}},{{11.0,5.0-(p1/2)},{11.0,-10.0}},{10.0-p2,5.5}},
(* Wg *){Polygon[{{10.0,20.0},{10.0,5.0+(p1/2)},{11.0,5.0+(p1/2)},{11.0,20.0}}],Polygon[{{10.0,-10.0},{10.0,5.0-(p1/2)},{11.0,5.0-(p1/2)},{11.0,-10.0}}],Disk[{10.0-p2,5.5},colrad]},
(* bb *){{1.0,9.0},{1.0,9.0}},
(* door *){{10.0,5.0+(p1/2)-0.5},{10.0,5.0-(p1/2)+0.5}}},
2,(* hallway with bulge *)
{(* W *){{{0.0,6.5},{7.0,6.5}},{{7.0,6.5},{10.0,6.5+3.0Tan[p1]}},{{10.0,6.5+3.0Tan[p1]},{13.0,6.5}},{{13.0,6.5},{20.0,6.5}},{{0.0,3.5},{7.0,3.5}},{{7.0,3.5},{10.0,3.5-3.0Tan[p1]}},{{10.0,3.5-3.0Tan[p1]},{13.0,3.5}},{{13.0,3.5},{20.0,3.5}}},
(* Wg *){Polygon[{{0.0,6.5},{7.0,6.5},{10.0,6.5+3.0Tan[p1]},{13.0,6.5},{20.0,6.5},{20.0,7.0},{13.0+0.5Tan[0.5p1],7.0},{10.0,6.5+0.5Sec[p1]+3.0Tan[p1]},{7.0-0.5Tan[0.5p1],7.0},{0.0,7.0}}],Polygon[{{0.0,3.5},{7.0,3.5},{10.0,3.5-3.0Tan[p1]},{13.0,3.5},{20.0,3.5},{20.0,3.0},{13.0+0.5Tan[0.5p1],3.0},{10.0,3.5-0.5Sec[p1]-3.0Tan[p1]},{7.0-0.5Tan[0.5p1],3.0},{0.0,3.0}}]},
(* bb *){{0.0,8.0},{4.0,6.0}},
(* door *){{15.0,6.0},{15.0,4.0}}}
]
(* generates a color that varies linearly from green (at input 0) to red (at input equal to max pressure) *)
colmap[x_,max_:1800]:=RGBColor[Min[(x/max),1],0.8(1-Min[(x/max),1]),0]
(* generates a normalized vector and distance for a given layout pointing from a specified wall to a specified coordinate *)
wallnorm[coord_,wall_]:=Module[{pt,far=1000},
Switch[lo,
0,(* single door *)
Switch[wall,
1,(* top wall *)
If[coord[[1]]<W[[1,1,1]],
(* left of wall *)
If[coord[[2]]≥W[[1,2,2]],
(* next to long side *)
pt={W[[1,1,1]],coord[[2]]},
(* near corner *)
pt=W[[1,2]];
],
(* right of wall *)
If[coord[[1]]>W[[2,2,1]],
(* too far to right *)
pt={far,5.0},
(* inside doorway *)
pt={coord[[1]],W[[1,2,2]]};
];
],
2,(* bottom wall *)
If[coord[[1]]<W[[1,1,1]],
(* left of wall *)
If[coord[[2]]≤W[[4,2,2]],
(* next to long side *)
pt={W[[4,1,1]],coord[[2]]},
(* near corner *)
pt=W[[4,2]];
],
(* right of wall *)
If[coord[[1]]>W[[5,2,1]],
(* too far to right *)
pt={far,5.0},
(* inside doorway *)
pt={coord[[1]],W[[4,2,2]]};
];
],
3,(* nothing *)
pt={far,5.0};
],
1,(* column and door *)
Switch[wall,
1,(* top wall *)
If[coord[[1]]<W[[1,1,1]],
(* left of wall *)
If[coord[[2]]≥W[[1,2,2]],
(* next to long side *)
pt={W[[1,1,1]],coord[[2]]},
(* near corner *)
pt=W[[1,2]];
],
(* right of wall *)
If[coord[[1]]>W[[2,2,1]],
(* too far to right *)
pt={far,5.0},
(* inside doorway *)
pt={coord[[1]],W[[1,2,2]]};
];
],
2,(* bottom wall *)
If[coord[[1]]<W[[1,1,1]],
(* left of wall *)
If[coord[[2]]≤W[[4,2,2]],
(* next to long side *)
pt={W[[4,1,1]],coord[[2]]},
(* near corner *)
pt=W[[4,2]];
],
(* right of wall *)
If[coord[[1]]>W[[5,2,1]],
(* too far to right *)
pt={far,5.0},
(* inside doorway *)
pt={coord[[1]],W[[4,2,2]]};
];
],
3,(* column *)
If[coord[[1]]>W[[5,2,1]],
(* too far to right *)
pt={far,5.0},
(* closest point on column *)
pt=W[[7]]+colrad(((coord-W[[7]])/(Norm[coord-W[[7]]])));
];
],
2,(* hallway with bulge *)
Switch[wall,
1,(* top wall *)
If[par1==0,
(* straight hallway *)
pt={coord[[1]],6.5},
(* angled hallway *)
If[coord[[1]]≤7.0||coord[[1]]≥13.0,
(* straight section *)
pt={coord[[1]],6.5},
(* middle *)
If[coord[[1]]<10.0,
(* left half *)
If[coord[[2]]>-Cot[par1]coord[[1]]+6.5+7Cot[par1],
(* left angled wall *)
pt={((7+(coord[[2]]-6.5)Cot[par1]+coord[[1]]((Cot[par1])^2))/(1+((Cot[par1])^2))),((6.5+(coord[[1]]-7)Tan[par1]+coord[[2]]((Tan[par1])^2))/(1+((Tan[par1])^2)))},
(* corner *)
pt={7.0,6.5};
],
(* right half *)
If[coord[[2]]>Cot[par1]coord[[1]]+6.5-13Cot[par1],
(* right angled wall *)
pt={((13 +(6.5-coord[[2]]) Cot[par1]+coord[[1]] ((Cot[par1])^2))/(1+ ((Cot[par1])^2))),((6.5+(13-coord[[1]]) Tan[par1]+coord[[2]] ((Tan[par1])^2))/(1+((Tan[par1])^2)))},
(* corner *)
pt={13.0,6.5};
];
];
];
],
2,(* bottom wall *)
If[par1==0,
(* straight hallway *)
pt={coord[[1]],3.5},
(* angled hallway *)
If[coord[[1]]≤7.0||coord[[1]]≥13.0,
(* straight section *)
pt={coord[[1]],3.5},
(* middle *)
If[coord[[1]]<10.0,
(* left half *)
If[coord[[2]]<Cot[par1]coord[[1]]+3.5-7Cot[par1],
(* left angled wall *)
pt={((7+(3.5-coord[[2]])Cot[par1]+coord[[1]]((Cot[par1])^2))/(1+((Cot[par1])^2))),((3.5+(7-coord[[1]])Tan[par1]+coord[[2]]((Tan[par1])^2))/(1+((Tan[par1])^2)))},
(* corner *)
pt={7.0,3.5};
],
(* right half *)
If[coord[[2]]<-Cot[par1]coord[[1]]+3.5+13Cot[par1],
(* right angled wall *)
pt={((13 +(coord[[2]]-3.5) Cot[par1]+coord[[1]] ((Cot[par1])^2))/(1+ ((Cot[par1])^2))),((3.5+(coord[[1]]-13) Tan[par1]+coord[[2]] ((Tan[par1])^2))/(1+((Tan[par1])^2)))},
(* corner *)
pt={13.0,3.5};
];
];
];
],
3,(* nothing *)
pt={far,5.0};
];
];
If[lo==1&&wall==3,
{((coord-pt)/(Norm[coord-pt])),If[EuclideanDistance[coord,W[[7]]]≤colrad,0,Norm[coord-pt]]},
{((coord-pt)/(Norm[coord-pt])),Norm[coord-pt]}
]
]
(*Manipulate[Show[Graphics[Wg],Graphics[{Text[wallnorm[pt,w],{15,9}]}],Graphics[{Red,Arrow[{pt,pt+2wallnorm[pt,w][[1]]}]}],PlotRange->{{0,20},{0,10}}],{{pt,{4.0,6.0}},Locator},{{w,1},{1,2,3}}]*)
(* restart *)
restart[layout_:0,np_:Np,p1_:par1,p2_:par2]:=Module[{},
{W,Wg,bb,door}=wallgen[layout,p1,p2];
Np=np;
par1=p1;
par2=p2;
rad=RandomReal[{radmin,radmax},np];
radsum=ConstantArray[0,{np,np}];
Do[
Do[
radsum[[i,j]]=radsum[[j,i]]=rad[[i]]+rad[[j]],
{j,1,i}],
{i,1,np}];
r=Transpose[{RandomReal[bb[[1]],Np],RandomReal[bb[[2]],Np]}];
If[layout==1,
(* check for collisions with column *)
Do[
While[EuclideanDistance[r[[i]],W[[-1]]]<colrad,
r[[i]]={RandomReal[bb[[1]]],RandomReal[bb[[2]]]};
],
{i,1,np}];
];
v=dvi=ConstantArray[0,{np,2}];
dij=ConstantArray[0,{np,np}];
nij=tij=dvijt=e0=fij=ConstantArray[0,{np,np,2}];
niw=tiw=fiw=ConstantArray[0,{Np,Nw,2}];
diw=ConstantArray[0,{Np,Nw}];
pressure=ConstantArray[0,{Floor[(cutoff/dt)]+1,Np}];
pos=ConstantArray[0,{Floor[(cutoff/dt)]+1,Np,2}];
pos[[1]]=r;
evac=tpressure=ConstantArray[0,Floor[(cutoff/dt)]+1];
generate[];
]
(* generate all positions for entire time horizon *)
generate[dummy_:0]:=Module[{step=1},
Do[
step++;
update[];
pos[[step]]=r;
pressure[[step-1]]=Table[Total[Table[Norm[fij[[i,j]]],{j,1,Np}]],{i,1,Np}];
evac[[step]]=Length[Select[r[[;;,1]],#≥11.0&]];
tpressure[[step-1]]=If[step<10,0,Min[Total[pressure[[step-1]]],1800Np]],
{t,dt,cutoff,dt}];
pressure[[-1]]=pressure[[-2]];
tpressure[[-1]]=tpressure[[-2]];
]
(* initialize *)
{W,Wg,bb,door}=wallgen[lo,par1,par2];
rad=RandomReal[{radmin,radmax},Np];
radsum=ConstantArray[0,{Np,Np}];
Do[
Do[
radsum[[i,j]]=radsum[[j,i]]=rad[[i]]+rad[[j]],
{j,1,i}],
{i,1,Np}];
r=Transpose[{RandomReal[bb[[1]],Np],RandomReal[bb[[2]],Np]}];
v=dvi=ConstantArray[0,{Np,2}];
dij=ConstantArray[0,{Np,Np}];
nij=tij=dvijt=e0=fij=ConstantArray[0,{Np,Np,2}];
niw=tiw=fiw=ConstantArray[0,{Np,Nw,2}];
diw=ConstantArray[0,{Np,Nw}];
pressure=ConstantArray[0,{Floor[(cutoff/dt)]+1,Np}];
pos=ConstantArray[0,{Floor[(cutoff/dt)]+1,Np,2}];
evac=tpressure=ConstantArray[0,Floor[(cutoff/dt)]+1];
pos[[1]]=r;
restart[];

Demonstration

Manipulate[Module[{step},
step=1+Floor[(t/dt)];
Framed[ListPlot[{Table[{4.8(tee/cutoff)+14.6,2.5((evac[[(tee/dt)+1]])/Np)+0.2},{tee,0,t,dt}],Table[{4.8(tee/cutoff)+14.6,2.5((tpressure[[(tee/dt)+1]])/(1800Np))+0.2},{tee,0,t,dt}]},PlotStyle->{Directive[Opacity[0.5],RGBColor[0,0.8,0]],Directive[Opacity[0.5],RGBColor[1,0,0]]},Joined->True,Prolog->{RGBColor[0.95,0.95,0.95],Rectangle[{14.5,0.1},{19.5,2.9}]},Epilog->Join[{Text["time = "<>ToString[PaddedForm[t,{3,1}]],{19,9},{1,0}]},Wg,{Gray},If[layout≠1,If[layout≥2,Line/@wallgen[layout,a1 Degree,p2][[1]],Line/@wallgen[layout,p1,p2][[1]]],{Line[wallgen[1,p1,p2][[1,1;;6]]],Circle[wallgen[1,p1,p2][[1,7]],colrad]}],Riffle[Table[colmap[pressure[[step,i]]],{i,1,Np}],Table[Disk[pos[[step,i]],rad[[i]]],{i,1,Np}]],{}],PlotRange->{{0,20},{0,10}},Axes->False,AspectRatio->(1/2)]]
],Style["animation controls",14,Bold],Row[{"",Button["apply changes",Module[{},
lo=layout;
v0=v00;
Ap=Ap0;
(*Aw=Aw0;*)
κ=κ0;
restart[layout,np,If[layout<2,p1,a1 Degree],p2];
lo=layout;
t=0.0;]]}],{{t,0.0,"time"},0.0,cutoff,dt,Animator,AppearanceElements->{"ProgressSlider","ResetButton","PlayPauseButton","FasterSlowerButtons"},AnimationRate->1.0},Delimiter,Style["room options",14,Bold],{{layout,0,"room layout"},{0->"single door",1->"column and door",2->"hallway with bulge"},ControlType->PopupMenu},PaneSelector[{0->Column[{Control[{{p1,1.5,"door width"},1.0,4.0,0.05,Appearance->"Labeled"}]}],1->Column[{Control[{{p1,1.5,"door width"},1.0,3.0,0.05,Appearance->"Labeled"}],Control[{{p2,2.0,"column spacing"},1.5,4.0,0.05,Appearance->"Labeled"}]}],2->Column[{Control[{{a1,30,"angle"},0,45,5,Appearance->"Labeled"}]}]},Dynamic[layout]],Delimiter,Style["simulation options",14,Bold],{{np,20,"crowd size"},{20,30,40,50}},{{v00,1.0,"pedestrian speed"},0.4,1.8,0.1,Appearance->"Labeled"},{{Ap0,2000,"pedestrian repulsion"},0,5000,100,Appearance->"Labeled"},{{κ0,240000,"sliding friction"},0,500000,1000,Appearance->"Labeled"},ControlPlacement->Left]