Social Force Model for Pedestrian Dynamics by Euler Method

The social force model is a model using Newtonian forces to describe the movement of individuals. As seen page 1

Each individual feels the following forces:

1. A driving force towards the goal

$$m_i\frac{v_i^0(t)\mathbf{e}_i^0(t)-\mathbf{v}_i(t)}{\tau_i}$$

where $$v_i^0$$ is the desired velocity of the individual, $$\mathbf{e}_i^0$$ is the direction of the goal, and $$\mathbf{v}_i$$ is the actual velocity of the individual. Here, $$\tau_i$$ is the relaxation time, meaning the time needed for an individual to reach the desired velocity.

1. A repulsive force coming from other individuals and obstacles $$\mathbf{f}_{ij}=\left\{ A_i \exp\left[\frac{r_{ij}-d_{ij}}{B_{ij}}\right] +kg(r_{ij}-d_{ij})\right\}\mathbf{n}_{ij}$$

where $$A_i$$ is the repulsion coefficient, $$B_{ij}$$ is the distance before the individual feels the repulsion, $$(r_{ij}-d_{ij})$$ is the distance between two individuals, $$k$$ is an elastic constant, and $$g$$ is a function depending on the distances between the two individuals. What $$kg(r_{ij}-d_{ij})$$ does, it pushes more the individual that is coming too close to the selected individual. $$\mathbf{n}_{ij}$$ is the direction of the repulsive force.

Finally, I get the following equation (1) using Newton second law:

$$m_i\frac{d\mathbf{v}_i}{dt}=m_i\frac{v_i^0(t)\mathbf{e}_i^0(t)-\mathbf{v}_i(t)}{\tau_i}+\sum\limits_{j,j\ne i}\mathbf{f}_{ij} \tag{1}$$

This equation is applied to each individual in the room, and my goal is to solve it.

To do so, I use the explicit Euler method, projecting the forces on the X and Y axis separately.

My reasoning is the following:

$$\ddot{x}[n]=\frac{\dot{X}[n+1]-\dot{X}[n]}{h}$$

where $$\dot{X}$$ is the speed projected on the X-axis, and $$\ddot{x}$$ is its derivative. Then, (2) is derived by dividing (1) by $$m$$

$$\ddot{x}=\frac{vi^0\cos{\theta}-\dot{X}\cos{\alpha}}{\tau_i}+\frac{h}{m}\sum\mathbf{f}_{ij}\cdot \hat{x} \tag{2}$$

Here, $$\mathbf{f}_{ij}\cdot \hat{x}$$ is the projection of the sum of repulsive forces the $$i$$th individual feels on the X-axis, $$\theta$$ is the angle of directions of the driving force, $$\alpha$$ is the angle of the direction of the individual, set to random at first.

By injecting (2), I get

$$\dot{X}[n+1]=\dot{X}[n]+\ddot{x}[n]h$$ $$\dot{X}[n+1]=\dot{X}[n]\left(1-\frac{h}{\tau_i}\right)\cos{\alpha}+\frac{h}{\tau_i}vi^0\cos{\theta}+\frac{h}{m}\sum\mathbf{f}_{ij}\cdot \hat{x}$$

Similarly, for $$Y-axis$$ by switching $$\cos$$ to $$\sin$$ in (2).

I get the actual positions coordinates by applying: $$X[n+1]=X[n]+\dot{X}[n]h$$

All that is left is to repeat those instructions for $$N$$ steps for $$X$$ people.

Here is my code:

from math import*
from numpy import*
from copy import*
import matplotlib.pyplot as plt
##Condition initiale
purpose=[50,50] #Coordonné du sommet du rectangle [600,300]
purpose_exit=[purpose[0]+1E3,purpose[1]+20] #La seconde sortie X pixel plus loin
rayon=5
number_of_particles = 1 #Nb d'individus généré aléatoirement
v_désiré=1
m=70 #masse individu
A=200 #coefficient de répulsion
B=0.08 #coefficient d'interraction
τ=0.5

##Taille de l'objectif
lenght_goal=rayon*2 #Longueur du rectangle (Horizontale)
width_goal=rayon*4 #Largeur du rectangle (Verticale)

##Particle object
class Particle: #Ensemble tel que tout les objets qui appartiennt à cette ensemble ont les mêmes caractéristiques
def __init__(self, x, y, v, angle, rayon, n, type): #Initialisation du blueprint de l'objet désigné "self"
self.x = x #position
self.y = y
self.rayon = rayon
self.couleur = (0, 100, 255)
self.thickness = 1
self.v = v #vitesse
self.v_d = v_désiré #Vitesse désiré
self.angle = angle #angle d'orientation de la vitesse
self.n = n #numéro individu
self.goal = True #Indicateur de comportement, True vise la sortie, False non
self.goal1 = False #Pour l'autre sortie
self.type = type
class Barrier:
def __init__(self,x1,y1,x2,y2,couleur):
self.x1 = x1
self.y1 = y1
self.x2 = x2
self.y2 = y2
self.couleur = (0,200,200)
self.thickness = 6
self.particules = []
l = sqrt((x2-x1)**2+(y2-y1)**2) #longueur du mur
dx=x2-x1
rayon_obs = rayon*2
#Ajout auto des obstacles aux barrières
if (x2-x1)!=0 and (y2-y1)!=0:
c=round((y2-y1)/(x2-x1)) #coefficient directeur arrondie due aux pixels => Entiers
d=round(y1-c*x1)
for i in range(x1,x2+1,rayon_obs):
self.particules+=[Particle(i,(c*i+d),0,0,rayon_obs,i,"obstacle")] #Pour l'affichage
elif (x2-x1)==0: #Droite verticale
for i in range(1,int(l/rayon_obs)):
self.particules+=[Particle(x1,y1+rayon_obs*i,0,0,rayon_obs,i,"obstacle")] #Pour l'affichage
elif (y2-y1)==0: #Droite horizontale
for i in range(1,int(l/rayon_obs)):
self.particules+=[Particle(x1+rayon_obs*i,y1,0,0,rayon_obs,i,"obstacle")] #Pour l'affichage
def but(self,purpose1): #Vecteur de direction
but_y=purpose1[1]-self.y #Différence de vecteurs résultant en la direction de sortie
but_x=purpose1[0]-self.x
return math.atan2(but_y,but_x)

def f_α_β(p1,p2): #force de répulsion individuelle
return A*exp((-d_α_β(p1,p2))/B),angle_α_β(p1,p2)

def d_α_β(p1,p2): #Calcul la distance entre l'individu α et β
dx = p1.x - p2.x
dy = p1.y - p2.y
return sqrt(dx**2+dy**2) #Correspond à la distance séparant 2 particules (la norme de cette dernière)

def angle_α_β(p1,p2): #Direction de la force de répulsion
dx = p1.x - p2.x
dy = p1.y - p2.y
tangent = math.atan2(dy, dx) #Angle de rotation de collision des particules
angle = tangent #Prise en compte du décalage +pi/2
return pi+angle

def Euler_s():
global Position,Vitesse,Angle_v,Ẋ,Ẏ
h=T/N #le pas
for n in range(len(t)): #Boucle de temps
Vitesse+=deepcopy([Vitesse[0]])
Angle_v+=deepcopy([Angle_v[0]])
Position+=deepcopy([Position[0]])
Ẋ+=deepcopy([Ẋ[0]])
Ẏ+=deepcopy([Ẏ[0]]) #Les valeurs seront écrasés ensuite, il s'agit d'un place-holder
for j in range(len(my_particles)): #Boucle d'individu α
if my_particles[j].type=="individu": #On n'applique pas Euler aux obstacles
for z in range(len(my_particles)): #Calcul de la composante de répulsion
Σ_X=0
Σ_Y=0
if z!=j: #L'individu j ne se repousse pas lui-même
A=f_α_β(my_particles[j],my_particles[z]) #1 appel par rapport à une particule
norme, β = A[0], A[1]
# Σ_X+=norme*cos(β)
# Σ_Y+=norme*sin(β)
# print(Σ_X,Σ_Y)
θ = but(my_particles[j],purpose)
α = Angle_v[n][j]
##Calcul nouvelle Vitesse
Ẋ[n+1][j]=Ẋ[n][j]*(1-(h/τ)*cos(α))+(h/τ)*v_désiré*cos(θ)+(h/m)*Σ_X #Composante X de la Vitesse à l'instant t+1, de la particule j
Ẏ[n+1][j]=Ẏ[n][j]*(1-(h/τ)*sin(α))+(h/τ)*v_désiré*sin(θ)+(h/m)*Σ_Y
Vitesse[n+1][j]=sqrt(Ẋ[n+1][j]**2+Ẏ[n+1][j]**2) #Norme de la vitesse
arc_tangent = math.atan2(Ẏ[n+1][j], Ẋ[n+1][j]) #Angle de vitesse
Angle_v[n+1][j]= arc_tangent #Prise en compte du décalage +pi/2#Angle de la vitesse
my_particles[j].angle=Angle_v[n+1][j]
##Calcul nouvelle Position
X_1=Ẋ[n][j]*h+Position[n][j][0]
Y_1=Ẏ[n][j]*h+Position[n][j][1]
Position[n+1][j]=[X_1,Y_1]
my_particles[j].x=X_1
my_particles[j].y=Y_1
return

##Historique des grandeurs physiques
N=10000 #Nb de points, où N>>>T
T=120 #Temps final
t=linspace(1,T,N) #Nb de pas de temps
# Position=[[(100,100),(100,120),(120,100),(120,120)]] #Chaque lignes de la matrices correspond à un instant
# Position_backup=[[[rayon+80,rayon],[rayon+20,rayon+20],[rayon+20,rayon],[rayon,rayon+20]]] #Chaque lignes de la matrices correspond à un instant
Position_backup=[[[100,0]]]
Position=[[]]
Vitesse=[[]]
Ẋ=[[]] #Composante de la vitesse
Ẏ=[[]]
Ẋ[0]=number_of_particles*[0]
Ẏ[0]=number_of_particles*[0]
Angle_v=[[]] #A déterminé en fn de la pos° des individus

##Génération des individus
my_particles = []
for n in range(number_of_particles): #Attribution des caractères aléatoires des particules
Position[0]+=[Position_backup[0][n]]
Vitesse[0]+=[0] #Vitesse initiale nulle
x = Position[0][n][0] #Lit le 1er élèments de la liste, indexé par le n° de l'indi, puis par x/y
y = Position[0][n][1]
v=Vitesse[0][n]
angle_v = math.atan2(y,x) #Direction du vecteur vitesse, tous repéré par rapport à l'origine (Voire comment fonction atan2)
Angle_v[0]+=[angle_v]
type="individu"
my_particles.append(Particle(x, y, v, angle_v, rayon, n, type))

Euler_s()
##Drawing
fig = plt.figure()
n=0
X=[]
Y=[]
for i in range(len(Position)):
X+=[Position[i][n][0]]
Y+=[Position[i][n][1]]
plt.plot(X,Y,'ro')
plt.plot(purpose[0],purpose[1],'gs')
plt.axis([0, 100, 0, 100])
fig.show()


You can directly copy it and run it for yourself. And that's where's my problem.

If I put an individual in some position, such as

Position_backup=[[(80,20)]]


, but his goal is at

purpose=[50,50]


he just won't go there.

I want that person to go towards that green goal. He should do it in a straight line.

I don't understand why it doesn't work.

Edit: I moved this topic from Stack Overflow on the advice of LutzL.

• Your code doesn't run. linscape is not defined. number_of_particles is not defined. rayon is not defined. In order to get help, you should post a minimum working example that demonstrates the issue. – Richard Feb 7 at 20:52
• That's strange, it works for me. Does the fact that I use miniconda3 makes it different for other users ? Furthermore, I think I cleary defined rayon and number_of_particles at the start of this code, whereas for linspace, I wrote from numpy import* which should give access to that function. – Jerome15 Feb 7 at 22:13
• It seems to be working now. One other question: I wonder if you could translate the comments into English? I suspect that will make them accessible to a wider audience. – Richard Feb 7 at 23:59
• I also see that you have global arrays called Position and also arrays of Particle. This duplication of data could easily lead to problems. I suggest dropping Position and keeping only particles. – Richard Feb 8 at 0:00
• Welcome to scicomp! A tip: You can use MathJax to typeset your mathematical formulas. This will make the question much easier to read. – Mauro Vanzetto Feb 9 at 11:18