A micro-macro hybrid model with application for material and pedestrian flow

Abstract: In this paper, a hybrid modeling approach for granular flow-like applications is presented. The approach allows to switch for a priori fixed points in time between the different levels of description which are the microscopic andmacroscopic scale, respectively. Based on the numerical discretization of the models, the switching procedure is able to interpret information on individual objects as density distributions and vice versa. In particular, the reverse direction, i.e. from a macroscopic to a microscopic perspective, requires the solution of a nonlinear least squares problem subject to further constraints. Simulation results are given and demonstrate the good performance of the algorithm in the case of material and pedestrian flow models.


Introduction
Multi-scale models are a common tool to simulate the behavior of granular materials in many engineering applications. Such models are for example used to simulate pedestrian flow (Cristiani, Piccoli, & Tosin, 2014;Etikyala, Göttlich, Klar, & Tiwari, 2014), granular flow (Cristiani, Piccoli, & Tosin, 2011), traffic flow on roads (Herty & Moutari, 2009;Moutari & Rascle, 2007) or material flow on conveyor belts (Göttlich, Hoher, Schindler, Schleper, & Verl, 2014). Multi-scale models combine ABOUT THE AUTHOR The Scientific Computing Research Group belongs to the Department of Mathematics at the University of Mannheim. It was founded in 2011 and directed since then by Prof. Göttlich. The research focuses on the mathematical modeling, numerical simulation and optimization of dynamic processes with applications in manufacturing systems, traffic flow and pedestrian dynamics and power grids.

PUBLIC INTEREST STATEMENT
Conveyor belts are commonly used in many industries to transport a large amount of goods. Typical applications include for example bottling, canning and packaging. We introduce a hierarchy of models enabling to consider different levels of simulation accuracy. Our intention is to develop a hybrid numerical framework which is able to switch between the scales at user-defined points in time. This allows for either a fine or coarse resolution of the underlying equations to detect collisions or identify free flow regimes. Based on the numerical approximations, the switching procedure is able to interpret information on individual goods as density distributions and vice versa. In particular, the reverse direction requires the solution of a restricted nonlinear least squares problem. Simulation results are presented to demonstrate the good performance of the algorithm in the case of material flow models. the microscopic scale, where the trajectories of single objects are considered separately (e.g. pedestrian flow (Helbing, Farkas, & Vicsek, 2000), granular flow (Cleary & Sawley, 2002;Cundall & Strack, 1979)), with the macroscopic scale, where the entity of objects is treated as a mass distribution (e.g. (Colombo, Garavello, & Lécureux-Mercier, 2012;). Typically, multi-scale models are derived from the microscopic scale via hydrodynamic limits. An overview on multi-scale models for crowd dynamics is given in (Bellomo, Piccoli, & Tosin, 2012). Hybrid multi-scale traffic flow models have been investigated in Herty & Moutari, (2009) and Moutari & Rascle, (2007). However, the transition between the different scales has been mostly analyzed from a theoretical point of view and there exist only a limited number of contributions, where simulation techniques have been applied to switch the dynamics (Herty & Moutari, 2009;Moutari & Rascle, 2007). Our intention is now to develop a numerical framework which is able to switch between the microscopic and macroscopic scale at fixed user-defined points in time. This allows for either a fine or coarse resolution of the underlying model equations depending on special situations (e.g. queuing effects).
For the presentation of the micro-macro hybrid model within this article, the formal derivation of limiting macroscopic equations is omitted and only briefly explained. The focus is on a hybrid switching method to convert the different scales without loosing information and properties of the model. In the case of material and pedestrian flow problems, we present the model equations, introduce suitable discretizations methods and explain how special model characteristics are incorporated in the implementation of the switching algorithm. Different to (Herty & Moutari, 2009;Moutari & Rascle, 2007), the switching procedure works directly in Eulerian coordinates and the problems are defined on a two-dimensional domain. Figure 1 illustrates the switching directions, i.e. from the microscopic to the macroscopic scale and vice versa. For the switch from the microscopic to the macroscopic representation, the idea is to place a Gaussian bell curve corresponding to the mass of each object. The total sum of all curves then leads to a density distribution with high density in dense regions and low density in regions, where the objects are widely spread. The transformation from the macroscopic to the microscopic scale is more involved since the precise positioning is in general not possible. In order to have the objects placed as dense as possible in regions of high density, the idea is to fill the geometry completely and to remove objects according to the density distribution. Therefore, the density is interpreted as a probability density and a nonlinear least squares problem is solved to meet the given density distribution as best as possible. An update algorithm ensures that the objects do not overlap and the placement is admissible.
The article is organized as follows. Section 1 is concerned with the general description of multiscale models. The model equations are presented and embedded into the context of hybrid switching. In Section 2, the model switch is applied to a material flow model ) and a pedestrian model (Helbing et al., 2000). Section 3 deals with the implementation of the multi-scale model and the hybrid switching algorithm. Finally, in Section 4, numerical examples are analyzed to evaluate the accuracy of the switching method.

Micro to Macro
Macro to Micro

Multi-scale models
In this section, we present the microscopic and macroscopic scale of the hybrid multi-scale model in a general framework as well as the transition between the levels. Furthermore, the switching idea is introduced based on some key characteristics.
Multi-scale models are commonly used to simulate the evolution of objects over time under certain conditions and inside a given geometry. However, assuming a high number of objects, there might be scenarios, where statistical values of the averaged behavior of objects such as flow (objects per time) or density (objects per area) are sufficient (see e.g. Göttlich et al., 2014)). The general multi-scale model considered here consists of two simulation levels: a microscopic one, described by Newton's equations of motion, and a macroscopic one relying on a hyperbolic conservation law in two space dimensions (Bellomo et al., 2012;Göttlich, Klar, & Tiwari, 2015).
The setting we are interested in are bounded domains with time-independent external forces. The latter can be exerted for example by the boundaries of the geometry, the friction of a conveyor belt (see Section 2.1) or a desired velocity for pedestrians (see Section 2.2). Other external forces (e.g. gravitational forces) are also possible. The objects are assumed to be uniform and indistinguishable, meaning that they have the same mass, shape and set of parameters. They are assumed to be incompressible and are not allowed to overlap.

From microscopic to macroscopic scale
The microscopic model is based on Newton's equations of motion (e.g. (Bellomo et al., 2012)). Being N the number of objects in the system, x i 2 R 2 ; i ¼ 1; . . . ; N the spatial coordinate, v i 2 R 2 the velocity of object i and m the mass of a single object, the equations are characterized by together with an initial condition (1:2) Equation (1.1a) describes the velocity of the objects and Equation (1.1b) includes all forces influencing their acceleration. The force term F : R 2 Â R 2 ! R 2 collects the interaction forces, whereas the force term G : R 2 Â R 2 ! R 2 describes exterior forces as for example interactions with boundaries or friction forces.
To derive the limiting macroscopic equations of (1.1), the following formal computations have to be done. We refer to Carrillo et al. (2009);Etikyala et al. (2014); Göttlich et al. (2015); Ruzhansky, Cho, Agarwal, & Area (2017) where such considerations are made for swarming, pedestrian and material flow models in a more detailed way.
We start with the Liouville equation where f ðNÞ ¼ f ðNÞ ðx 1 ; :::; x N ; v 1 ; :::; v N Þ describes the distribution function in the phase space. We integrate over dΩ 1 ¼ dx 2 Á Á Á dx N Â dv 2 Á Á Á dv N and make use of Equation (1.1b). We consider the chaos assumption f ð2Þ ðx 1 ; v 1 ; x 2 ; v 2 Þ ¼ f ð1Þ ðx 1 ; v 1 Þf ð2Þ ðx 2 ; v 2 Þ, set F 12 ¼ Fðx 1 À x 2 ; v 1 À v 2 Þ and let N go to infinity, then the resulting equation for the distribution function f ðx; v; tÞ is the mean field equation We define the density Then, the temporal derivative of the density ρ leads to the continuity equation and the computation of the temporal derivative of ρu yields the momentum equation (1:5) A monokinetic closure function f ,ρðx; tÞδ uðx;tÞ ðvðtÞÞ is used to close the system, i.e. fluctuations are neglected. Moreover, we ignore the velocity dependence meaning that the elastic normal forces dominate the friction forces. So Equation (1.5) reduces to the hydrodynamic limit equation (1:6) which describes the system together with the continuity Equation (1.4).
To reduce the system to a scalar limit equation, a quasi-stationary approach is considered, where the temporal and spatial derivatives are set to zero. We get Gðx; uÞ ¼ ÀðF Ã ρÞðxÞ: If the forces F and G are known, we get an expression for u, the so-called closure velocity (analogous to (Göttlich, Knapp, & Schillen, 2017)) that is of the form (1:7) where the first part depends explicitly on the density ρ and the second term only depends on the spatial variable x. We plug the closure velocity u Ã into the continuity Equation (1.4) to get the scalar limit equation which defines the macroscopic model.
The macroscopic model is obviously given by a hyperbolic conservation law in two space dimensions (Aggarwal, Colombo, & Goatin, 2015). The density distribution is denoted by ρðx; tÞ with spatial coordinate x 2 R 2 and time variable t 2 R þ . The velocity function corresponds to the closure velocity (1.7). Note that the model is nonlocal since u Ã ðρ; xÞ does not only depend on ρðxÞ but on the density distribution ρ over the whole geometry. The boundaries of the domain are denoted by @Ω W and the inner normal vector is called nðxÞ. The model Equation (1.8) together with the closure velocity (1.7) and additional initial and boundary conditions can be expressed as follows: hu Ã ðρ; xÞ; nðxÞi ¼ 0; x 2 @Ω W : (1:9d) Equation (1.9a) describes evolution of density governed by the hyperbolic conservation law, see Equation (1.8). The initial data is given by (1.9b), the inflow condition is stated by (1.9c) and the condition that no mass is going through the boundaries is ensured by (1.9d).

Hybrid switching method
The goal of the hybrid switching method is to provide a computational framework to change the scale from the microscopic description (1.1) to the macroscopic (1.9) and vice versa. The switch is directly based on the Eulerian representation, contrary to the situation considered in (Herty & Moutari, 2009;Moutari & Rascle, 2007), and works for a priori defined fixed points in time.
As one might imagine, the transformation from the microscopic to the macroscopic level is straightforward since all necessary information is available. Averaged quantities such as density or flux can be computed in a straightforward way and are sufficient to describe the approximate behavior of the system over time. However, the lack of missing information is the more challenging problem while transforming the system back. Therefore, the main task is to find a method to place the objects inside a given domain such that the given density distribution matches the position of objects in the microscopic scale. For application purposes, the transformation from the macroscopic to the microscopic representation is necessary whenever the trajectory of single objects is needed.
The consistent relation between the amount of mass and the number of objects is ensured by with m tot being the total mass in the system, r the radius of an object and N again the total number of objects. Note that this number is independent on the switching direction and the derivation is as follows: If we have a domain of maximal density ρ max , the corresponding object distribution is as dense as possible which is a hexagonal configuration (Degond, Ferreira, & Motsch, 2017). Then, the volume corresponding to one object corresponds the volume of a hexagon with inner circle radius r.
The area of such an hexagon is A hex ¼ 2 ffiffiffi 3 p r 2 . The volume for one object therefore is (1:11)

Switch from microscopic to macroscopic scale
The transformation from the microscopic to the macroscopic representation is performed such that the original position of objects can be interpreted as a density distribution. Thus, the idea is to spend a two-dimensional Gaussian bell curve for each object representing the density distribution and sum them all up. Since the volume of a two-dimensional Gaussian function with variance σ 0 is normalized and the volume per object has to be equal to V object , we have to multiply by the corresponding volume as done in Equation (1.11). The x-resp. y-coordinates of the objects are given by the vectors X o and Y o and the continuous density distribution for one object i is then ( 1:12) where X o;i and Y o;i is the i-th component of the vector X o and Y o , respectively. Summing the densities for each object (1:13) we end up with the desired continuous density distribution.

Switch from macroscopic to microscopic scale
The switch becomes more involved when the system is transformed from the macroscopic to the microscopic representation. This is due to missing detailed information on individual objects at the macroscopic level. To remedy this drawback, we need to solve an optimization problem to find a good guess for the microscopic density interpretation. To do so, the density distribution is considered as a normalized probability distribution, meaning that the maximal density is ρ max ¼ 1: Since it is not possible to compare the placement of objects to the macroscopic density distribution ρ mac , the object distribution has to be transformed back to a density distribution called ρ app , cf. (1.13). Given the position of objects X o ; Y o f g respecting the boundaries and the non-overlapping conditions, we can set up an optimization problem of the following form: where Ω ad is the admissible interior domain such that objects placed here do not interact with the boundaries. The optimization problem (1.14) is in fact a nonlinear least squares problem subject to the constraints that the total numbers of objects is met (1.14b), the objects are placed inside the admissible domain (1.14c) and do not overlap (1.14d). For the numerical implementation, the constraints (1.14c) and (1.14d) are ensured by an update algorithm and are not considered directly in the solution of the optimization problem.
A crucial point is to find an admissible start solution ρ app : Therefore, we require that high-density regions, where ρ max is attained, are packed as dense as possible such that the geometry is filled completely and then remove objects until the number N ¼ mtot V object is reached. First, objects, where the density is lower than a certain tolerance value , are removed. For any other placement, remove the object with a probability of ð1 À ρðxÞÞ, where x is the center of the corresponding object. This means that in each possible position, the objects are placed in a Bernoulli distributed way with probability equal to the density. In general, the number of objects is not equal to N after this procedure. So if the number of objects is too large, remove objects in low density regions. Otherwise, if the number of objects is too low, refill positions in high-density regions.
Further details on the numerical solution of (1.14) and the update algorithm can be found in sub section 3.1.
In the next section, we present two applications to demonstrate the switching idea. They mainly differ in the evaluation of the force terms in particular the exterior forces.

Material flow model
We briefly introduce the microscopic and macroscopic material flow model. The models are described in detail in Göttlich et al. (2014Göttlich et al. ( , 2015 and the formal derivation of the macroscopic via the hydrodynamic limit is given in Göttlich et al. (2015).

Microscopic model
In the considered application case, circular-shaped parts of radius r are transported on a conveyor belt. The arising exterior forces are the interaction forces between parts and boundaries as well as the friction force induced by the conveyor belt. Analogous to Göttlich et al. (2015), we have to specify the force terms in (1.1) as F being the interaction force between the parts and G the friction force of the belt combined with wall forces. The interaction force F in Equation (1.1b) can be split into a tangential part f t and a normal part f n : ( 2:1) where H denotes the common Heaviside function with HðÁÞ ¼ 1 if k x k < 2r and 0 otherwise. The normal force component itself can be split into an elastic repulsive part f n el and a dissipative part f n diss : The repulsive force is modeled by a spring damper model with the normal spring constant k n while the dissipative force is equipped with the normal viscous coefficient γ n . The tangential force is given by with Coulomb friction coefficient μ and tangential viscous damping coefficient γ t . The tangential The friction force f fric describes the interaction between a part and the surface of the conveyor belt The interaction force between the parts and the boundaries of the conveyor belt are modeled in the same way as the interaction force of the parts. The set containing the indices of walls and boundaries is called W.
The exterior forces are the sum of the interaction forces with the boundaries and the friction forces between the parts and the conveyor belt with f n w;i and f t w;i the normal and tangential forces exerted by the walls on part i.
The complete microscopic material flow model is obtained by inserting these force terms into the general microscopic model (1.1) and reads with F and G as in (2.1), (2.2) and initial conditions

Macroscopic model
The static space-dependent velocity gðxÞ in Equation (1.7) represents the friction force of the conveyor belt and the force exerted by the boundaries of the given setting, see Figure 2a. The resulting velocity field is shown in Figure 2b. The dynamic velocity fðρÞ is needed to deal with situations of maximal density ρ max .
For the formulation of the macroscopic model as in (1.9a) we set where the static velocity field v stat is as in Figure 2b. The nonlocal dynamic vector field v dyn should be active if the density exceeds a maximal density ρ max and parts start to pile up. Below the maximal density, v dyn is inactive, i.e. free flow regime, and the parts are transported with the velocity given by v stat . Thus, the dynamic vector field includes the interaction between parts weighted by the factor > 0 and models the formation of congestion in front of an obstacle by the collision operator IðρÞ: v dyn ðρÞ ¼ Hðρ À ρ max Þ Á IðρÞ; (2:5) where H denotes the common Heaviside function and IðρÞ ¼ À Ñðη Ã ρÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1þ k Ñðη Ã ρÞ k 2 2 q : (2:6) The negative of the gradient of the convolution means that if there is a pile of maximal density, the mass will be pushed in a region with lower density. Inside a pile with maximal density, this force will be zero. Additionally, we have the initial condition ρðx; 0Þ ¼ ρ 0 ðxÞ; x 2 R 2 (2:7) and the boundary condition  hu Ã ðx; ρÞ; nðxÞi ¼ 0; x 2 @Ω W : The computational domain is denoted by Ω and the boundaries are @Ω. The inflow boundary at the left in Figure 2 is called @Ω inflow . The walls and the obstacle are denoted by @Ω W , even if the obstacle is not at the boundary of the computational domain. The normal vector is called n(x).

Microscopic model
In the following sections, the pedestrian model derived by Helbing (Helbing et al., 2000) is briefly reviewed and a corresponding macroscopic model is stated. We assume the pedestrians to be of circular shape with radius r and to act in the same way in each direction, meaning that the forces are rotationally invariant. We specify again the force terms F for the interaction and G for the desired velocity of pedestrians as well as the forces exerted by the walls. We set F ij to characterize the interaction force between pedestrian i and j in (1.1b). The latter is modeled as with hðyÞ ¼ HðyÞ Á y and H the Heaviside function. The constants A; B; k are positive, whereas κ can be equal to zero if tangential force effects are neglected. The distance of pedestrians i and j is given by d ij ¼k x i À x j k with x i and x j the centers of mass, n ij ¼ ðx i À x j Þ=d ij is the normalized vector pointing from pedestrian j to pedestrian i, t ij ¼ ðÀn ij Þ is the tangential direction and Δv t ji ¼ ðv j À v i Þ Á t ij .
The interaction force F iw between wall w 2 W and pedestrian i is modeled as with d iw being the distance of pedestrian i to wall w, n iw the normalized normal vector of wall w and t iw the tangential vector of wall w. Combined with the component of the desired velocity where τ is the reaction time and v 0 ðxÞ is the desired velocity depending on the pedestrians location, we set the exterior force as Compared to the material flow model, where the main ingredient is the friction of the conveyor belt, the exterior force for the pedestrian model is mainly driven by the desired velocity of pedestrians. Note that this leads to a different numerical treatment within the switching algorithm. Inserting these terms into the system (1.1), we have the scaled microscopic pedestrian model

Macroscopic model
For the macroscopic pedestrian model, we have to determine the parts of the velocity u Ã ðρ; xÞ ¼ f ðρÞðxÞ þ gðxÞ, see (1.7), as þ khð2rÀ k x À y kÞ xÀy kxÀyk ρðyÞdy; (2:11) v 0 ðxÞ is the desired velocity depending on the pedestrians location and where h is defined before. The constants A; B; k and the radius of the pedestrians r are the same as in the microscopic model. Here, dðxÞ denotes the distance of x and the interacting wall.

Numerical discretization
In this section, the numerical discretization of the microscopic (1.1) and the macroscopic model (1.9) are presented. The discretized models are then used to set up the hybrid switching method introduced in subsection 1.2.
We set dt the time step of the discretization. The discretized microscopic model reads (3:1b) together with the initial condition It is solved via a classical Runge-Kutta method of fourth order. If the specific application is known, a priori estimates can be performed for the force terms F and G. Then, the time step can be chosen appropriately and the need of an adaptation of the step size during the computation can be avoided.
For the discretization of the macroscopic model we denote by Δx 1 , Δx 2 the space step sizes and Δt the time step size. The grid points are defined as ( 3:2) The space is subdivided into squares S i;j ¼ x We apply a finite volume scheme, where the density is considered as a piecewise constant function ρðx; tÞ ¼ ρ k i;j 2 R for x 2 S i;j ; k 2 1; . . . ; N t f g : Due to (1.7), we have that the velocity parts f ðρÞ and gðxÞ only depend on ρ or x, respectively. This allows for a separate numerical computation in terms of the numerical fluxes i; j ðρ k iÀ1; j ; ρ k i; j ; x iÀ 1 2 ; j Þ (3:3b) where the density-dependent part is given by i;j ρ k i;j ; ρ k iþ1;j ; x iþ 1 2 ;j ¼ ρ i;j g ð1Þ ðx iþ 1 2 ;j Þ g ð1Þ ðx iþ 1 2 ;j Þ ! 0 ρ iþ1;j g ð1Þ ðx iþ 1 2 ;j Þ g ð1Þ ðx iþ 1 2 ;j Þ 0: ( We remark that f and g depend on the application we have in mind. In vertical direction, the computation of the numerical fluxes F i;j and G i;j follows analogously.
A dimensional splitting method is used to reduce the original two-dimensional problem into the solution of one-dimensional problems only. The discretized system is theñ

Hybrid switching algorithms
The numerical discretizations are now used to present an implementation of the hybrid switching method, see subsection 1.2.

Microscopic to macroscopic representation
As we have seen, the given geometry is partitioned into grid cells for the numerical computation. Let us denote by x c and y c vectors containing the x-resp. y-coordinate of the grid cell centers.
Corresponding to (1.13), the density in grid cell i can be calculated via with D j being the density corresponding to object j, see Equation (1.12). The variance σ 0 > 0 has to be chosen such that the numerical support of the bell function is larger than the actual domain of one object, otherwise we obtain peaks for densely packed objects and not a uniform-density distribution. In doing so, the implementation of the microscopic to macroscopic direction is straightforward as expected.

Macroscopic to microscopic representation
The optimization problem (1.14) can be interpreted as a nonlinear least squares problem if the constraints (1.14c) and (1.14d) are omitted (cf. discussion in subsection 1.2) and treated in an additional routine. Therefore, the computation of the macro to micro switching direction can be split into the following steps: (1) Find admissible start solution ρ app : The pseudocode to find an admissible start solution is specified in algorithm 1. The algorithm can be run several times to improve the configuration of the start solution since a random choice of object placement is used.
(2) Solve nonlinear least squares problem (1.14a)-(1.14b): The nonlinear least squares problem is solved with the MATLAB function lsqnonlin https://de.mathworks.com/help/optim/ug/ lsqnonlin.html and requires the solution P of algorithm 1 as initialization. The numerical solution of the nonlinear least sqaures problem is denoted by Q.
(3) Solve problem including the constraints (1.14c)-(1.14d): The main idea of the so-called update algorithm 2 is to move the objects i with coordinates PðiÞ from the optimized solution QðiÞ such that they do not overlap any more. This means in particular if object i interacts with object jÞi or with the wall, i.e. jSðiÞ À SðjÞj 2 2r; jÞi or QðiÞ‚Ω ad , we set SðiÞ ¼ ðPðiÞ þ SðiÞÞ=2 and check again for interactions until PðiÞ À SðiÞ j jis less than a predefined tolerance δ. The final configuration satisfying (1.14c)-(1.14d) is denoted by SðiÞ.

Algorithm 1 Find admissible start solution
Input: Density ρ mac , mass m, radius r, # objects N, tolerance , admissible region Ω ad Output: Solution P containing the coordinates of the object centers 1: Compute a hexagonal grid structure with grid point distance 2r and grid points g i : = (x i , y i ) 2: for i = 1: # grid points do 3: if g i ∉Ω ad then 4: Delete grid point g i 5: end if 6: end for 7: for i = 1: # remaining grid points do 8: Determine corresponding density value ρi = ρ mac (x i , y i ) 9: if ρ i < then 10: Delete grid point i 11: else 12: Delete grid point i with probability (1 − ρi) 13: end if 14: end for 15: if N > # remaining grid points then 16: Put grid point with index argmaxi{ρ i |g i has been deleted} until N is met 17: else if N < # remaining grid points then 18: Delete grid point with index argmini{ρ i |g i not deleted} until N is met 19: end if 20: return Remaining grid points as matrix P corresponding to the placed objects centers  the MATLAB algorithm lsqnonlin, see Figure 3c, where we directly observe that the nonoverlap condition is injured, and a figure for the final solution after applying the update algorithm 2, see Figure 3d. Concluding, the latter is a combination of the start and optimized solution. We observe that in regions, where the density is high, the objects are placed as dense as possible, whereas in regions, where the density is far from the maximal density, objects get stuck.
Note that the runtime growths linearly in N for algorithm 1 and the solution of the optimization problem. The runtime for the update algorithm 2 growths cubically in N. However, the solution of the optimization problem using MATLAB takes several hours and dominates the performance of the update algorithm for the test cases we present in Section 4. if sðiÞ‚Ω ad then 7: interaction = 1 8: end if 9: for j = 1: N do 10: if |S(i) − S(j)| 2 ≤ 2r, j ≠ i then 11: interaction = 1 12: end if 13: end for 14: if interaction = 1 then 15: if |P (i) − S(i)| < δ then 16: STOP, S(i) = P (i), go back to step 3 17: end if 18: S(i) = (P (i) + S(i))/2 and go back to step 6 19: end if 20: end for 21: end for 22: return Final solution S

Numerical results
The application examples in Section 2 are now analyzed from a numerical point of view. For each scenario, implementation details on the model are given and numerical examples are explained in detail.

Results for the material flow model
The numerical results of the switch are investigated for a setting as presented in Figure 2. This experimental setup has been originally introduced in Göttlich et al. (2014), where the microscopic and macroscopic material flow models have been validated against real data. The idea is now to use the same data to test the performance of the switch.
As described in Section 3, the simulation of the microscopic model is done with a classical Runge-Kutta method of fourth order which is applied to the discretized microscopic model (3.1). The discretization of the force terms stated in Section 2.1 is straightforward. The parameters are characterized by the bottom viscous damping coefficient γ b ¼ 0:5, the normal viscous coefficient γ n ¼ 0:1, the viscous damping coefficient γ t ¼ 0:16, the normal spring constant k n ¼ 200, the bottom friction coefficient μ b ¼ 0:5, the mass m ¼ 0:01, the radius r ¼ 0:012, the gravitational constant g ¼ 9:81, the Coulomb friction coefficient μ ¼ 0:17, the time step size dt ¼ 0:005 and time horizon T ¼ 4. The number of parts is N ¼ 192: Considering the macroscopic material flow Equations (2.8), we have to evaluate the numerical fluxes for the static velocity and the dynamic velocity, see (2.4), according to (3.3), where the first component of the interaction operator I (cf. (2.6)) is discretized as being a grid cell. The discretization of the second component with respect to the x 2 -direction is computed analogously. We refer to Göttlich et al. (2014) for more details. We choose a Gaussian convolution kernel η of type As grid sizes of the macroscopic model, we have Δt ¼ 0:0025 and Δx 1 ¼ Δx 2 ¼ 0:008. The final time horizon is again T ¼ 4 and the switch is evaluated at time t switch ¼ 1: To validate the results of the switching procedure from micro to macro, the evolution of parts including the switch of the system is compared to real data, cf. . We measure how many mass or parts leave the system behind the obstacle, i.e. the outflow over time. The corresponding diagram is shown in Figure 4.
The maximal outflow error, meaning the maximal difference of the outflow of the model including the switch and the real data outflow, is about 6%. The maximal outflow error of the macroscopic model without any switch is approximately 5% (see ), and thus indicates that the switch performs quite well.
The reverse situation, i.e. from macro to micro, is depicted in Figure 5. Here, the left picture shows the particle distribution resulting from the update algorithm 2. Even if the switch from the macroscopic to the microscopic representation is the more difficult direction, the maximal outflow error is now smaller compared to Figure 4. This is due to the higher accuracy of the microscopic model. The maximal outflow error, again compared to the real data outflow, is about 3:8%.

Results for the pedestrian flow model
We consider an experimental setting with a room which has to be evacuated. Figure 6a shows the room geometry as well as the static velocity consisting of the desired velocity v 0 ðxÞ and the wall forces ∑ W F W ðxÞ, see Figure 6b.
For the numerical implementation, almost the same parameters as in Helbing et al. (2000) are used. We choose the number of pedestrians as N ¼ 81, the mass m ¼ 80, the radius r ¼ 0:25, the desired velocity v 0 ¼ 1, the acceleration time τ ¼ 0:5, the constants A ¼ 2000; B ¼ 0:08 and the parameter k ¼ 1:2 Á 10 5 . Since during the derivation, the tangential forces are neglected, we set κ ¼ 0 to directly compare the microscopic and macroscopic results. The final time horizon is T ¼ 50 and the time step to solve the Runge-Kutta scheme is dt ¼ 10 À4 .
The discretization of the force terms stated in Section 2.2 is straightforward, see (3.1) such that a classical Runge-Kutta method of fourth order is used again.  (a) Static velocity field (b) Velocity field induced by the walls Figure 6. Evacuation of a room with one exit.
Corresponding to the discretization of the macroscopic Equation (3.4), the discretized version of the flux terms (3.3) is being a grid cell. The function h is defined as before.
For the external force term, it is We choose nearly the same parameters as in the microscopic case, i.e. only the parameters A ¼ 1; B ¼ 0:1 and k ¼ 27 have to be varied, to match the results of the pure microscopic simulation. As grid sizes, we have Δx 1 ¼ Δx 2 ¼ 0:1 and Δt ¼ 0:001.
Similar to the switch for the material flow model, a switch is implemented for the pedestrian model at time t switch ¼ 10: To compare the qualitative behavior of the models, the outflow at the door is measured. Since there is no real data available this time, the outflow of the model with switch is compared to the outflow without any switch, contrary to the material flow case. Figure 7 shows the density distribution from the microscopic to macroscopic representation and the outflow behavior over time. As we can observe, the outflow of the switched model is close to the simulation results obtained from the microscopic simulation only. The maximal outflow error is about 5:3%.
The results for the macro to micro switch can be found in Figure 8. Considering the outflow diagram, the results of the switched model are compared to the macroscopic simulation. The maximal outflow error is about 3:2%, which is within an acceptable range.

Computation times
Finally, we present the computation times for our switch implementation. All computations have been performed on a PC equipped with 16 GB Ram, Intel(R) Core(TM) i5-4690 CPU@3.50 GHz. Tables  1 and 2 show the computation times for both switching directions (micro to macro and vice versa) and the two applications material and pedestrian flow. The same discretizations as introduced formerly are used. However, we remark that in the macro to micro direction the radius needs to be adapted in order to match the total mass (1.10).
The times for the three steps in the macro to micro direction are listed separately. As we can see, the most time-consuming part is the solution of the optimization problem (1.14) which takes several minutes up to almost one hour for N ¼ 200 objects, whereas the computation of the start solution (algorithm 1) as well as the update algorithm 2 takes only a couple of seconds.
The total computation time increases exponentially with the number of objects. Nevertheless, it might happen that an optimum is found quite fast, e.g. for N ¼ 150 pedestrians in Table 2. Note that the total computation times differ due to the different discretizations and geometries for the material and pedestrian flow models, cf. discussions in sub sections 4.1 and 4.2.
The computation time of the microscopic model (without any switch) increases quadratically in the number of objects N . Therefore, the simulation of the macroscopic model  provides a less costly alternative and good results. A switch to the microscopic representation is for instance useful if the user is interested in a better resolution of the underlying dynamics, even if the optimization step takes some time.