Gradient-based parameter calibration of an anisotropic interaction model for pedestrian dynamics

We propose an extension of the anisotropic interaction model which allows for collision avoidance in pairwise interactions by a rotation of forces \cite{arXiv:1912.04234} by including the agents' body size. The influence of the body size on the self-organization of the agents in channel and crossing scenarios as well as the fundamental diagram is studied. Since the model is stated as a coupled system of ordinary differential equations, we are able to give a rigorous well-posedness analysis. Then we state a parameter calibration problem that involves data from real experiments. We prove the existence of a minimizer and derive the corresponding first-order optimality conditions. With the help of these conditions we propose a gradient descent algorithm based on mini-batches of the data set. We employ the proposed algorithm to fit the parameter of the collision avoidance and the strength parameters of the interaction forces to given real data from experiments. The results underpin the feasibility of the method.


Introduction
Mathematical modelling of pedestrian dynamics is of practical benefit in civil engineering [23,26,27] for example in the design of complex architectures, e.g., stadiums, city centers, shopping malls, or for the management of large public events like festivals, concerts, pilgrimages, or manifestations [10].Capturing both, the individual and collective behaviours in pedestrian dynamics, is rather complex [2,7].Many different approaches have been proposed in the literature: for example models based on magnetic forces proposed by S. Okazaki and S. Matsushita in which pedestrians are modeled as magnetic charges in a magnetic field [33]; the gas-kinetic model which treats pedestrians as molecules in liquefied gas [22]; cellular automata [4,5,13]; models incorporating anticipative, rational behaviour [1,11,12] and (smooth) sidestepping [15,34].Another group of pedestrian models has emerged from the pioneering work on social forces [20] and can be classified as force-based [9,34] and the overview given in [8].
Most of these models share the property of being able to reproduce collective features such as lane formation in counterflow scenarios and travelling waves in crossing flows.Moreover, they can be used to study evacuation scenarios.On the other hand they strongly differ in their description.Indeed, some models, for the example the class of forces-based models have a sound mathematical description and allow for a statement in terms of a closed system of ordinary or partial differential equations.Others have a rather algorithmic structure because they require the solution of optimization problems to estimate for example the time-to-collision in very iteration.For the latter, a rigorous mathematical study of well-posedness is difficult.
Naturally, the modelling process is followed by an optimization or calibration procedure.For pedestrian dynamics the optimization of buildings, evacuation routes and traffic safety or the minimization of the occurence of high densities are of special interest [28,37].Moreover, with the collection of data from real world experiments grew the interest in parameter fitting for the different pedestrian models [3,17,18].This work is in a similar spirit.First, we extend the anistropic model proposed in [34] by incorporating a body size for the interacting agents.This induces another dimension of volume exclusion in the model and aims at making the model more realistic.We study the influence of the body size on the formation of lanes and travelling waves as well as the fundamental diagram of the dynamics numerically and provide a rigorous study of the well-posedness of the interaction dynamics with and without body size employing standard ODE theory.The second part of the article is concerned with the rigorous derivation of an gradient-based parameter calibration algorithm.We begin with the derivation of the first-order optimality system and propose a mini-batch gradient-descent algorithm for the calibration problem.
In more detail, the article is organized as follows: the anisotropic model for pedestrian dynamics is extended by including the agents' body size in Section 2.Moreover, we study the influence of the body size on the collective behaviour and the fundamental diagram in there.Section 3 is devoted to the parameter calibration problem.We begin with the statement of the problem, investigate its well-posedness and derive the corresponding first-order optimality conditions.Finally, the iterative gradient descent algorithm based on mini-batches for the parameter calibration is proposed in Section 4, where we show results obtained with this algorithm.We conclude with a summary of the results and an outlook on future projects.

Microscopic model with body size
We include the body size in the anisotropic model proposed in [34] as follows: let us consider a second order equation of motion with N ∈ N agents.Their positions and velocities are denoted by Moreover, the agents are assumed to have a body diameter d > 0. This leads to the following interaction dynamics where a pairwise interaction force between the agents i and j.The rotation matrix M (v i , v j ) changes the direction of the interaction force.It reads In addition, the model includes a relaxation parameter τ > 0 which controls the adaption of the current velocity v i towards the given desired velocity The rotation of the force vectors induced by the matrix M models a collision avoidance behaviour of the agents.The direction of the collision avoidance is controlled by the sign of the parameter λ.For λ > 0 agents move to the right, to avoid a collison, for λ < 0 the movement is directed to the left.See [34] for further details.
For notational convenience, the solution of the system is expresses by the vectors Remark 1.We can easily include obstacles or walls in the model, by describing them as artifical agents with fixed positions and fixed velocities and adding an additional interaction term similar to the interaction of the agents in (1b).

Well-posedness
In this section we study the well-posedness of the dynamics given in (1).We make the following assumptions on the interaction force K (d, x i , x j , v i , v j ) with i, j ∈ {1, . . ., N }.Assumption 1.The interaction forces K (d, x i , x j , v i , v j ) are locally Lipschitz and globally bounded w.r.t. the positions x i , x j and the velocities v i , v j .Assumption 2. The gradients of interaction forces ∇K (d, x i , x j , v i , v j ) exist, are locally Lipschitz and globally bounded w.r.t. the positions x i , x j and the velocities v i , v j .
Remark 2. Note that the first assumption is necessary to show to well-posedness of (1) while we need the second assumption later on to obtain the well-posedness of the calibration problem.
A key step to derive the well-posedness of the system is to establish the Lipschitz property of the right-hand side.In particular, the rotation of the force vector is of interest.

Lemma 1. For the rotation of the force term with
Proof.We introduce the short hand notation In the following we omit the dependence on d.Hence we rewrite the left-hand side of the Lipschitz condition as We estimate the first term on the right-hand side of (3) as where Analogously, we derive dM dv j .Each element of ∇M (v 2 + s(v 1 − v 2 )) is bounded, for a detailed proof of the boundedness see Appendix B. Note that K is globally bounded by Assumption 1.Moreover, K is locally Lipschitz by Assumption 1 which allows us to estimate the second term on the right-hand side of (3).Altogether, this proves the existence of the Lipschitz constants L 1 and L 2 as desired.
Theorem 1 (Existence and Uniqueness).Let Assumption 1 hold.Further we assume Then system (1) admits a unique solution Proof.On the basis of Assumption 1 and Lemma 1, the result can be directly obtained with the help of Picard-Lindelöf theorem.

Body size
In the previous discussion the body size is an abstract parameter.To give more details, we consider a variation of the Morse potential [14] leading to the interaction potential leading to the forces K given by Remark 3. Note that the forces depending on the body size as given in (5) satisfy the assumptions of Theorem 1.Hence, we have existence and uniqueness of solutions for the case with body size d > 0 as well.

Numerical studies
For the numerical studies of the model we draw random initial positions with uniform distribution in the domain and set initial velocities with respect of the desired direction of motion.We set the initial velocity to the desired velocity.Then we solve (1) with a variant of the leap frog scheme [34].Indeed, the relaxation terms are solved implicitly and the interaction is solved explicitly as given by where, ∆t denotes the step size of the time discretization.

Influence of the body size
In the following we provide some numerical results showing the influence of the body size on lane formation in the corridor and crossing scenario, respectively.The first experiment simulates the movement of two oncoming streams of pedestrians along a spacious corridor.The group of blue agents moves from left to right with desired velocity w blue = (0.7, 0) T , wereas the red group of agents moves from right to left with desired velocity w red = (−0.7,0) T .We consider N blue blue and N red red agents.Hence, the total number of pedestrians in the corridor is N = N blue + N red .The initial positions of the pedestrians x(0) = x 0 and their initial velocities v(0) = v 0 are illustrated in Figure 1a.
To assure that the pedestrians do not leave the scenario, we add reflective and periodic boundary conditions.In the corridor case the black lines (top and bottom) in Figure 1a show reflective boundaries.We model the avoidance of wall contact, by reflecting the velocity vector of an agent that would step outside of the domain in the next time step.The behavior of reflection from the wall is the same as in [34].The light blue lines illustrate periodic boundaries.Blue agents leaving the domain at the boundary on the right, enter again from the left.Analogously for the red agents.In the second scenario we consider two groups of pedestrians at a crossing.Here, the blue group of pedestrians moves from left to right with desired velocity w blue = (0.7, 0) T , and red group of pedestrians moves from bottom to top with desired velocity w red = (0, 0.7) T .In total, there are N = N blue + N red pedestrians.The initial positions and initial velocity vectors are presented in Figure 1b.There, light blue and green lines represent periodic boundaries for blue and red agents respectively.The bottom and top green lines are inflow and outflow for red agents, left and right blue lines are inflow and outflow for blue agents.Altogether, this results in Algorithm 1 for the state system (1).

Algorithm 1: state system
Input: initial positions x 0 and velocities v 0 of pedestrians, model parameters while t < T do update the positions and velocities of the agents using ( 6) check reflective boundary condition: if an agent is about to leave through a reflective border, its velocity is changed in the direction pointing inwards check periodic boundary condition: if an agent leaves through a periodic boundary, its position is adjusted due to the periodic boundary conditions t = t + ∆t end Result: trajectories and velocities at the time steps

Study for different body sizes
To analyse the simulation results for different body sizes, we fix values for the force parameters, and desired velocities of each pedestrian.The parameters are chosen to satisfy the stability ranges of the interaction force discussed in [14].In fact, in the range R/A > 1 and r/a < 1 the interaction force K is repulsive in a short-range and attractive in a long-range.This allows the distance between pedestrians to be maintained.Even if we include body size into the interaction force, it remains attractive in a short-range and repulsive in a long-range.In all simulations we see the formation of so-called traffic lanes.This formation seems to be independent of the choice of the random initial positions and velocities.It is interesting to note that even though every pedestrian is guided by simple rules for movement and interaction, phenomena arise that go beyond the behaviour of single pedestrians.Such phenomena of the self-organization are manifested in many multiagent systems [19].They were reported in many articles concerning the movement of pedestrian flows [31,36], which speaks in favour of the proposed model.Moreover, we   As the body size, ratio of repulsion and attraction amplitudes, and size of the corridor have similar effects in terms of volume exclusion, we suspect from these studies that the volume exclusion is the main driver of the lane formation process.

Fundamental diagram
Often fundamental diagrams are employed to analyse crowd motion models [27,36].Main objective is the relationship of speed and density [29,32,36] which we study for the bi- directional flow simulated with the model described in Section 2 on the corridor with 17m length and 4m width, see Figure 5 for an illustration.The density approximation is realized with the help of Voronoi diagrams as proposed in [6,32].Initially, agents move with their desired walking speed until they slow down (or speed up) due to interaction forces.Most interactions take place in the centre of the domain, it is therefore the focus of our interest.To approximate the density we construct Voronoi cells based on the positions of agents.The Voronoi diagram allows separating the area into computational grids (polygons) based on the triangulation of the computational area, in particular, the Delaunay triangulation [30].Every point on the Voronoi diagram has a region that is closer to it than any other [32].In Figure 6a we see Voronoi cells of all agents.Figure 6b    We employ the method proposed in [6,32] to compute the density as follows where A i (t) gives the Voronoi cell area for agent i, and ρ xy is the density distribution of the space.Voronoi cells are computed using the Python, Scipy module, with the Voronoi and ConvexHull methods.Note that the velocities are given by the integration of (1).
In Figure 7 we illustrate the relationship of density and speed at different time points.At the beginning, the agents move with the desired speed in regions with lower density leading to weak interaction forces.When the two groups meet, we see that as the density increases and at the same time the speed decreases as expected.The body size seems to have minor influence on this effect.This ends our analysis and numerical study of the model.In the following sections we are concerned with its parameter calibration based on real data [25].We begin with the statement of the calibration problem and analyse its well-posedness.

Analysis of the parameter calibration problem
In this section, we state the parameter calibation problem and analyse the existence of minimizers.Then we lay the theoretical ground for the formulation of a gradient-descent method by deriving the first-order optimality conditions, the existence of adjoint states and finally the identification of the gradient of the reduced cost functional.This prepares the formulation of a gradient descent-based calibration algorithm that will be employed in the numerical section to fit the parameters to real data from the BaSiGo experiment carried out in Düsseldorf in 2013 [6,25] in the next section.
Let us begin with the statement of the parameter calibration problem.Table 1 shows all model parameters.For the calibration we focus on the scaling factor λ involved in the collision avoidance process and the force strengths A and R. For fixed 1 > 0, we define the set of admissible parameters and we want to find u := (λ, A, R) ∈ U ad such that the model trajectories fit the real trajectories best.We thus consider the cost functional with x data given trajectory data from experiments.The first term measures the distance of the trajectories resulting from the model to the real trajectories from the data.The second term of the cost functional penalizes the distance of the parameters to some given reference parameters u ref .
In case there are no reference values available, we set σ 2 = 0. Note that the calibration problem is constrained by the ODE system without boundary conditions.The boundary conditions are only incorporated in the numerical simulations to reflext the domain of the experiments appropriately.
In the following we consider the spaces Y and U given by Note that both, Y and U are Hilbert spaces and U ad ⊂ U is closed.For notational convenience we define the state vector y = (x, v) ∈ Y and the state operator where F (y, u) is the vector containing the right-hand side of (1a) and (1b), respectively.

Well-posedness of the parameter calibration problem
The proof of the existence of an optimal parameter set for the calibration problem will be based on the following to Lemmata which are concerned with the boundedness of the states with respect to the control parameters and the weak continuity of the state operator e.
Lemma 2 (Boundedness).Let w i ∈ C([0, T ], R 2 ) for all i = 1, . . ., N and Assumption 1 hold and suppose the interaction forces are given by (5).For given u ∈ U ad there exists C > 0 depending only on the body size d and such that the solution y ∈ Y with e(y, u) = 0 satisfies Proof.We begin with the estimate for y L 2 (0,T,R 2 ) .It holds Hence, using Gronwall Lemma we obtain |y(t)| 2 ≤ |u| 2 e Ct and integration over [0, T ] yields y L 2 (0,T,R 2 ) ≤ C 1 |u|.For the time derivatives, we find Integration over [0, T ] leads to d dt y L 2 (0,T,R 2 ) ≤ C(1 + |u|).The two estimates together give the result.Proof.We need to show that e(y k , u k ) e(ŷ, û) as k → ∞, which can be reformulated as follows.
For any given test function ϕ ∈ C 1 c (Z), we need to obtain the following convergence property lim k→∞ e(y k , u k ), ϕ dt → e(ŷ, û), ϕ dt.
We estimate Clearly, the first integral tends to zero for k → ∞ by the weak convergence of y ∈ Y.The first term in the second integral tends to zero by continuity of the map λ → M λ (v) and the second term of the second integral tends to zero by the continuity of the interaction force with respect to A and R. Altogether, this yields the desired result.
Theorem 2. There exists at least one solution (y * , u * ) ∈ Y × U ad to (P).
Proof.The cost functional J is bounded from below and the state system is well-posed, so there exists m = inf (y,u)∈Y ×U ad J(y, u).
Let (u k ) ∈ U ad be a minimizing sequence.The sequence (u k ) ⊂ U ad is bounded and by reflexivity of U it has a weakly convergent subsequence (not relabeled) with limit û.By Lemma 2 we obtain the boundedness of (y k ) and, again by reflexivity, the existence of ŷ such that u k û in U ad and The weak continuity of e(y, u) shown in Lemma 3 implies Hence ŷ is the solution to the state equation with parameters û.By the weak lower semicontinuity of the norm, we obtain which allows us to conclude that (ŷ, û) is a minimizer of (P).
Remark 4. Because of the nonlinearity of the state equation we cannot expect uniqueness of the optimal control.
Having the existence of an optimal solution, we proceed with the derivation of the first-order necessary conditions, which will form the basis of the identification algorithm.

First-order necessary conditions
We introduce the dual pairing e(y, u), (ξ, η) where ξ, η are the Lagrange multipliers in the Banach space that represent the space of the adjoint states.Here, ξ = (ξ x , ξ v ), η = (η x , η v ) and (ξ, η) ∈ Z.The space of states Y and the set of controls U are defined at the beginning of in Section 3.
We formally derive the first-order optimality system with the help of the Lagrangian corresponding to the constrained parameter calibration problem given by

Derivation of the adjoint system and the optimality condition
Solving dL(y, u, ξ, η) = 0 yields the first-order optimality condition [21].We begin with the computation of the directional derivatives of ( 10) with respect to the state y.For notational convenience we consider x and v separately.First, we obtain Gâteaux derivatives of the objective function Now we derive the directional derivatives w.r.t. the positions of agents for the second part of Lagrange functional Using integration by parts, we obtain ) by the radial symmetry of K, we obtain Similarly, we obtain the derivative in direction [h v i ].Indeed, by using integration by parts we find To simplify the notation we introduce the operator d v i M * (v i , v j ) resulting from matrix reformulations, see Appendix A for more details.We get Moreover, the directional derivatives w.r.t. the control are given by Remark 5.For the specific choice (5), as we choose the interaction forces for the numerical results, the derivatives read

Existence of adjoint states
The proof of the existence of adjoint states is based on Corollary 1.3 in [21], which we give in Appendix C for completeness.
), i = 1, . . ., N be given, the Assumption 1 and 2 hold and u * ∈ U ad be an optimal solution of Problem (P) and let y * ∈ Y such that e(y * , u * ) = 0. Then there exist an adjoint state p * = (ξ * , η * ) ∈ Z * such that the following optimality conditions hold Proof.We check requirements given in Assumpion 3: (A2) We first note that U, Y, Z are Banach spaces.Further, J is of tracking type and therefore Fréchet differentiable [35].We are left to show Fréchet differentiability of the state operator e(y, u) : In Section 3.2.1 we computed the first variations of e(y, u) with w.r.t.y and u by . There, continuity of linear terms follows directly from the definition.We show the continuity for nonlinear terms: Let, the sequence y k := (x k , v k ) ⊂ Y have the limit y k such that x k → x and v k → v as k → ∞.Using Assumption 2 we show continuity of nonlinear terms in (11).Indeed, for the i-th component, it holds where L 1 , L 2 are Lipschitz constants.Analogously, we show the continuity for the nonlinear terms of (12) using Assumption 1 and Appendix B. These, yields to the continuity of the state operator e(y, u) w.r.t.y.The continuity w.r.t.u, can be concluded from the continuity of the λ → M, and the linearity of the force term w.r.tA and R. Altogether, this proves the continuity of e as desired.
(A3) By Theorem 1 the state system e(y, u) = 0 has an unique solution y = y(u) ∈ Y for all u ∈ V ⊂ U a neighbourhood of U ad .
(A4) We have to show that e y (y(u), u) ∈ L(Y, Z) has a bounded inverse for all u ∈ V ⊃ U ad .We can write d y e(y, u)[h] in general form where d y F (y, u) is integrable in t over every finite interval I ⊂ [0, T ] thanks to Assumption 2.
We consider d y e(y, u)[h] = r for arbitrary r ∈ Z * .By Carathéodory's existence theorem, we get for every r ∈ Z * a unique solution [16], namely h = d y e(y, u) −1 [r].
Using d y e(y(u), u), we obtain and with the help of Gronwall's Lemma we get the boundedness of the inverse of d y e(y(u), u).
Altogether, the requirements of Assumption 3 are satisfied, and thus the Proposition 1 (see Appendix C) yields the existence of an adjoint state.
Remark 6.Note that assuming more regularity of the adjoint states, for instance Z = Y , we may establish the strong formulation of the adjoint system given by supplemented with the terminal conditions ξ 1,i (T ) = 0, ξ 2,i (T ) = 0.The strong form will be employed for the numerical results in Section 4.

Gradient of the reduced cost functional
To minimize the objective function we aim to apply gradient descent algorithm.In order to determine the gradient we define the control-to-state operator F : U → Y , and introduce the reduced cost functional as Ĵ(u) := J(F(u), u).Using e(y, u) = 0, we obtain 0 = d y e(F(u), u)[dF(u)] + d u e(F(u), u)).
Taking the derivative of the Lagrangian with respect to the state y, we get d y e(y, u) * ξ = −d y J(y, u).
With these, we can compute the Gâteaux derivative of the reduced cost functional in the direction h u ∈ U , and obtain Note that we have already computed d u J(y, u) in Section 3.2.1.This allows us to identify the gradient of the reduced cost functional as

Calibration algorithm and results
To calibrate the control parameters we use real data from the BaSiGo experiment carried out in Düsseldorf in 2013 [6,25].In the following we denote the trajectories from experimental data by For the corridor case, the data we take from file "bi_corr_400_a_02.txt" in [25].These show the positions of the pedestrians in the domain Ω = [−6, 6] × [0, 4.2], over time t ∈ [0, 150] seconds.For the crossing case, the data we take from file "CROSSING_90_E_2.txt" in [25].This file provides the positions of the pedestrians in the crossing corridors ]), over time t ∈ [0, 283] seconds.However, for the calibration procedure, we use only 8-second intervals for both scenarios.

Numerical schemes and steepest descent algorithm
In general, the nonlinearities make it rather difficult to solve the optimality system all at one.We therefore opt for an iterative approach to compute the gradient of the calibration problem.Indeed, we first solve the state system (1) as considered in Section 2.1.1.Then, we integrate the adjoint system (15) with the help of a second-order Runge-Kutta method backward in time.Here, we use the same time steps as for the state problem and transform the time via s = T − t to recover an initial value problem.With the state solution and the adjoint solution, we calculate the gradient using (16), where the integral is approximated with the trapezoidal rule.
We apply a steepest descent algorithm to update the control parameters in every iteration The initial positions of the agents x 0 coincide with the initial positions of the experimental data, which are distributed in the domain [−6, 6] × [0, 4.2].As the initial velocity of the pedestrians, we set the average velocity from the experimental data.We probably induce some error here, as we do not have the exact values from the data.
As initial guess of the control parameters we take u 0 = (0, 0, 40) for both experiments.We set the regularisation parameters in the cost functional to σ 1 = 1 and σ 2 = 0, we therefore do not need to choose reference values for the control.In this setting we perform 100 iterations with Algorithm 2. The step size parameter is set to β k = (20, 4000, 4000) to account for the different ranges of the control values.
Figure 8a illustrates the decrease of the cost functional for the corridor case.It starts from approximately 6.51 in the first iteration and terminates around 5.11 in the last iteration.In Figure 8b we see the evolution of the cost for the crossing scenario from 9.03 to 7.34.We note that the cost values of the corridor are smaller than the ones of the crossing case.This indicates that the proposed model captures effects of counter-flow better than effects of crossing flow.In addition, we detect that the evolution of the cost is smoother for the crossing case.This could be for the same reason.
The parameters of the calibration procedure reached to optimal values λ = −0.17,A = 5.56, and R = 29.23 for the simulation in the corridor.For simulation at the crossing optimal values reached λ = −0.18,A = 10.41, and R = 29.02.Interestingly, the optimal λ is negative and in the same range in both scenarios.This indicates that the pedestrians in the data set have a tendency to move to the left to avoid collisions.Moreover, the repulsion strengths are similar for both scenarios, but the values for the attraction strengths differ.The difference of the attraction strengths may arise from the post-interaction behaviour of the model.Simulations as reported in [34] show that two interacting agents move on with slightly shifted positions after an collision avoiding interaction.We suspect that this has more impact on the results in the crossing than the corridor scenario.So far, however, this is only a conjecture and needs to be proven or discarded by further investigations.

Conclusion and outlook
We extended the anisotropic interaction model proposed in [34] by including a body size and thus additional volume exclusion effects.Numerical studies indicate that the body size has influence on lane formation process.In fact, it seems that a smaller body size leads to a higher number of lanes formed and vice versa.Moreover, we investigated the fundamental diagram of the dynamics and found that higher densities lead to lower velocities.This observation was expected from other experiments and underlines the feasibility of the approach.Due to the ODE formulation of the model, we were able to analyse the model in terms of well-posedness and further to rigorously derive a gradient-based descent algorithm for a calibration problem using real data.The optimal scaling parameters for the collision avoidance and the repulsion strength turn out to match very well for the two tested scenarios.For the attraction parameter we find different values for the two scenarios.We suspect that this is related to the post-collision behaviour of the model.
To prove or discard this conjecture is interesting future work.Moreover, the rigorous analysis of stationary states like the lanes in the corridor case and the travelling waves in the crossing case is planned for future investigations.

A Tensor
Here, we present the details of the mathematical manipulations made to tensor d v i M (v i , v j ) in the derivation of d v i M * (v i , v j ).Note that d v i M (v i , v j ) is a tensor with dimension 2 × 2 × 2, given by where , for v i = 0, v j = 0, 0, else .
The elements of the tensor d v i M (v i , v j ) in ( 18) we denote as 12) is the transposed tensor after swaping its' axes, as given below

Figure 1 :
Figure 1: Initial positions and initial velocity vectors for the two scenarios with parameters N = 80, N blue = 40, N red = 40, d = 0.2.

6 Figure 2 :
Figure 2: Simulation results in the corridor by different body size of pedestrians at time T = 35.On each simulation we fix parameters: A = 5, R = 20, a = 2, r = 0.5, λ = 0.25.Desired velocities for red and blue agents are w red = (−0.7,0) T and w blue = (0.7, 0) T respectively.Time step in the Leap-Frog Scheme is ∆t = 0.00625.

Figure 2
Figure 2 shows simulation results of the corridor scenario for different body sizes.The results indicate a relation between the body size and the number of lanes formed.The smaller the body size, the more lanes are obtained.The parameters used for the simulation are reported in the caption of the figure.Similar results are found for the crossing scenario in Figure 3. Again, the smaller the body size, the more lanes are formed.In all simulations we see the formation of so-called traffic lanes.This formation seems to be independent of the choice of the random initial positions and velocities.It is interesting to note that even though every pedestrian is guided by simple rules for movement and interaction, phenomena arise that go beyond the behaviour of single pedestrians.Such phenomena of the self-organization are manifested in many multiagent systems[19].They were reported in many articles concerning the movement of pedestrian flows[31,36], which speaks in favour of the proposed model.Moreover, we

5 Figure 3 :
Figure 3: Simulation results at the crossing by different body size of pedestrians at time T = 35.On each simulation we fix parameters: A = 5, R = 20, a = 2, r = 0.5, λ = 0.25.Desired velocities for red and blue agents are w red = (0, 0.7) T and w blue = (0.7, 0) T respectively.Time step in the Leap-Frog Scheme is ∆t = 0.00625.

Figure 4 :
Figure 4: Simulation results in the corridor by different force parameters at time T = 35.On each simulation body diameter of the agents are fixed: d = 0.5.Time step in the Leap-Frog Scheme is ∆t = 0.00625.
shows the same cells, but only the ones in the region of interest Ω = [−2, 2] × [0, 4].Outside of the domain Ω agents are involved in less interactions and walk approximately with their desired velocities.The model parameters in these simulation are set λ = 0.07, A = 0.2, R = 2, a = 1, r = 0.1.

Table 1 :
Model parameters λ scaling factor (collision angle) τ relaxation parameter (desired velocity) A attractive force amplitude R repulsive force amplitude a attractive force range r repulsive force range w desired velocity d body size of pedestrianTo study the well-posedness we use the following notion of optimality:Definition 1.We call u ∈ U ad optimal, if it is a solution to the optimization problem