Automatic model construction for the behavior of human crowds

Applied Soft Computing 56 (2017) 368–378

Contents lists available at ScienceDirect

Applied Soft Computing journal homepage: www.elsevier.com/locate/asoc

Automatic model construction for the behavior of human crowds Jinghui Zhong a,b,∗ , Wentong Cai b , Michael Lees c,d , Linbo Luo e a

South China University of Technology, School of Computer Science and Engineering, Guangzhou, China Nanyang Technological University, School of Computer Science and Engineering, Singapore c University of Amsterdam, Computational Lab, Amsterdam, Netherlands d National Research University ITMO, St. Petersburg, Russia e Xidian University, School of Cyber Engineering, Xi’an, China b

a r t i c l e

i n f o

Article history: Received 3 August 2016 Received in revised form 26 January 2017 Accepted 11 March 2017 Available online 29 March 2017 Keywords: Agent-based modeling Crowd modeling and simulation Genetic programming Gene expression programming Symbolic regression

a b s t r a c t Designing suitable behavioral rules of agents so as to generate realistic behaviors is a fundamental and challenging task in many forms of computational modeling. This paper proposes a novel methodology to automatically generate a descriptive model, in the form of behavioral rules, from video data of human crowds. In the proposed methodology, the problem of modeling crowd behaviors is formulated as a symbolic regression problem and the self-learning gene expression programming is utilized to solve the problem and automatically obtain behavioral rules that match data. To evaluate its effectiveness, we apply the proposed method to generate a model from a video dataset in Switzerland and then test the generality of the model by validating against video data from the United States. The results demonstrate that, based on the observed movement of people in one scenario, the proposed methodology can automatically construct a general model capable of describing the crowd dynamics of another scenario in a different context (e.g., Switzerland vs. U.S.) as long as that the crowd behavior patterns are similar. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Agent-based modeling is a popular and powerful crowd modeling paradigm that has a range of applications in digital game design, military training, crowd evacuation planning, etc. [1–6]. In the agent-based paradigm, individuals are modeled as autonomous agents that can perceive information and make decisions based on individual behavioral rules. A fundamental issue in all agentbased modeling is how to design suitable behavioral rules so that the simulation can match the desired crowd behaviors (e.g., those observed in videos). Once the correct behavioral rules are obtained, the crowd simulation model can be used for prediction and performing “what-if” analysis. However, designing correct behavioral rules to generate a desired crowd behavior remains a challenging task because an agent’s behavior is affected by its environment as well as the surrounding agents. Typically, this is done in a fairly adhoc manner, and relies on expert intuition. However, intuition can often be misleading as crowds are complex systems and as such the relationships between the individual behavioral rules and the collective crowd dynamics are often unintuitive.

∗ Corresponding author. E-mail address: [email protected] (J. Zhong). http://dx.doi.org/10.1016/j.asoc.2017.03.020 1568-4946/© 2017 Elsevier B.V. All rights reserved.

There are a number of existing models that have taken a scientifically Popperian approach to understanding crowd dynamics. A theory is developed which can then be falsified against real data. Helbing et al. [7,8] suggested using virtual force to guide pedestrians’ motions. They designed mathematical equations to formulate the attractive and repulsive forces from goals, other pedestrians and obstacles. Yamaguchi et al. [9] designed behavioral rules using an energy function that considers multiple personal and social factors such as damping and grouping. Similarly, Pellegrini et al. [10], Scovanner and Tappen [11] also designed energy functions by considering different personal and social factors to determine movements of pedestrians. These manual methods usually require a model calibration process to tune coefficients of the initialized formulas, so that the model dynamics can match the desired crowd dynamics [9,12,10,11,13,14]. However, if the initialized formulas are not well-defined, the calibration will not be able to generate the desired simulation results. In addition, the model calibration is also a time consuming and tedious process. Recently, more Baconian, or data-driven approaches, have been proposed to automate the generation of crowd behaviors. The key idea is to learn examples (mostly in terms of state-action pairs) from video data. These examples are then used as behavioral rules in the simulation to generate agents’ movement [15–20]. In the computer vision community, several approaches have been proposed that learn global motion patterns (e.g., certain

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

movement curves) of crowds from video data [21–25]. The global motion patterns can be integrated with a crowd simulation model to generate realistic crowd behaviors. However, the knowledge learned by these methods (e.g., examples and global motion patterns) are typically highly scenario specific where training data is available. Once the application scenario changes, the learned knowledge may become incapable of generating the desired crowd dynamics. In this paper and in previous work [26] we aim to develop a way of merging the Baconian and Popperian philosophies, by developing ways to automate model (or rule) generation. Effectively we allow the computer to find sets of plausible theories (models) that we can read and understand and perhaps manipulate later. In summary the core question of this paper is: “Can we design an automatic method to learn generic behavior rules from video data?” Here the term “generic” means that the behavioral rules are applicable not only to the scenario where training data are obtained from but also to other new scenarios. Also the notion of rule is critical, whereas a pure data-mining approach may be able to derive a model (e.g., trained neural network), the notion of a rule implies that humans can understand and perhaps modify the model. Specifically, this paper focuses on simulating crowd behaviors in an area with one or multiple source regions (SRs) and destination regions (DRs). Each pedestrian in the given area keeps moving toward its own DR and leaves the area via the DR. Such a scenario is quite common in public places, e.g., train stations, airports, etc. We assume that the environment contains obstacles and that pedestrians move toward their DRs casually, following certain curved paths (e.g., “S” and “C” shape curves). The movement curves are determined by specific environment factors such as the positions of obstacles and the destinations. Our goal is to automatically learn a generic behavioral rule (labeled as  ) from video data which describes the way pedestrians move toward their DRs while avoiding obstacles. In the proposed methodology, an abstract dual-layer architecture is used to model individual behaviors. The bottom layer adopts the commonly used social force model (SFM) [8] to model the collision avoidance behaviors of pedestrians. Meanwhile, in the top layer, a behavioral rule that defines velocity fields is used to determine the strategic or goal-directed movement of pedestrians. A velocity field is a set of position-velocity pairs that describe the velocities of pedestrians in the simulation regions [27,28,20]. Based on the dual-layer modeling architecture, the problem of finding a suitable behavioral rule to generate the desired crowd behaviors is formulated as a symbolic regression problem. By using a genetic programming (GP) [29] based method to solve the symbolic regression problem, a behavioral rule that leads to the desired simulation results can be obtained automatically. To validate its effectiveness, the proposed methodology is applied to learn  from video data of a scenario in Switzerland. The learned  is then used to generate crowd dynamics in a new scenario in United Stated. The results have demonstrated that the proposed methodology is capable of distilling suitable behavioral rule that can generate realistic crowd dynamics in the new scenario. This paper extends two of our previous papers on data-driven crowd modeling [26,20]. In the first paper [26], we have conducted a preliminary research on learning goal selection rules of agents from simulation data. This paper extends [26] by addressing a complementary problem – to distill generic behavioral rules to define the velocity fields that can be used to drive agents toward their goals. In the second paper [20], we have proposed a data-driven method to learn behavioral patterns (e.g., the goal selection patterns and the velocity fields) from videos so as to generate realistic crowd dynamics. However, the behavioral patterns learned in [20] are only applicable to the specific scenario where training data is obtained from. In this paper, we aims to distill generic behavior

369

rules that can be applied directly to different scenarios that have similar crowd behavior patterns. In general, the major contributions of the paper are as follows. First, we have formulated the problem of distilling suitable crowd behavior rule from crowd data as a symbolic regression problem. Second, we have developed a genetic programming based framework to solve the formulated problem. The proposed framework is generic such that the solution obtained from one scenario can be applicable to other scenarios that have similar crowd behavior patterns. In addition, two real world datasets have been used to validate the effectiveness and generality of the proposed framework. The rest of the paper is organized as follows: Section 2 describes the related work. Section 3 presents the outline of the proposed framework. Section 4 describes the procedure of applying Self-Learning Gene Expression Programming (SL-GEP) to distill behavioral rules from video data. Section 5 presents the simulation studies. Finally, Section 6 presents the conclusions and discussions.

2. Related work Over the past decades, various methods have been proposed to help automate the process of modeling human crowds. These works can generally be classified into two groups. The first group focuses on automatic parameter calibration. The objective of this class of methods is to use optimization algorithms such as genetic algorithms (GA) to tune parameters of a model so that the simulations match the desired crowd behaviors (e.g., those observed in the video). For example, Johansson et al. [12] proposed an evolutionary algorithm (EA) to calibrate the parameters of the social-force model (SFM) [7]. They utilized the microscopic motions of pedestrians such as the moving speed and direction to evaluate the quality of simulation results. Wolinski et al. [13] suggested using EAs (e.g., GA with a local search) to calibrate crowd models such as SFM and RVO2 based on video data. They defined several metrics (e.g., speeddensity relationships) for simulation evaluation. Pellegrini et al. [10] proposed a dynamic model named Linear Trajectory Avoidance (LTA) model for multi-target tracking. The LTA model consists of an energy function with six parameters to estimate the next velocities of pedestrians. A GA is used to optimize the parameters of the LTA model based on video data so that the prediction errors of the LTA model are minimized. Yamaguchi et al. [9] proposed an agent-based behavioral model to improve tracking performance. They designed an energy function that encodes various personal, social and environmental factors to determine the next velocity of pedestrians. They used a simplex algorithm to learn optimal parameters of the energy function based on training data extracted from videos. Similarly, Scovanner and Tappen [11] proposed a new energy function with 17 parameters to predict movements of pedestrians. They utilized a modified Newton’s method to optimize model parameters so that the simulated movement match tracks in the videos. In [30] a density-based evaluation method was proposed to evaluate the quality of simulation output. Differential evolution (DE) was adopted to calibrate parameters of two modified SFMs so that the simulated crowd density matched real crowd density calculated based on video data. In [14], an enhanced DE algorithm that included sensitivity analysis and a local search procedure was proposed to improve the efficiency of parameter calibration. The main difference between the work described in this paper and the existing ones described above, is that these existing works focus on using data to calibrate parameters of models that are designed by humans; whereas in this paper we aim to automatically distill behavioral rules from data. The second group of methods focuses on automating the generation of crowd behaviors from data. The core idea is to derive typical

370

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

motion patterns from data, and then to use these derived patterns to drive agents’ movement. The existing approaches differ in the definition of typical and the way to determine which pattern is to be used during the simulation to drive agent’s movement. For example, Lee et al. [15] proposed to use a hierarchical model to generate behaviors of agents. They used a finite state machine to control the transitions between high-level activities such as joining or leaving a group, while using state-action pairs to control the low-level motions of agents. Each state-action pair determines the future velocities for agents under particular situations. Their work focused on learning state-action pairs from video to control the low-level motions of agents. Lerner et al. [17] proposed a data-driven crowd modeling method which requires extracting a database of patterns from video in advance. Each pattern is a short trajectory segment. During the simulation, when the trajectory segment of an agent is exhausted or there is a significant change in the current situation, the agent will query a similar pattern to extend its trajectory. Zhao et al. [19] proposed a data-driven crowd modeling approach, which also uses patterns extracted from videos to update the velocities of agents. Musse et al. [16] proposed a data-driven crowd modeling approach for simulating human behaviors in normal life. They developed an automatic method to extract trajectories from videos of non-dense scenarios. Then the extracted trajectories are used to learn the desired velocities of agents in different regions. Eunjung et al. [18] proposed a data-driven crowd modeling approach which can construct a large-scale crowd that looks similar to the input video, in terms of crowd formation and moving styles. In their model, a database of formations and trajectory segments are firstly extracted from videos. Then, each agent iteratively selects a trajectory segment from the database to extend its moving trajectory so that it dose not collide with other agents and the formation constraint is also satisfied. In the computer vision community, there has been an increasing amount of work on learning global motion patterns of crowds from video data [21–25,31]. For example, Zhou et al. [31] proposed to use Linear Dynamic System (LDS) to model the movement curves of pedestrians in the video data. They utilized an Expectation Maximization algorithm to learn parameter settings of the LDS so that the movement curves offered by the LDS can match those in the video. The data-driven approaches above hope to design crowd behavior by making many observations, and that through some statistical analysis of many observations, general behaviors can be generated. This approach can work well, but it typically isn’t general enough to apply in different scenarios, or different environments. Also, as there is no underlying model (or theory) nothing is really understood about why the crowd moves as it does. The approach we adopt in this paper also uses data, but we try to formulate a general model (or theory) from the data. That is we aim to learn a generic behavioral rule, rather than creating a database of examples for different scenarios. The hope is that with this

more Popperian approach, the behavioral rule distilled is a general description of behavior and as such can be applied to different new scenarios. 3. Data-driven agent-based crowd modeling framework This section proposes a generic data-driven agent-based framework to generate realistic crowd behaviors. First, the general structure of the proposed framework is described. Then, the mechanism of automatically learning behavior rule from video data is described. 3.1. Data-driven agent-based crowd modeling framework The procedures of the proposed data-driven modeling framework are illustrated in Fig. 1. First of all, a navigation rule is distilled from crowd data by using GP. Then, the distilled rule is utilized in a dual-layer simulation model to determine the desired moving direction of an agent. In the dual-layer simulation model, the top layer adopts a “Sense-Think-Act” paradigm to determine the next moving directions of pedestrians. Environmental features, such as the distance to the nearest obstacle, are used as inputs to the navigation rule. The navigation rule determines a moving direction using a mapping function  that maps the input environmental features to a moving direction. Once the moving direction is calculated by  , it is combined with a desired speed to determine a desired next way point for the agent. This desired next way point is then used as the input to the bottom layer which modifies the velocity according to the SFM [8] for collision avoidance. 3.2. Learning  from video In this section, the problem of finding  from real world data is formulated as a symbolic regression problem. Specifically, we consider the data extracted from videos as individual trajectories, each of which is represented by a vector of points. Our objective is to identify a suitable  based on these extracted trajectories. First, the extracted trajectory data are transformed into a series of input–output pairs: (I1 , 1 ), . . ., (IN , N )

(1)

where N is the total number of input–output pairs; Ii represents the input feature values of the ith pair, and  i is the corresponding next moving direction. In our approach, the modelers simply have to specify the features that are important to the agents (individuals) when determining their next moving directions. That is, the modelers, must define what terms should appear in the mapping function  . In this paper, we consider three input features for function  . As illustrated in Fig. 2, the first feature is the distance to the goal (i.e., |AC|), denoted as dg . The second feature is the distance to

Fig. 1. The proposed data-driven modeling framework.

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

Fig. 2. The three input features under consideration.

the nearest obstacle (i.e., |AB|), denoted as do . The third feature is  and AB  (i.e., ∠CAB), denoted as  go . It should the angle between AC be noted that there are some trade-offs when defining the input features. On one hand, if important features are not considered, the algorithm will be not able to find a satisfying solution. On the other hand, considering unnecessary input features will increase the search space and may reduce the search efficiency of the algorithm. In general, the approach should try to find the minimum set of input features that has the most significant impact on individuals’ decisions. Once a particular  is determined, it is possible, and advisable, to perform sensitivity analysis for all parameters, to assess the significance in the generated model. To generate the input–output pairs from the trajectory data, the first step is to divide the simulation region into a discrete grid with each grid cell having a size of  × . The value of  determines the number of input–output pairs to be generated. Generally, if the extracted trajectories include frequent direction changes, then a smaller  should be adopted so that more input–output pairs are used to describe the movement curves of pedestrians. In the second step, each trajectory is represented by a small number of segmentation points. The segmentation points are extracted from the original data points in a given trajectory, so that the length of each segment between √ any two consecutive points is equal to or slightly larger than 2 (i.e., the diagonal length of each grid cell). By doing so, we ensure that each grid cell contains at most one segmentation point for each trajectory. The following vector of segmentation point can then be used to represent the ith trajectory: Ti = s1i , s2i , . . ., ski  i

(2)

where sji is the jth segmentation point of the ith trajectory and ki is the total number of segmentation points of the ith trajectory. In the third step, all trajectories are classified according to their destination regions (DRs) so that trajectories in the same class have the same DR. Based on the classification results, we can estimate the goal of each trajectory, which is useful for calculating dg and  go . In this study, our previously proposed K-means method [20] is used to accomplish this task. In the last step, a set of input–output pairs for each DR are generated based on the trajectories of individuals whose destinations are associated to that DR. For simplicity, we denote the set of segmentation points of the trajectories associated to the ith DR as Ci , the center point of grid cell j as pj , and the position of the ith DR as gi . To generate a set of input–output pairs based on Ci , the grid cells are checked one-by-one. Without loss of generality, we suppose there are  segmentation points (not including the last segmentation points of the trajectories because these points do not have their next segmentation points) ∈ Ci falling into the jth grid cell, and the subsequent segmentation points of these  points are s1 , s2 , . . ., and s respectively. Then, we can estimate the next

371

way point the trajectories lead to based on the average of these segmentation points (denoted as s¯ ). If  > 0, a new input–output pair is generated for the grid cell, where the input feature values are calculated based on pj and the output is the angle from vector AC to vector AD, where A = pj , C = gi , and D = s¯ . In this way, the input–output pairs generated based on all Ci form the input–output pairs for training. Based on the N input–output pairs, we hope to find a generic mathematical formula to describe the relationships between the input features and the output. The mathematical formula is composed of the input features (known as terminals in GP) and functions such as +, −, and tan. Similar to the definition of input features, the functions to be used are chosen by the modeler. Generally, if important functions are missing, the algorithm will be unable to find satisfying results. As the optimal  is unknown, the necessary functions to construct  are unknown. The major factors to affect the choice are the types of the terminals and outputs. In this paper, the terminals are distance and angle values, while the output is an angle value.  aims to map distance and angle values to an angle value. Thus, we adopt functions that are commonly used to calculate distance (e.g., +, −, ×, /) and angle values (e.g., tan, and atan) as building blocks, i.e., the final function set is {+, −, *, /, tan, atan}. With the above descriptions, the rule identification problem can be transformed into the following optimization problem: Given a set of input–output pairs, a feature set {dg , do ,  go } and a function set {+, −, *, /, tan, atan} as building blocks, the task is to find the optimal mapping function  * that minimizes the prediction error:



∗

 = argmin 

N 2 ( (Ii ) − i ) i=1

N

(3)

where N is the total number of input–output pairs,  (Ii ) is the output of  for the ith data, and  i is the output obtained based on the real data. The above problem can be viewed as a symbolic regression problem which is commonly solved by using GP and its variants [32]. In Section 4, we will describe how to apply a recent published GP variant named SL-GEP [29] to solve the above symbolic regression problem. 4. Distil behavioral rules by using SL-GEP The SL-GEP [29] is a new GP variant for automatic function generation. Compared with the traditional GP and GEP, the SLGEP is easy to implement and contains fewer control parameters. The SL-GEP has been shown to be more efficient than the GP and the GEP on a range of problems, including the symbolic regression problems. Hence, in this paper, the SL-GEP is adopted to distill a suitable mapping function  to fit the generated input–output pairs. In this section, the chromosome representation of SL-GEP is described at first by considering the problem defined in Section 3. Then, a brief introduction of the algorithm framework of SL-GEP is presented. The readers are referred to [29] for the details of the SL-GEP. 4.1. Chromosome representation Mathematical formula and logical rules are commonly expressed using two types of symbols: function and terminal. A function (e.g., +, tan) can map its inputs to an output value, while a terminal (e.g., variable dg and constant ) represents a value. In the SL-GEP, a fixed-length string with two parts is used to represent computer programs. The first part is the main program that provides the final output. The second parts are evolvable subfunctions which can be used in the main program. Both the main program and subfunctions are represented using a gene expression technique proposed by Ferreira [33]. For each gene expression string, there

372

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

4.2.1. Step 1 – initialization This step randomly generates NP chromosomes as the initial population. Each chromosome in the population (named a target vector) is represented by a vector of symbols: Xi = [xi,1 , xi,2 , . . ., xi,n ]

(8)

where n is the number of symbols in each chromosome. To set the value of xi,j , a feasible type is firstly randomly selected among {function, sub-function, terminal, input argument}. Then, xi,j is set to a random value of the selected feasible type. For example, suppose the selected feasible type is terminal, then xi,j can be set to dg , do , or d . 4.2.2. Step 2 – mutation The second step performs a mutation operation on each target vector to generate a mutant vector. The SL-GEP adopts the “DE/current-to-best/1” mutation operation as can be expressed by Yi = Xi + F · (Xbest − Xi ) + ˇ · (Xr1 − Xr2 ).

Fig. 3. A typical example chromosome and its decoded expression trees. (A) An example behavioral rule represented in the chromosome. (B) The decoded expression trees of the main program and the sub-function. (C) The final decoded expression trees of the whole chromosome.

is an equivalent expression tree. The gene expression string can be obtained by a breadth first traversal of the expression tree. Fig. 3 shows an example behavioral rule represented in the chromosome of SL-GEP, which can be expressed as: [G, atan, G, G, do , dg , dg , do , do , ∗, tan, ∗, b, a, b, b]

(4)

where {+, −, *, /, tan, atan} is the function set, G is a sub-function, {do , dg ,  go } is the terminal set, and {a, b} is the set of input arguments. For this chromosome, the main program part is {G, atan, G, G, do , dg , dg , do , do }. As shown in Fig. 3, by using a breadth first traversal scheme, the main program can be decoded as G(atan(G(dg , do )), G(do , dg ))

(5)

The sub-function part is { *, tan, *, b, a, b, b}. By using the same traversal scheme, the sub-function can be decoded as G(a, b) = tan(b) ∗ (a ∗ b)

 = tan(tan(dg ) ∗ (do ∗ dg )) ∗ (atan(tan(do ) ∗ (dg ∗ do )) ∗ (tan(dg ) ∗ (do ∗ dg )))

where Xbest is the best-so-far solution in the population, and Xr1 and Xr2 are two random individuals in the population. Since Xi , Xr1 , Xr2 , and Xbest are vectors of discrete symbols, the numerical operators in the traditional “DE/current-to-best/1” mutation operation such as · and + are not applicable. In the SL-GEP, these numerical operators are redefined so as to evolve solutions in a discrete search space. By using the newly defined operators, each element of Yi is assigned by mutating the corresponding element in Xi . The mutation probability ϕ is calculated by the following process. Firstly, the values of F and ˇ are randomly set by: F = rand(0, 1)

(10)

ˇ = rand(0, 1)

(11)

where rand(a, b) returns a random value uniformly distributed within [a, b]. Then, for each element (xi,j ) of the target vector, a mutation probability is calculated by ϕ = 1 − (1 − F · where

(7)

(xbest,j , xi,j )) ∗ (1 − ˇ ·

(xr1,j , xr2,j ))

(12)

(a, b) is defined as

 (a, b) =

(6)

Inserting (6) into (5), we then obtain the final expression of  as:

(9)

1,

if a = / b

0,

otherwise

(13)

The probability that yi,j has the same value as xi,j is 1 − ϕ, and the probability that yi,j is set to a new value (called mutant value) is ϕ. If yi,j ∈ sub-function, it will be set to a random value as done in the initialization step. Otherwise, the type of yi,j is set first by considering the frequencies of feasible types appearing in the population. Feasible types that appear more often in the population are more likely to be selected. Once the feasible type of yi,j is determined, yi,j is then set to a random value of the selected feasible type.

4.2. Algorithm framework With the above chromosome representation, the SL-GEP firstly generates a population of random chromosomes. Each chromosome represents a candidate solution (i.e., a candidate  ). Then, genetic operators such as mutation and crossover are performed to iteratively evolve the population until meeting the termination condition. Finally, the best chromosome in the population is decoded to obtain the best distilled rule, which can then be applied in the crowd simulations. Specifically, the procedure of SL-GEP contains the following four steps.

4.2.3. Step 3 – crossover In the crossover step, each target vector Xi is crossover with its mutant vector Yi to generate a trial vector Ui :



ui,j =

yi,j ,

if rand(0, 1) < CR

xi,j ,

otherwise

or

j=k

(14)

where CR with CR = rand(0, 1) is the crossover rate, k is a random integer between 1 and n, ui,j , yi,j and xi,j are the jth variables of Ui , Yi and Xi respectively.

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

373

Fig. 4. Scenarios for training and testing and the feature set under considerations. (A) Scenario for training. (B) Scenario for testing.

4.2.4. Step 4 – selection The selection operation selects the fitter solution between each pair of the target and trial vector to form a new population:



Xi =

Ui ,

if f (Ui ) ≤ f (Xi )

Xi ,

otherwise

(15)

where f(X) is the fitness evaluation function which calculates the prediction error of X (see (3)). Steps 2–4 are repeated until the maximum number of generations is reached. 5. Simulation studies 5.1. Application of our method In this paper, we focus on modeling the walking behavior of human in normal cases where people walk casually. More specifically, we aim to construct a human behavior model to generate crowd movements in an area with one or multiple SRs and DRs. The pedestrians move casually from SRs to DRs by following certain paths of motion, while at the same time avoiding collisions with nearby obstacles and other pedestrians. Our objective is to automatically learn generic  from data (e.g., video feeds) which can be applied  in other scenarios where people also walk casually. To achieve this goal, we first apply the proposed method to learn  from a video of walking pedestrians in Switzerland (training dataset). Fig. 4(A) shows the structure of the walking region.1 The training dataset contains 367 trajectories of pedestrians in the Eidgenössische Technische Hochschule (ETH) Zürich. The trajectory points are annotated at 2.5 fps. The goal is to learn the following rule  based on the trajectory data:  =  (dg , do , go )

(16)

where  describes the desired moving direction of an agent and dg , do , and  go are the three terminals that used to construct  (see Fig. 2). Then, we test whether the extracted  is general enough to be applied to different scenarios. To do this we use a testing dataset from New York Grand Central Station, as shown in Fig. 4(B), which is more complex and involves more entrances and exits.2 We do not generate a new rule from the testing dataset, instead we test

1 The related video and dataset of the training scenario can be downloaded from: http://www.vision.ee.ethz.ch/datasets/index.en.html. 2 The video and extracted trajectories of the testing scenario can be downloaded from: http://www.ee.cuhk.edu.hk/∼xgwang/grandcentral.html.

the same rule extracted from the training dataset and see if it can generate similar crowd behavior as the New York video data. The experiment is served as a proof-of-concept of the generality of the proposed method.3 We use the SL-GEP as described in Section 4 to learn  . In SLGEP, the population size is set to 100 and the maximum number of generations is set to 20,000. It should be noted that the maximum number of generations is set based on two factors: the convergence state of the algorithm and the maximal time budget. Since the fitness evaluation is not expensive (0.2 ms), we used a relatively large maximum number of generations (i.e., 20,000). With this number of generations, we observe that the best-so-far fitness value generally becomes stable based on our trial runs. Other parameters are set the same as in [29]. We perform SL-GEP for 30 independent runs with different random seeds on a computer with Intel(R) Core(TM) i7-7500U [email protected] 2.90 GHz, and 8G memory. 5.2. Algorithm analysis First, we analyze the computational cost of the proposed method. In general, the total computational cost of the proposed method is determined by two major factors: (1) the total number of fitness evaluations (i.e., maximum generation * population size), and (2) the computational cost of each fitness evaluation. With the experiment settings given in Section 5.1 and the results of 30 SLGEP runs, the average computational time of each SL-GEP run is 556 s. It should be noted that the training process is an offline procedure and the solution provided by SL-GEP (i.e.,  ) is used in the simulation of new scenarios. During the simulation, the time taken by each agent to evaluate  is insignificant. Thus, the computation time of the proposed method is acceptable. Next, we investigate the impacts of function sets. To achieve this, we investigate the performance of the proposed method with different function sets. Algorithm A1 uses function set {+, −, *, /, tan, atan}, A2 uses three basic functions (i.e., {+, −, * }), A3 uses four basic functions (i.e., {+, −, *, /}), A4 uses eight basic functions (i.e., {+, −, *, /, tan, atan, exp, log}), and A5 uses ten basic functions (i.e., {+, −, *, /, tan, atan, sin, cos, exp, log}). Fig. 5 shows the evolution curves of the best fitness verse the number of generations. It can be observed that the algorithms with smaller function set (i.e., A2 and A3) perform much worse than A1. This results demonstrate that tan and atan are useful to construct the final  , because the behavior rules in this scenario are basic geometry. This result also indicates that if useful

3 The source code of the proposed method, training and testing data sets, simulation results, and related videos can be downloaded from: https://www.dropbox. com/sh/z8eos0qt63pqsv5/AAAgK4SeXdH22EBszA25fiB a?dl=0.

374

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

the 30 runs is expressed in Eq. (17).  = atan(((G0 (G0 (do , go ), dg )/(atan(go ) − dg )) − (G0 (dg , do )/ tan(go ))))/(G0 (do , (go − atan(dg ))) + do )

Fig. 5. The best fitness values verse the number of generations.

functions are not considered, the performance of the algorithm will reduce significantly. On the other hand, A4 and A5 perform similarly to A1. This may be because that new functions such as exp and log are not very useful to construct the final  . Thus, increasing the size of function set may not always be helpful to improve the search efficiency. It is worth noting that if problem-specific knowledge is available, then we can manually define user-defined functions as elements of the function set. These user-defined functions may be useful to improve the search efficiency.

5.3. Training results In the following parts, we use the best solution found by A1 for analysis and testing. Specifically, the best solution found by A1 in

(17)

where G0 (t1 , t2 ) = atan((t2 * t2 )) is a subfunction of  . To demonstrate the effectiveness of the rule, we plot the velocity fields generated by Eq. (17) which has the best fitness value. We assign each DR with a velocity field so as to capture the movement curves of pedestrians that move toward the DR. Specifically, the velocity field assigned to the ith DR is denoted as Fi = {p1 , vi1 , p2 , vi2 , . . ., < pC , viC }, where p1 , . . ., pC are C representative positions and vij is the velocity at pj for the ith DR. The representative positions are obtained by evenly dividing the entire region into discrete grids, with each grid having a size of 1m × 1m. The centers of grid cells that contain no obstacle are used as representative positions. To calculate vij , we first use  to obtain a desired moving direction ji at pj . Note that the positions of the goals and obstacles are known in advance, we can calculate the three feature values (i.e., dg , do , and  go ) based on the position of pj . Thus, we can obtain the desired moving directions of all representative points by using  . Then, the direction of vij is set equal to ji and the length of vij is set equal to a preferred speed v0 (e.g., ||vij || = v0 = 1 m/s). In this way, we can obtain velocity fields of all DRs. By plotting the velocity fields of all DRs in a 2D space, we can visually evaluate whether the velocity fields describe plausible motion. Fig. 6 illustrates two example velocity fields generated by Eq. (17) and those generated based on the trajectory data. It can be observed that the velocity fields generated based on the trajectory data are incomplete because there aren’t trajectories covering the entire region. However, Eq. (17) is able to generate velocities of the entire region. Instead of pointing to the destinations directly, the velocities generated by Eq. (17) form realistic smooth movement curves which can

Fig. 6. Examples of velocity-fields generated by the distilled rule in the training scenario. (A) Velocity-field of DR1 generated based on the trajectory data. (B) Velocity-field of DR1 generated by Eq. (17). (C) Combination of velocity-fields in (A) and (B). (D) Velocity-field of DR3 generated based on the trajectory data. (E) Velocity-field of DR3 generated by Eq. (17). (F) Combination of velocity-fields in (D) and (E).

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

375

Fig. 7. The outputs of varying inputs generated by Eq. (17). (A) dg is fixed as 10 m. (B) do is fixed as 3 m. (C)  go is fixed as /6.

guide agents to avoid the obstacles in the environment and reach their destinations with relatively short distances. Furthermore, the comparison results in Fig. 6(C) and (F) demonstrate that the velocity fields generated by Eq. (17) can resemble those generated based on the trajectory data. A distinct advantage of the proposed method over other datadriven modeling approaches is that the proposed method is capable of providing explicit formulas to describe the relationships between input features and the outputs. By analyzing the explicit formulas, we can gain more insights of the physical system. For example, to investigate the relationships between features and the output, we calculate the outputs of Eq. (17) by fixing one feature and varying the other two features. Fig. 7 shows the outputs of Eq. (17) with varying inputs. In Fig. 7(A),  go has larger influence on the output when do is small (e.g., do < 4). If do is large enough (e.g., do > 4), the output becomes relatively stable and small. This indicates that agents will directly set the moving direction toward the goal if there is no immediate obstacle nearby, which is consistent with human’s steering strategy in practice. The contours in Fig. 7(B) and (C) indicate that when dg is large enough (e.g., dg > 3), the output becomes mainly dependent on  go and do . This means that if the goal is far away, how agent moves is mainly affected by the distance to the near-by obstacle and angle between the obstacle and the goal. To find out which feature has the largest influence on the crowd dynamics, we further perform a sensitivity analysis on Eq. (17) by using Morris’s method [34]. The sensitivity analysis results indicate that  go (i.e., the angle between the nearest obstacle and the destination) has the largest influence. This implies that our method can also be used to identify what are the most critical decision factors that affect the agent’s decisions.

method [20]. If a grid cell does not contain sufficient segmentation points of trajectories, then the estimated objective velocity at the representative point is considered not accurate enough and is considered invalid. In this study, a grid cell contains less than 15 segmentation points is considered invalid. For the ith DR, the error at representative point pj is calculated by: i = I(ωji ) · |ωji − ji |

where ωji is the objective moving direction and ji is the moving direction given by  , I(ωji ) returns one if ωji is valid and zero otherwise. The value of  error is calculated by:

M C error =

i=1

To test the generality of the extracted rule, we apply it to the testing dataset. The testing dataset contains 40,000 trajectories which are extracted from a 33-minute video. The trajectory points are annotated at 24 fps. First, we plot the velocity fields generated by Eq. (17) in this new scenario. Fig. 8 demonstrates two example velocity fields generated by Eq. (17). It can be observed that the distilled rule is able to provide complete velocity fields over the entire region. The velocities generated by the distilled rule are not directly pointing to the destinations. Instead, they form “C” shape curves that avoid the obstacles in the central region of the environment. According to the shapes of the velocity fields, the velocity fields generated by the distilled rule look very similar to the objective results. To quantitatively measure the performance of the distilled rule in the new scenario, two performance metrics are adopted. The first metric (denoted as  error ) measures the error of velocity fields. To calculate  error , we first estimate the objective velocity field of each DR based on the real trajectory data by using the pattern learning

j=1

M C i=1

i j

(19)

I(ωji ) j=1

where M is the total number of DRs, and C is the total number of representative points. The second metric (denoted as error ) measures the error of simulated crowd density distribution. To calculate error , we first calculate the real density distribution based on all trajectories extracted from the video (i.e., the accumulated density during the whole observed period). Then the simulation model with Eq. (17) is run using the same number of time steps as in the video data. Other crowd statistics such as the arrival rate of pedestrians in each ER are set the same as in video data. The C representative positions are used to estimate the crowd density distribution. The accumulated local density of the representative point pi is computed by [27]: 1  exp R2 m

5.4. Applying  to new scenario

(18)

j

pi =

k=1

 −

d(pi , k) R2

2

 (20)

where m is the number of points in all trajectories under consideration, d(pi , k) returns the distance from the kth trajectory point to pi , and R (here R = 2) is a scaling parameter. In this way, a trajectory point that is closer to the representative point will have a larger contribution to its local density value. Note that when calculating the crowd density from the video, m is set to be the number of points in all extracted trajectories. Otherwise, m is set to be the number of points in all trajectories generated by the simulation (the simulation time step is set equal to the annotated rate). Based on Eq. (20), we can calculate the error of crowd density at pj by:  j = | j − j0 |

(21)

where j is the simulated crowd density at the representative point pj , and j0 is the crowd density calculated based on video data. We perform the simulation model for 30 independent runs with different random seeds and use the average crowd density of the 30

376

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

Fig. 8. Examples of velocity-fields generated by the distilled rule in the testing scenario. (A) Velocity-field of DR1 generated based on trajectory data. (B) Velocity-field of DR1 generated by Eq. (17). (C) Combination of velocity-fields in (A) and (B). (D) Velocity-field of DR6 generated based on trajectory data. (E) Velocity-field of DR6 generated by Eq. (17). (F) Combination of velocity-fields in (D) and (E).

simulation results to estimate j . Then, the value of error is calculated by

C error =

j=1

j

compared to the other two approaches. The average values of i j

 j

(22)

C

To demonstrate the effectiveness of the distilled  , we compare it with two other navigation strategies. The first strategy simply set the desired moving directions of agents to be those pointing to the agents’ destinations. The second strategy is the shortest path strategy, which is commonly used to generate navigation paths of pedestrians in crowd simulation [35]. In the second strategy, the shortest path from each representative point to each DR is calculated in advance. Then, agents near the representative point will use the shortest path associated to the representative point to determine their desired moving directions. To calculate the shortest path from a representative point to a DR, we first constructed a connected graph where the nodes of the graph are the DR and all representative points. Each node can only connect to nearby nodes within a certain range (r). The nodes (i.e., the representative points) are used as candidate way-points of the shortest path. Then, based on this connected graph, we can obtain the shortest path from the representative point to the DR by using the Dijkstra algorithm. In this study, r is set to 5 m, which can generate satisfying paths that can avoid the obstacle in the center region. To facilitate description, we denote the simulation models that use the first strategy, the second strategy, and the distilled rule  as M1, M2, and M3, respectively. Fig. 9(A) shows the distribution of i , where the x-axis represents the value range of i (with j

ranges and less counts on high error ranges for both i and  j

j

an interval length of 0.1), and the y-axis represents the count of i that falls into the value range. Fig. 9(B) shows the correspondj

ing results of  j where the interval length of  j is set to 400. It can be observed that M3 generally has more counts on low error

for M1, M2, and M3 are 0.274, 0.288, and 0.257 respectively. Meanwhile, the average values of  j for the three models are 1483, 1473, and 1318 respectively. Since the testing scenario is relatively regular and does not contain many obstacles, real pedestrians are likely to move directly toward the destination in most regions, thus M1 performs better than M2 in terms of i . M2 performs slightly j

better than M1 in terms of  j . This may be because that M2 can provide more realistic velocities in dense regions (e.g., those near the center region) than M1. For both metrics, our proposed M3 always performs the best among the three methods. Furthermore, we perform a Wilcoxon signed rank test to investigate the significant difference between results of M3 and those of the other two methods. For i , the P-values are 5.02E−13 and 2.2E−16 for M1 j

and M2 respectively. For  j , the P-values are 0.0067 and 0.0075 for M1 and M2 respectively. The Wilcoxon signed rank test results indicate that M3 performs significantly better than M1 and M2 in terms of i and  j at 95% confidence level. These results demonj

strate that using the proposed approach, the behavior rule learned from one scenario can be effectively applied to describe the crowd dynamics in a totally different scenario. 6. Conclusions and discussions This paper has proposed an automatic methodology to extract generic behavioral rules from video for agent-based crowd modeling. Most importantly we have shown that this automatic approach to modeling can extract rules from a dataset in Switzerland and then be applied to generate realistic behavior for a crowd in the United States. We believe that this approach may be an effective

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

377

Fig. 9. Distributions of i and  j . (A) Distributions of i . (B) Distributions of  j . j

method for systematic and automated creation of agent-based models in general. The proposed methodology is capable of learning a navigation rule that in turn defines velocity fields that can be used to drive pedestrians toward their destinations. In our approach we, as the modelers, simply have to tell the algorithm what features are important to the agents (individuals) when determining their desired moving direction. That is, we must define what terms should appear in the function  . It should be noted that if important features and functions are not considered, the algorithm will not be able to find a satisfying solution. On the other hand, if unnecessary input features (or functions) are considered, the search efficiency of the algorithm may reduce. In our simulation study, the behavior model is trained from the video data where people walk casually. If the walking behaviors of people (e.g., those in rush hours or in emergent situations) are much different from those in training data, the trained model may not be directly applicable. In these cases, we need to retrain a new model for the new scenario. The retraining efforts include generating new training data, defining new terminal set and function set for the specific scenario, and re-running the proposed method based on the new training data, terminal and function sets. In addition, if there are significant changes in the volume and behavior of the crowd in the training video, the proposed methodology may not be able to generate satisfying results. It would be an interesting research topic to design an automatic method that can distill different behavior rules from videos that have different behavior patterns in terms of density and speed. In future, we would try to do some research on this direction and study the impacts of the function and terminal set for different scenarios. The proposed methodology is not limited only to generating navigation rules. It can also be used to generate other behavioral rules of crowd models or other behavior models (e.g., to discover the flocking rules of birds and the behavioral rules of schooling fish). The key idea is to generate a set of input–output pairs from the real trajectory data and then utilize GP-based methodology to search for optimal explicit formulas that describe the relationships between inputs and outputs. However, an important question for the modeler is how to determine the key input features (or the variables) and functions that play important roles in the final derived behavior rule (i.e.,  ). The proposed method has a number of potential practical applications. For example, we can apply the proposed method to predict the trajectories of pedestrians, which is useful for crowd tracking. Besides, by comparing the real trajectories with predicted trajectories, the proposed method can also be used to identify

j

abnormal behaviors of pedestrians. Applying the proposed method to such real world practical applications is a promising future research work. In addition, considering the rule complexity as another optimization objective and utilize multi-objective optimization techniques to find better and more concise behavioral rules can be another interesting research direction.

Acknowledgments The research reported in this paper is supported by National Natural Science Foundation of China (Grant Nos. 61602181 and 61502370), Tier 1 Academic Research Fund (AcRF) under project Number RG23/14, Natural Science Basic Research Plan in Shannxi Province of China (Program No. 2016JQ6003) and China 111 Project (Program No. B16037). M. Lees would like to acknowledge the Russian Science Foundation project #14-21-00137 “Supercomputer modeling of critical phenomena in complex social systems”.

References [1] W. Shao, D. Terzopoulos, Autonomous pedestrians, in: Proceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, ACM, 2005, pp. 19–28. [2] D. Thalmann, S.R. Musse, Crowd Simulation, Springer, London, 2007. [3] N. Gaud, F. Gechter, S. Galland, A. Koukam, Holonic multiagent multilevel simulation: application to real-time pedestrian simulation in urban environment, in: Proceedings of the 2007 International Joint Conference on Artificial Intelligence, 2007, pp. 1275–1280. [4] L. Luo, S. Zhou, W. Cai, M.Y.H. Low, F. Tian, Y. Wang, X. Xiao, D. Chen, Agent-based human behavior modeling for crowd simulation, Comput. Anim. Virtual Worlds 19 (3–4) (2008) 271–281. [5] J. Tsai, N. Fridman, E. Bowring, M. Brown, S. Epstein, G. Kaminka, S. Marsella, A. Ogden, I. Rika, A. Sheel, et al., Escapes: evacuation simulation with children, authorities, parents, emotions, and social comparison, in: Proceedings of the 2011 International Conference on Autonomous Agents and Multiagent Systems, IFAAMS, 2011, pp. 457–464. [6] S. Zhou, D. Chen, W. Cai, L. Luo, M.Y.H. Low, F. Tian, V.S.-H. Tay, D.W.S. Ong, B.D. Hamilton, Crowd modeling and simulation technologies, ACM Trans. Model. Comput. Simul. 20 (4) (2010) 1–35. [7] D. Helbing, P. Molnar, Social force model for pedestrian dynamics, Phys. Rev. E 51 (5) (1995) 4282. [8] D. Helbing, I. Farkas, T. Vicsek, Simulating dynamical features of escape panic, Nature 407 (6803) (2000) 487–490. [9] K. Yamaguchi, A.C. Berg, L.E. Ortiz, T.L. Berg, Who are you with and where are you going? in: Proceedings of the 2011 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2011, pp. 1345–1352. [10] S. Pellegrini, A. Ess, K. Schindler, L. Van Gool, You’ll never walk alone: modeling social behavior for multi-target tracking, in: Proceedings of the 2009 IEEE International Conference on Computer Vision, IEEE, 2009, pp. 261–268. [11] P. Scovanner, M.F. Tappen, Learning pedestrian dynamics from the real world, in: Proceedings of the 2009 IEEE International Conference on Computer Vision, vol. 9, IEEE, 2009, pp. 381–388.

378

J. Zhong et al. / Applied Soft Computing 56 (2017) 368–378

[12] A. Johansson, D. Helbing, P.K. Shukla, Specification of the social force pedestrian model by evolutionary adjustment to video tracking data Adv. Complex Syst. 10 (Suppl 02) (2007) 271–288. [13] D. Wolinski, S.J. Guy, A.-H. Olivier, M. Lin, D. Manocha, J. Pettré, Parameter estimation and comparative evaluation of crowd simulations, in: Computer Graphics Forum, vol. 33, Wiley Online Library, 2014, pp. 303–312. [14] J. Zhong, N. Hu, W. Cai, M. Lees, L. Luo, Density-based evolutionary framework for crowd model calibration, J. Comput. Sci. 6 (2015) 11–22. [15] K.H. Lee, M.G. Choi, Q. Hong, J. Lee, Group behavior from video: a data-driven approach to crowd simulation, in: Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, ACM, 2007, pp. 109–118. [16] S.R. Musse, C.R. Jung, J. Jacques, A. Braun, Using computer vision to simulate the motion of virtual agents, Comput. Anim. Virtual Worlds 18 (2) (2007) 83–93. [17] A. Lerner, Y. Chrysanthou, D. Lischinski, Crowds by example, in: Computer Graphics Forum, vol. 26, Wiley Online Library, 2007, pp. 655–664. [18] E. Ju, M.G. Choi, M. Park, J. Lee, K.H. Lee, S. Takahashi, Morphable Crowds, vol. 29, ACM, 2010, pp. 1–10. [19] M. Zhao, S.J. Turner, W. Cai, A data-driven crowd simulation model based on clustering and classification, in: Proceedings of the 2012 IEEE/ACM International Symposium on Distributed Simulation and Real Time Applications, IEEE, 2013, pp. 125–134. [20] J. Zhong, W. Cai, L. Luo, Y. Haiyan, Learning behavior patterns from video: a data-driven framework for agent-based crowd modeling, in: Proceedings of the 2015 International Conference on Autonomous Agents and Multi-agent Systems, IFAAMS, 2015, pp. 801–809. [21] S. Ali, M. Shah, A Lagrangian particle dynamics approach for crowd flow segmentation and stability analysis, in: Proceedings of the 2007 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2007, pp. 1–6. [22] W. Ge, R. Collins, R. Ruback, Vision-based analysis of small groups in pedestrian crowds, IEEE Trans. Pattern Anal. Mach. Intell. 34 (5) (2012) 1003–1016. [23] D. Lin, E. Grimson, J. Fisher, Learning visual flows: a lie algebraic approach, in: Proceedings of the 2009 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2009, pp. 747–754.

[24] D. Lin, E. Grimson, J. Fisher, Modeling and estimating persistent motion with geometric flows, in: Proceedings of the 2010 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2010, pp. 1–8. [25] R. Mehran, B.E. Moore, M. Shah, A streakline representation of flow in crowded scenes, in: Proceedings of the 2010 European Conference on Computer Vision, Springer, 2010, pp. 439–452. [26] J. Zhong, L. Luo, W. Cai, M. Lees, Automatic rule identification for agent-based crowd models through gene expression programming, in: Proceedings of the 2014 International Conference on Autonomous Agents and Multi-agent Systems, International Foundation for Autonomous Agents and Multiagent Systems, 2014, pp. 1125–1132. [27] D. Helbing, A. Johansson, H.Z. Al-Abideen, Dynamics of crowd disasters: an empirical study, Phys. Rev. E 75 (4) (2007) 046109. [28] X. Jin, J. Xu, C.C. Wang, S. Huang, J. Zhang, Interactive control of large-crowd navigation in virtual environments using vector fields, IEEE Comput. Graph. Appl. (6) (2008) 37–46. [29] J. Zhong, Y.-S. Ong, W. Cai, Self-learning gene expression programming, IEEE Trans. Evol. Comput. 20 (1) (2015) 65–80. [30] J. Zhong, W. Cai, Differential evolution with sensitivity analysis and the Powell’s method for crowd model calibration, J. Comput. Sci. 9 (2015) 26–32. [31] B. Zhou, X. Wang, X. Tang, Understanding collective crowd behaviors: learning a mixture model of dynamic pedestrian-agents, in: Proceedings of the 2012 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2012, pp. 2871–2878. [32] J. McDermott, D.R. White, S. Luke, L. Manzoni, M. Castelli, L. Vanneschi, W. Jaskowski, K. Krawiec, R. Harper, K. De Jong, et al., Genetic programming needs better benchmarks, in: Proceedings of the 2012 Genetic and Evolutionary Computation Conference, ACM, 2012, pp. 791–798. [33] C. Ferreira, Gene Expression Programming: Mathematical Modeling by an Artificial Intelligence, vol. 21, Springer, 2006. [34] M.D. Morris, Factorial sampling plans for preliminary computational experiments, Technometrics 33 (2) (1991) 161–174. [35] M. Sung, M. Gleicher, S. Chenney, Scalable behaviors for crowd simulation, in: Computer Graphics Forum, vol. 23, Wiley Online Library, 2004, pp. 519–528.