A Harmony Search Algorithm for the Reproduction of Experimental Data in the Social Force Model

Crowd dynamics is a discipline dealing with themanagement and flow of crowds in congested places and circumstances. Pedestrian congestion is a pressing issue where crowd dynamics models can be applied. The reproduction of experimental data (velocitydensity relation and specific flow rate) is a major component for the validation and calibration of such models. In the social force model, researchers have proposed various techniques to adjust essential parameters governing the repulsive social force, which is an effort at reproducing such experimental data. Despite that and various other efforts, the optimal reproduction of the real life data is unachievable. In this paper, a harmony search-based technique called HS-SFM is proposed to overcome the difficulties of the calibration process for SFM, where the fundamental diagram of velocity-density relation and the specific flow rate are reproduced in conformance with the related empirical data. The improvisation process of HS is modified by incorporating the global best particle concept from particle swarm optimization (PSO) to increase the convergence rate and overcome the high computational demands ofHS-SFM. Simulation results have shownHS-FSM’s ability to produce near optimal SFMparameter values, whichmakes it possible for SFM to almost reproduce the related empirical data.


Introduction
Massive congestion in pedestrian environments has been a pressing concern lately.This issue stems from an increase in population growth rate, where among the serious problems linked to it are crowd stampedes, accidents, and other disasters.Congestion has therefore motivated researchers to find solutions, where studies have focused on offering good pedestrian facilities or understanding incorrect pedestrian behavior and correcting them [1].
May [2] basically divides pedestrian dynamics studies into two levels.The first level is the macroscopic studies, which is concerned with macroscopic behaviors of the whole crowd.Macroscopic variables resulting from these behaviors are density (), mean speed (), and pedestrian flow () [3].The second level is the microscopic studies, which is more concerned with the detailed interactions among the pedestrians, such as avoiding collisions, deviations, acceleration and deceleration, and their effect on motion [1].Compared to the macroscopic pedestrian studies, the consideration of the detailed pedestrian interactions in the microscopic studies enables the researchers for proposing models characterized with more comprehensive pedestrianflow control and better movement-quality improvement [1].Famous examples of microscopic crowd dynamics models are the cellular automata model [4][5][6][7][8] and the social force model (SFM) [9][10][11][12][13][14][15][16][17][18][19].
For the purpose of validation of such microscopic crowd dynamics models, a variety of methods was adopted.These methods include the introduction of the self-organization phenomena and the reproduction of experimental data obtained by experimental and empirical studies.The experimental data, for instance, are the specific flow rate (number of persons crossing an opening per unit of time and width) and the fundamental diagrams (a graphical representation as shown in Figure 1 used to describe the relationship between Density (p/m 2 ) Figure 1: The diagram for the estimated speed-density relations for unidirectional pedestrian traffic flow in empirical studies.Fruin [20], Sarkar and Janardhan [21], and Older [22] used a linear speeddensity relation, while Weidmann [23] used a double S-bended curve.
two of macroscopic variables mentioned above, such as the mean speed and the density ( = ()).The fundamental diagrams of unidirectional and bidirectional streams (among which are Fruin [20], Sarkar and Janardhan [21], Older [22], and Weidmann [23]) form essential tools for the assessment of the models; whether, it can describe the pedestrian stream appropriately with respect to the empirical studies.The main concern is the calibration of these models.A common method of the calibration process involved with the models' assessment is to optimize some model parameters and model components long enough; until, the simulations fit the fundamental diagram, such as that of Weidmann [23], or fit the specific flow rate such as that proposed by Parisi et al. [24].
In this paper, we focus our attention on proposing new technique for calibrating essential parameters in the social force model (SFM).The model has been rendered as one of the most vital models in microscopic studies because of its successful introduction of self-organization phenomena of pedestrian dynamics in normal and panic situations [14][15][16].The most important feature of the SFM is its representation of pedestrians' motivations in terms of the other objects (pedestrians and obstacles) surrounding them as social forces.The sum of these forces is implemented in a Newtonian equation, which, in turn, determines the acceleration of the pedestrian's motion.
In [18,19,[24][25][26][27], calibrating the SFM parameters A rep , B rep , and , which represent the strength of the social forces, the repulsive distance range of these forces, and the angular perception parameter, respectively, is a challenging task for the reproduction of the fundamental diagram.Johansson et al. [25] proposed an evolutionary optimization algorithm to determine optimal parameter specifications for the SFM.Zainuddin and Shuaib [26] adopted a different approach where the model of the parameter B rep is introduced as a function in terms of density.However, although the reproduction of the fundamental diagram was successful in [25,26] without changing the principles of the SFM, the values of the adopted parameters do not help in reproducing the flow rate in normal evacuation situations within the range that lied from 1.25 to 2 p/m/s as stated in [24].
In this paper, a harmony search-based optimization technique named HS-SFM is proposed to overcome the difficulties of the calibration process for the SFM, wherein the fundamental diagram of velocity-density relation and the specific flow rate are reproduced in conformance with the related empirical data.Harmony search (HS) is a relatively new population based on metaheuristic algorithm, which has obtained excellent results in the field of combinatorial optimization [28].With the intention of increasing the convergence rate and overcoming the high computational demands of the proposed HS-SFM algorithm, the memory consideration, which is the most important operator in the improvisation process of HS, is modified.The global best particle concept from PSO is used as a new memory consideration operator to mimic the best harmony among the HM vectors, which in turn is used to construct the new harmony.In addition, a new multicriteria objective function is proposed to include the abovementioned two main factors that affect the goodness of the solution with respect to SFM.
The organization of this paper is described as follows.Section 2 provides an introduction to the social force model.Section 3 introduces readers to the harmony search algorithm.A detailed description of the proposed HS-SFM algorithm is provided in Section 4. The simulation study of the proposed algorithm is presented in Section 5. We conclude our findings in Section 6, including remarks regarding future works.

The Social Force Model
The representation of the pedestrians' motivations as social forces implemented in a Newtonian equation for introducing the pedestrians' motion is an essential feature of the SFM.An extension by Helbing et al. [13] incorporated the physical forces arising in the case of contact among the pedestrians into the model.The model is characterized by reproducing the self-organized phenomena of pedestrian dynamics in normal and panic situations [13][14][15][16].
The main equations of the model are ), which were modeled as linear functions in [13,14], in analogy with the granular forces; ⃗  ,object () is analogous with ⃗   () but with regard to objects such as walls and columns.
The model of the repulsive social force is where A rep and B rep are constant parameters that denote the strength and the repulsive distance range of the corresponding force, respectively; ⃗   is the normalized vector which points from the object  to the pedestrian ;   and   are the sum of the radius of  and  and the distance between the centers of  and , respectively;   () denotes the angle between the direction ⃗   () = ⃗ V  ()/‖ ⃗ V  ()‖ of motion and the direction − ⃗   , (i.e., cos   () = − ⃗   ⋅ ⃗   ); and  is an angular perception parameter to model the effect of the perception of individual  to those who are behind him on the magnitude of the force.Its value belongs to the range [0, 1] and specifies the shape of the anisotropic force field (territorial area) around individual .

Harmony Search Algorithm
Geem et al. [28] developed HS to solve optimization problems.Ever since its inception, HS has attracted many researchers to develop HS-based applications for a variety of optimization problems [29].HS imitates the natural phenomenon of musicians' behaviors when they collectively tune the pitches of their instruments to achieve a fantastic harmony as measured by aesthetic standards.It is an effective metaheuristic algorithm that can explore the search space of the given dataset in a parallel optimization environment, where each solution (harmony) vector is generated through intelligent exploration and exploitation of a search space.It is thus considered a population-based algorithm with local search-based aspects.The following provides a description of HS.The standard HS algorithm consists of five steps as follows.
Step 1 (HS initialization and optimization problem parameters).In the case of HS, optimization is treated as a minimization (or maximization) problem.Basically, minimize (or maximize) () subject to   ∈   ,  = 1, 2, . . ., , where () is the objective function;  is the set of decision variables (  );  is the set of the possible range of any decision variable,   = {  (1),   (2),   (3), . . .  (), }; () is the number of possible ranges for each decision variable; and  is the number of decision variables.These parameters of the HS are then initialized, which form the harmony memory size (HMS), harmony memory consideration rate (HMCR), pitch adjustment rate (PAR), and the number of improvisations (NI).
Step 2 (harmony memory initialization).The harmony memory (HM) is a matrix of solutions with size HMS, as shown in (8).In this step, the solutions are randomly constructed and rearranged progressively according to their objective function values Step 3 (improvise a new harmony).In this step, the HS generates (improvises) a new harmony vector, ), based on the three operators of memory consideration, random consideration, and pitch adjustment.In memory consideration, the values of the new harmony vector are randomly inherited from the historical values stored in the HM with the probability of HMCR.The value of the decision variable  NEW

𝑖
) based on the historical values stored in the HM.The usage of the HM is similar to the step where the musician uses his or her memory to generate an excellent tune.This important step ensures that good harmonies are considered as the elements of the new harmony vectors.Since it is a cumulative process, other values are not chosen depending on memory consideration but are chosen according to their possible ranges through random consideration with a probability of 1-HMCR.This step is called randomization, which increases the diversity of the solutions and drives the system further to explore various diverse solutions so as to attain the global optimum.Furthermore, the new vector components that are selected out of the memory consideration operator are examined to be pitch adjusted with probability of PAR ∈ [0, 1] as The variable  is an arbitrary distance bandwidth used to improve the performance of HS.The PAR parameter simulates the music by "changing the frequency", which means generating a slightly different value for the new harmony vector component.This parameter explores more solutions in the search space for that purpose.

𝑁
), replaces the worst harmony in the HM, only if its fitness value (measured in terms of the objective function) is better than that of the worst harmony.
Step 5 (check the stopping criterion).Termination occurs when the maximum number of improvisations NI is reached.
The following section describes the proposed HS-SFM algorithm.

The Harmony Search-Based Social Force Model (HS-SFM)
This section presents the proposed HS-SFM algorithm in detail and shows how HS is implemented to find the optimal values of SFM's parameters.

Harmony Memory Initialization.
Each harmony memory vector represents a candidate of SFM parameter values.To initialize the HM with feasible solutions, each parameter value is randomly generated from its valid range.In this study, the valid ranges are as follows: (i) strength of the social repulsive force ( rep ): [0, 3000]; (ii) repulsive distance range of the corresponding force ( rep ): [0, 2]; (iii) angular parameter (  ): [0, 1].
For example, the vector [1250, 0.67, 0.45] is a candidate harmony memory vector that can be randomly generated from the given ranges.
In order to discover the goodness of the generated harmony vector with regard to SFM, the fitness function for each harmony vector is calculated (as will be discussed in Section 4.3) and saved in the harmony memory.After the harmony memory is constructed and the fitness function for each vector is calculated, the HS-FSM algorithm initiates an iterative search for the optimal values for the given SFM parameters.This search is governed by improvisation rules, which are described in the following section.

Improvise a New
Harmony.The new harmony vector is a new candidate SFM parameter value vector, and the values of this new harmony vector is generated depending on the HS's improvisation rules.In this study, however, memory consideration, which is the most vital operator in HS, is modified.In the standard HS, most of the variables in the new harmony are randomly selected from the other vectors stored in HM.In this study, we adopt the supporting mechanism concept from PSO [30] to select a promising value for the variables from the vectors stored in the HM.This supporting concept mimics the best harmony among the HM vectors to construct the new harmony.The inspiration behind this modification is so that HS-SFM converges faster as compared to the standard HS algorithm.In standard HS, convergence to the optimal parameter values is smoothly achieved.In this scheme, the number of simulations needed to reach the optimal state of HS-SFM is reduced.It is worth noting that although having different convergence rates, this modification does not affect the quality of the solutions for both the HS and HS-SFM algorithms [31].
Based on this modification, the values of the new harmony vector components will be selected from the best harmony memory vector stored in HM, with the probability of HMCR.On the other hand, the value of the components of the new harmony vector is selected from the possible range with a probability of 1-HMCR.The following equation summarizes these two steps, that is, memory consideration and random consideration Furthermore, the new vector components selected out of the memory consideration operator are examined to be pitch adjusted with the probability of PAR as in (11).In this study, we used the standard process of the pitch adjustment in HS as in (9).The bandwidth parameter, , is responsible for the slight modification of the selected component, which will hopefully achieve better solutions by deeply exploiting the search space.Consider Next, the fitness function is computed for the new generated harmony memory vector as described in Section 4.3.Then, the new vector is compared with the worst harmony memory solution using a fitness function.If it is better than the existing solution, the new vector is included in the harmony memory and the worst harmony is discarded.This process is repeated until the maximum number of iterations NI is reached.Finally, the best solution among the maximum value of fitness function of each HM solution vectors is selected as the best solution vector.

Evaluation of Solutions.
In order to measure the goodness of each harmony memory vector, which in turn measures the goodness of the SFM parameter values, a new fitness function is proposed that includes the two main factors that affect the goodness of the solution with respect to SFM.These factors are (1) the reproduction of the fundamental diagram and (2) the reproduction of the flow rate.
The basic model of our proposed fitness function is where where the function  fundemental denotes the absolute deviation between the mean velocity V sim () resulting from our simulation and the corresponding mean velocity V weid () resulting from (15), which is an estimation by [23].The density  in ( 15) is calculated as the number of pedestrians within a specific area at any given moment.The optimal value for  fundemental is zero, which means that the mean velocity V sim () produced from our simulation fits the corresponding mean velocity V weid () produced from (15).Therefore, the improvisation process in each iteration tries to reach the optimal value for this factor, which is zero.The measurement flow rate is defined as the number of pedestrians passing through a section per unit time.The allowed range for a specific flow rate is [low flow, up Flow], which is determined according to the relevant empirical studies mentioned in [24].Specifically, these studies state that the various values of flow rates in normal evacuation situations lie within the range of 1.25 to 2 p/m/s.According to (14), the flow rate produced from our simulation that falls within the optimal range results in  evacuation to be equal to zero, which is the optimal value for this factor.Accordingly, this also indicates that the generated harmony vector can produce a good solution for SFM and may also reproduce the fundamental diagram of [23].Resultantly, the minimization of  evacuation and  fundemental to reach the optimal value of  HS-SFM , which is equal to zero, is desired and any harmony vector solution that can reach this value is considered as the optimal solution for the HS-SFM algorithm.The pedestrians' radius: uniformly distributed within the range [0.25, 0.30] m The SFM's parameters  = 1.2 × 10 5 kg/sec 2  The strength of the contact (pushing) force  = 2.4×10 5 kgm −1 s −1 The strength of the friction V 0 = 1.34 m/s The preferred velocity  = 0.5 s The reaction time The

Simulation
Setting.This section presents the simulation results of the proposed HS-SFM algorithm.All simulations were conducted using MATLAB version R2010a.Each simulation comprises two scenarios (or parts), which were conducted simultaneously.The physical environment in the first scenario used for the reproduction of the fundamental diagram is that used in [26,32].The positions of all simulated pedestrians were simultaneously randomly initialized in the simulated area, which represents a horizontal unidirectional walkway, and their motions were directed horizontally towards their destination.Analogous to the characteristics of pedestrians in the experimental studies performed by Weidmann [23], all pedestrians had the same preferred direction, with preferred speeds varying between 0.97 m/s and 1.65 m/s.For the reproduction of the specific flow, the second scenario was for the evacuation process.Similar to the simulations performed by Parisi et al. [24], the physical environment was a 20-meter × 20-meter room, with one, 1.2 width, exit.The simulated pedestrians were randomly initialized inside the room with the verification that there was no overlap between any two pedestrians or a pedestrian with a wall.The simulated pedestrians evacuated from the environment were 100, 200, 300, and 400.Table 1 shows the values of pedestrians' parameters and the parameters of the HS-SFM used in this simulation.

Simulation Results.
The validation of each simulation being done by HS-SFM algorithm is done by a comparison between those results obtained from HS-SFM algorithm and those obtained from the fundamental diagram, developed by Weidmann [23], and the value of the flow rate for an evacuation process, obtained by Parisi et al. [24].
The simulations were run 400 times, and a snapshot of the results is shown in Table 2.The highlighted rows in Table 2 show that our HS-SFM algorithm can find a near optimal SFM parameter value and can almost reproduce the fundamental diagram as estimated by [23].Furthermore, the flow rate produced by our simulation is within the optimal range as stated in [24].
For further clarification of the results obtained from the HS-SFM algorithm and the results shown in Table 2, an arbitrary case can be considered and described.The case is simulation number 38, where the optimal values obtained by HS-SFM are  rep = 230.85, rep = 0.67, and  = 0.76, and the fundamental diagram obtained from these parameters is shown in Figure 2 with a comparison with the fundamental diagram estimated by [23].It can be seen that in most densities (i.e., 1, 2, 3, and 6), HS-SFM produced a fundamental diagram as in [23].However, some differences between the simulated and Weidmann fundamental diagrams appeared when the density is equal to 4 and 5.This affected the objective function  fundemental , whose result was equal to 0.07.Our future work is to study how to improve the performance of HS-SFM to overcome this shortcoming.For the flow rates produced for this case study, where the number of pedestrian varied from 100 to 400 pedestrians in a room/simulation, all of them are in the optimal range as stated in [24] and consequently the value of  evacuation was equal to zero.In total, the HS-SFM objective function was equal to  HS-SFM = 0.07 + 0, which is an acceptable solution.We can summarize our findings of the optimal SFM's parameter values by setting an optimal range for each parameter as seen in Table 3. Experimentally, the optimal range for the  rep is set to [120, 235],  rep is set to [0.41, 0.98], and  is set to [0.53, 0.76].
It is important to note that firstly the resulting fundamental diagram shown in Figure 2 conforms to the diagrams shown in Figure 1, which exhibits the same linear shape.Secondly, in contrast to the results of the computer simulation models, with these optimal vector values, the mean speed of the simulated pedestrians in the case of jam density ( > 4) does not vanish.Such behavior is in conformance with the experimental and relevant field studies stated by Zhang et al. [33].Thus, our result justifies the reliability of our approach for the calibration of the parameters.
Thirdly, the range of the  rep obtained in our study is consistent with the literature, that is, [10,13,18].For the purpose of comparison, in [25], the best fit curve to the Weidmann fundamental diagram was established when the repulsive distance range parameter has a significantly high value,  rep = 3.22.Such a high value would cause social forces between pedestrians to enforce longer distances between their locations.As a result, while performing simulations for the evacuation process with the consideration of this value, the densities of the clogging areas on exits were lesser as compared to those stated in the literature and observed in reality.This result is a main reason for a low specific flow rate.
Fourth, the repulsive distance parameter in [26] was modeled as a function of density, and thereby, the values of the parameter were more consistent with the literature, that is, [10,13,18].However, the value of the angular parameter was significantly high,  = 0.9, which granted the pedestrians Speed (m/s) almost full perception for those who were behind.From psychological point of view, such results are unrealistic as pedestrians (who have high behind perception) are exposed to pushing social forces from the pedestrians behind them.Therefore, the pedestrians become competitive as described in [24].Due to this competitiveness, the specific flow rate is higher than the range stated in the literature.In our results, the range of the angular perception,  ∈ [0.53, 0.76], is significantly moderate, which results in lower flow rate than what was reproduced in [26].

Conclusions and Recommendations
The SFM is considered one of the most distinguished representative microscopic dynamics models that present solutions to pedestrians' congestion problems.In order to improve its performance, we proposed the HS-based algorithm HS-SFM that offers near optimal parameter values for the SFM.Simulations conducted showed the ability of HS-SFM to find near optimal parameter values that can help in the reproduction of both the fundamental diagram as in [23] and flow rate as in [24].
Our future work will focus on improving the ability of HS-SFM to reduce the difference between the mean velocity V sim () resulting from our simulation and the corresponding mean velocity V weid () resulting from (15).The second direction of our future work will focus on how to decrease the computational time required by HS-SFM.
As the microscopic dynamics models have been exposed to many amendments, such as the incorporation of new submodels or refining the models itself, the need for a technique that can offer near optimal parameter values is necessary.We foresee that more effort for the development of these techniques is necessary by studying and potentially improving their behaviors in exploring the search space of the dynamics models parameter values.We also recommend the hybridization of techniques such as harmony search with other metaheuristic optimization algorithms in order to improve search capabilities to find the optimal parameter

4 , 𝑎 NEW 5 ,
. ..), are consecutively chosen in the same manner within the probability of HMCR ∈ [0, 1].The HMCR parameter is the probability of selecting one value from the decision variable, ( NEW fluctuation source of the pedestrian's acceleration is randomly assigned to each individual The HS-SFM parameters HMCR = 0.95 Harmony memory consideration rate HM = 20 Harmony memory size PAR = 0.75 Pitch adjustment rate NI = 500 Number of improvisation [10]repulsive and attractive motivations inside pedestrian  against and towards pedestrian , respectively[10]; they were modeled as exponential functions with different values of the parameters and opposite directions.The second type is the physical forces ( ⃗