0
$\begingroup$

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.enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Richard Feb 7 at 20:52
  • $\begingroup$ 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. $\endgroup$ – Jerome15 Feb 7 at 22:13
  • $\begingroup$ 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. $\endgroup$ – Richard Feb 7 at 23:59
  • $\begingroup$ 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. $\endgroup$ – Richard Feb 8 at 0:00
  • $\begingroup$ Welcome to scicomp! A tip: You can use MathJax to typeset your mathematical formulas. This will make the question much easier to read. $\endgroup$ – Mauro Vanzetto Feb 9 at 11:18

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.