Optimizing human activity patterns using global sensitivity analysis

Implementing realistic activity patterns for a population is crucial for modeling, for example, disease spread, supply and demand, and disaster response. Using the dynamic activity simulation engine, DASim, we generate schedules for a population that capture regular (e.g., working, eating, and sleeping) and irregular activities (e.g., shopping or going to the doctor). We use the sample entropy (SampEn) statistic to quantify a schedule’s regularity for a population. We show how to tune an activity’s regularity by adjusting SampEn, thereby making it possible to realistically design activities when creating a schedule. The tuning process sets up a computationally intractable high-dimensional optimization problem. To reduce the computational demand, we use Bayesian Gaussian process regression to compute global sensitivity indices and identify the parameters that have the greatest effect on the variance of SampEn. We use the harmony search (HS) global optimization algorithm to locate global optima. Our results show that HS combined with global sensitivity analysis can efficiently tune the SampEn statistic with few search iterations. We demonstrate how global sensitivity analysis can guide statistical emulation and global optimization algorithms to efficiently tune activities and generate realistic activity patterns. Though our tuning methods are applied to dynamic activity schedule generation, they are general and represent a significant step in the direction of automated tuning and optimization of high-dimensional computer simulations.


Introduction
High-dimensional computer models for simulating real world phenomena have many variables and present a difficult challenge in understanding the relationship between input and output. Known as the curse of dimensionality, a full space analysis of the nature of input-output relationships is NP-complete, scaling exponentially as s n , where s is the number of sample values for each of the n input variables (Rabitz and Aliş 1999). This paper presents an efficient method for determining these inputoutput relationships in high-dimensional models using a combination of global optimization and global sensitivity analysis. We demonstrate our method using a model of human activity and movement.
Human activity and movement patterns are complex and notoriously difficult to model (Berry et al. 2002). Large variations in movement patterns stem from demographic, geographic, and temporal differences. Quantifying the effects of these differences on human activity/schedules provides a difficult but important challenge (González et al. 2008). Realistic human activity and movement models are fundamental components for agent-based infrastructure simulations. These models use human activity patterns to simulate complex systems including epidemics (Eubank et al. 2004;Colizza et al. 2007;Mniszewski et al. 2008;Stroud et al. 2007), traffic (Kitamura et al. 2000), and natural disaster response (Pan et al. 2007). Despite their importance, models typically simplify the complexity of human movement and rely on estimates such as static activity patterns. The static approach results in a Groundhog Day-like effect, where every person performs the same activities day in and day out according to a fixed schedule. Since the schedule cannot be modified based on exogenous events, the schedule will inevitably repeat over some finite time scale.
The level of realism required in a model of natural phenomenon depends upon the scenario being modeled and the questions being addressed (Burton 2003;Burton and Obel 1995). In epidemic modeling, capturing emergent human behavior is crucial for accurately forecasting the spread of disease and the impact of mitigation strategies. Similarly, for modeling disaster response during a natural or man-made event, understanding people's activities before and after the event will help emergency responders allocate resources. Finally, supply and demand modeling of various utilities (e.g., water, electricity, and communications) depends on the population's activities as they move throughout the day. Therefore, capturing realistic activity patterns can help improve modeling efforts and save lives during emergencies.
We have built on the previous body of activity pattern research. Germann et al. presented a study analyzing mitigations for a pandemic influenza in the United States (Germann et al. 2006). In their study, 12-hour schedules were cycled to direct the activities of seven different mixing groups consisting of work, school, day care, play group, neighborhood, neighborhood clusters, and communities. Paleshi et al. performed a similar study that featured stricter mixing patterns according to four coarsely-defined demographic groups: preschool children (ages <1-4), school children (5-18), adults , and seniors (≥65) (Paleshi et al. 2011). In Stroud et al. (2007) and Mniszewski et al. (2008), epidemic simulations rely on static schedules with individuals cycling through nine different activities. Additionally, in contrast to other studies, individuals temporarily deviate from their schedules when ill, and parents stay home with sick children. Weekday schedules were further distinguished from weekends and holidays by replacement of work or school with home for a portion of the population in a study of social contact patterns and their effect on the spread of disease ; it was found that the lack of weekday and weekend activities can greatly overestimate the impact of disease spread. Brockmann et al. moved beyond the realm of static-based schedules by considering random walks as a proxy for human movement based on trajectories of almost 500 thousand dollar bills (Brockmann et al. 2006). González et al. studied the paths of 100 thousand mobile phone users and showed that humans do not behave randomly; rather, they follow simple reproducible spatial patterns (González et al. 2008). All of these models neglect basic human characteristics based on desire, need, and importance that can impact and change schedules accordingly (e.g., getting sick may force a person to go home early from work, or car maintenance may preclude shopping). A realistic human activity and movement model needs to dynamically take these basic human traits into account (Macy and Willer 2002).
We use the Dynamic Activity Simulator (DASim), previously known as Activity-Sim (Galli et al. 2009), that incorporates activity utility and priority to develop schedules for a population of individuals. DASim generates schedules that give each individual close to the maximal utility that complements their priorities for activities. This allows one to design population schedules by specifying priorities and utilities of a variety of activities for any number of demographic groups. Moreover, new schedules can be generated dynamically during a simulation.
Once these schedules are dynamically generated, it is not immediately apparent if they are realistic for a population. In actual populations, we expect demand hours (i.e., the total number of people participating in an activity aggregated over one hour) for certain activities, such as grocery shopping or working, to be stable on any given weekday. For recreation activities or hospital visits, we expect daily demand hours to fluctuate, with possibly a more stable amount of demand hours on a monthly or quarterly timescale. In this way, regularity of demand hours can be required in population's schedules to classify traits of certain activities, thus adding realism to the dynamic schedule generation. We propose to quantify an activity schedule's regularity using the sample entropy (SampEn) statistic (Richman and Moorman 2000). That is, the SampEn of the time series associated with DASim output is used to dynamically adjust schedules to be consistent with regular and irregular activity patterns. By tuning SampEn, one can design schedules comprised of activities that occur with a desired level of regularity.
Tuning the SampEn statistic for a schedule can be posed as a high-dimensional optimization problem. Global sensitivity analysis can be used to reduce the dimensionality of the optimization problem by targeting the input parameters in DASim that control the majority of variation in SampEn. The sensitivity analysis was carried out efficiently through the use of Bayesian Gaussian process regression. Once a low-dimensional set of influential parameters is discovered, a global optimization scheme, harmony search (HS) , is used to tune SampEn and there- fore adjust the regularity of activities in a schedule. We demonstrate that reducing the search space for HS to only influential parameters results in a more efficient search.

Dynamic activities model
DASim is a dynamic parallel agent-based discrete event movement and activity simulator. DASim requires two components to generate schedules: (1) a population with demographic characteristics, and (2) locations with geographic coordinates. DASim can use any population and location data, but the synthetic population we use is based on U.S. census data 1 and includes various demographic characteristics such as age, gender, income, and status (e.g., worker, student, and stay home). In addition, each person has a household consistent with the census data. Locations are derived from the Dun & Bradstreet business directory database, 2 which include addresses and business type. Businesses can be aggregated in a geographic area and may include multiple business types such as a shopping mall. DASim integrates all this information to generate realistic schedules according to the person's preferences and needs. Activities are defined based on the scenarios of interest. For example, they can be general (e.g., home, work, school, shop, social recreation), more specific (e.g., sleep, personal care, breakfast, lunch, food shopping, morning work, afternoon work), or mixed. Subsets of activities are stratified based on different demographic characteristics such as age, school and worker status, and/or gender. Some examples include children (0-5 years old), youth (6-18 years old), workers (19-64 years old), and seniors (65+ years old). In DASim, each demographic group is assigned an activity set comprised of various allowed activities as demonstrated in Table 1. Each activity in each set has associated constraints, a utility function, and a priority function. These controls provide the ability to finely tune activities for each population group.
DASim's utility and priority functions govern activity benefit and importance, respectively. In practice, utility functions influence activity duration, and priority functions influence the order in which activities occur. Utility increases up to a limiting or maximum useful duration. Priority indicates how often an activity is scheduled given Table 2 Parameters used in DASim's utility and priority functions. Each parameter in DASim is shown with a brief description of what it influences and the ranges over which will allow it to vary in our optimization. Figure 1 shows how varying these parameters changes the sigmoid function in practice. Refer to Joh et al. (2001) for more details on these parameters Parameter Description Range The α parameter determines the location of the sigmoid function along the X-axis. In the case of utility, larger α values result in longer-duration activities. Larger α parameters in priority functions result in less frequent activities.
[0, 86400] The β parameter determines the slope of the sigmoid function. For utility functions, larger β values limit the range of desirable activity durations. For priority functions, larger β values increase how rapidly an activity becomes important once enough time has elapsed since last execution.
[0, 1] The γ parameter determines the inflection point in the curve. It defines how quickly the sigmoid's concave region transitions to the convex region. For utility functions, larger γ values result in a more drastic concave to convex transformation, which in turn reduces activity duration range. Larger γ values in priority functions result in a larger range of time between activities.
[0, 10] the longest possible time between activity executions. So, the utility is a function of activity duration, d ≥ 0, while priority is a function of activity start time, t ≥ 0. Utility (U ) and priority (P ) functions are represented in DASim by the sigmoid function presented in (Joh et al. 2001), where α {u,p} , β {u,p} , and γ {u,p} are activity-specific parameters that determine the function's offset, slope, and inflection point, respectively. Table 2 describes in more detail how these parameters affect utility and priority. For a more detailed analysis of these parameters and what they mean in practice, the reader is referred to Joh et al. (2001). U and P vary over the interval [0, 1]. To change a dynamically-generated schedule in DASim, we vary the six parameters (α u , β u , γ u , α p , β p , γ p ) for each activity and for each demographic group. Figure 1 demonstrates sample utility and priority functions for several activities using different parameter sets. A schedule is defined as a set of activities, where each activity has a specified minimum and maximum duration, start and end window, utility and priority functions, maximum travel time, and probability the activity will be performed on a weekend. Activities are scheduled in windows of time (e.g., a 24-hour window means that each person schedules activities 24 hours in advance). A schedule s * is generated by maximizing an objective function that balances the utility of an activity against the priority of all activities and the time it takes to travel to each activity in order to rank schedules, Utility is a function of the duration of an activity. Activity utility is a representation of how desirable or beneficial an activity is. In this example, the maximum utility derived from sleeping is at about 8 hours. This means that sleeping for 8 hours is more beneficial than sleeping for only 6 hours, but there is no benefit to sleeping beyond 8 hours.
(b) Presents what the worker demographic group's priority functions might look like. Priority is a function of the time since an activity was last executed. Activity priority is a representation of how important the activity is. In this example, it becomes important for a person to sleep once it has been about 15-20 hours since they last slept. Additionally, a person needs to work once it has been approximately 15 hours since they last worked. The sigmoid function used and graphs shown here are adapted from data presented in Joh et al. (2001) Here, s is a schedule in the set of all possible schedules S, N is the number of activities in schedule s, U a i (d i ) is the utility of activity a i of duration d i , C is the priority multiplier, B is the number of all possible activities from which the agent can choose, P a i (t r ) is the priority of activity a i at time t r , D is the travel time multiplier, T T a max is the maximum travel time for activity a, and T a i is the travel time for activity a i . The two parameters in the objective function, C and D, weigh the importance of the priority function and travel time constraints, respectively. C and D are global parameters and apply equally to all activities for all demographic groups. The three parameters in each of the utility and priority functions are local parameters set on a per-activity basis.
Schedules are designed using the local search metaheuristic (Lourenço et al. 2003), similar to the method used in Joh et al. (2001). The local search algorithm iteratively adds new activities to a schedule or randomly selects an operator to ap- Table 3 Operators used by DASim when designing schedules. The local search metaheuristic randomly selects operators in order to create a new valid schedule from the current schedule. The objective function then compares the current schedule against this new schedule, as described in Fig. 2 Operator Description Delete-IncreaseDuration Delete an activity and increase the durations of the activities surrounding the deleted activity.

AdjustDuration
Change the duration of two consecutive activities.

Substitute
Replace an existing activity with another.
DecreaseDuration-Insert Decrease the duration of two consecutive activities and insert a new activity.
Append Add a new activity.
ply to the schedule from the operators presented in Table 3. Activities are selected randomly from a set specific to each demographic group with probability weighted by priority. Activity duration is chosen using the specified time constraints. Travel distance is calculated relative to the location of the previous activity and is calculated as Euclidean distance, not as road or travel distance. To calculate travel time, we divide the travel distance by the average speed (fixed at 16 m/s). The objective function is used to evaluate proposed schedule changes. The schedule for a time window is complete when full (i.e., when there is no unaccounted-for time in the individual's scheduling window) and a fixed number of optimization iterations have been completed. In our experiments, we use 10 iterations of the local search algorithm during the optimization step. A larger number of optimization steps allows local search to design slightly better schedules, but this comes at the cost of increased compute time. Figure 2 presents a diagram describing the local search process.
In this study, we concentrate on a randomly-generated 10-person test population. Each of the 10 people in the test population is allowed to create schedules from an activity set comprised of two activities. The first activity is allowed to be between 1 and 24 hours long (allowing for a variety of short-or long-duration activities, such as personal care, shopping, and medical appointments). The second activity is set to be between 4 hours and 10 hours (forcing longer-duration activities, such as work, home, and sleep). The weekend factor for both activities is 1.0 (indicating that the activities are equally likely to occur during the weekend as they are during the week). The maximum travel time for each activity is fixed at 2 hours. Activities are allowed to start and end at any point during the day.

Sample entropy
Certain human activities occur with a high degree of regularity (e.g., working, going home), while others occur more erratically (e.g., medical treatment, social recreation) (Bhat et al. 2004;Kitamura and Hoorn 1987;Kitamura et al. 2006;Schlich and Axhausen 2003). Here, we develop a procedure to choose DASim parameters (α {u,p} , β {u,p} , γ {u,p} , C, D) that ensure spontaneity or regularity in an activity. We use the sample entropy (SampEn) statistic to detect regularity in a time series associated with a schedule. The loop begins with a (possibly empty) initial schedule. Using the local search metaheuristic described in the text, a new valid schedule is constructed. The objective function compares the new schedule to the old schedule, and the better schedule is then used. This process continues until the scheduling window is full and the optimization timer is exhausted SampEn was first introduced by Richman and Moorman (Richman et al. 2004;Richman and Moorman 2000) in response to Pincus' seminal work on approximate entropy (ApEn) (Pincus 1991). Entropy quantifies the amount of order or disorder in a system. Ordered systems yield low entropy while disordered or chaotic systems yield high entropy. For a time series, this usually means that a low entropy system will have repeated changes or will remain constant, while a high entropy time series will have unpredictable changes that are highly variable. ApEn was originally developed to analyze regularity in medical and biological time series, specifically neonatal heart rates. It is still commonly used in medical literature (Goldberger et al. 2002;Hornero et al. 2005Hornero et al. , 2006Pincus and Goldberger 1994;Varela et al. 2003) and has also been applied to a variety of other fields including finance (Pincus and Kalman 2004) and human factors engineering (McKinley et al. 2011). SampEn improves on ApEn in several ways; most notably, it is a less biased statistic and requires about half the computing time (Richman and Moorman 2000).
SampEn computes the conditional probability that if a finite time series repeats itself within a tolerance r for m points, then it will also repeat itself for m + 1 points, without allowing self-matches (Lake et al. 2002). Small values of SampEn (values close to zero) indicate signal regularity (i.e., an ordered system), while relatively larger values indicate less regularity (i.e., a more disordered system). SampEn is still a comparative measure; there is no single threshold above which we may say that any arbitrary signal is irregular. It must be judged relative to the problem being addressed.
In our simulations, SampEn is used to quantify regularity of demand hours for activities on an hourly basis (i.e., m = 1). It is common practice to set r equal to some fraction of the standard deviation (σ ) of the data being analyzed, allowing measurements on datasets with different amplitudes (Richman and Moorman 2000); thus, we set r = 0.2σ , where σ is computed from DASim's demand hours output. We use the SampEn implementation written in C provided by PhysioNet. 3

Global sensitivity analysis
We perform a global sensitivity analysis on the SampEn values computed from 12week DASim simulations with respect to the input parameters for the priority and utility functions. DASim outputs demand hours on an hourly basis for each activity, which represent the total number of people participating in an activity aggregated over one hour. For our 12-week simulation period, DASim outputs 2,016 demand hour data points. Figure 3 shows a one-week sample of DASim output (168 demand hour data points). Note how regularity is evident for home and work activities on a 24-hour cycle.
Note that we will also be taking multiple samples of θ parameter sets to construct a statistical model of the SampEn. From M samples of θ parameter sets, we form the Fig. 3 Sample DASim output. DASim outputs "demand hours" on an hourly basis. Demand hours represent the total number of people participating in an activity at a given time. This example shows output for four activities over the course of two weeks for the Minneapolis-Saint Paul region in Minnesota. The DASim parameters used for this figure were selected by hand to agree with a desired schedule of activities (e.g., regularity, peaks and valleys of a certain size, etc.). Some activities (e.g., home and work) occur with obvious regularity while others (e.g., social recreation) occur more sporadically. We use the sample entropy statistic to quantify regularity M × 14 sample matrix Θ whose rows are the samples of the θ parameter sets. We will then use the notation Θ i,j to refer to the j th parameter in the i th sample with i = 1, 2, . . . , M and j = 1, 2, . . . , 14. A single subscript will refer to a row of Θ, so Θ i is the i th sample parameter set, i = 1, 2, . . . , M.
The dynamic scheduling and SampEn computation define the function with entries corresponding to each activity. We calculate Sobol-Saltelli sensitivity indices (Oakley and O'Hagan 2004;Saltelli 2008) for SampEn 1 (Θ) and SampEn 2 (Θ); here, we explain this process for just one of these. We denote the scalar Se, for one activity, without index, as Se = SampEn n (Θ) for n = 1 or 2 (this is done for brevity in the following formulas). First, we specify an allowable range for each of the parameters, θ j ∈ [θ − j , θ + j ], and consider θ j as a uniformly distributed random variable on [θ − j , θ + j ]. This makes the SampEn for each activity a random variable with variance determined by each of the ranges of θ j and its dependence on each of these variables.
We compute first order Sobol-Saltelli sensitivity indices, defined as: where V (Se) denotes the variance and E(Se|θ j ) denotes the expectation of the conditional random variable Se|θ j . In the variance of the conditional expectation, V (E(Se|θ j )), the expectation integral is taken over all variables except θ j , with the j th variable fixed, and the variance is an integral over just θ j . These sensitivity indices represent the fraction of the variance in Se that is attributed to variation in θ j . An equivalent interpretation of S j is the expected fraction by which the variance in Se will be reduced, if the value of θ j is fixed. The S j rank the importance of each variable, θ j , in terms of how much change in Se is present when θ j is varied within [θ − j , θ + j ]. However, the first order indices do not provide a complete ranking of parameter importance when simultaneous variation in sets of variables is allowed (Homma and Saltelli 1996). To quantify importance of a parameter while accounting for its interaction with other variables, we calculate total effect sensitivity indices: The S T j represent the expected fraction of the variance in Se remaining, if all parameters except θ j are fixed. This then accounts for how the remaining variance due to θ j can change, if θ ∼j is fixed at different values.
To rank the importance of each variable with respect to the variation in Se, we examine the entire set (Saltelli 2008): The sensitivity indices have some desirable properties when applied to ranking parameters with respect to their influence on the variance of an output. If a variable does not influence the function at all, S j = 0, and if a variable does not have any interaction with the other variables, S j = S T j (Sobol 2001). In all situations, we have (Sobol 2001): Regardless of the utility of these sensitivity indices, they can be difficult to interpret since they are dependent on the distribution of the input parameters. Changing the interval for the parameter θ j , [θ − j , θ + j ], changes the indices S j and S T j . Moreover, since this interval affects V (Se), changes to the interval θ j may affect the sensitivity indices of other parameters. This is due to the global nature of the Sobol-Saltelli sensitivity indices and may cause interpretation difficulties due to parameter interdependencies.
A traditional Monte Carlo approach to compute the sensitivity indices is computationally expensive due to the repeated/iterated terms such as V (E(Se|θ j )). A variety of approaches have been suggested to bring down the computational cost (Homma and Saltelli 1996;Marrel et al. 2009;Oakley and O'Hagan 2004;Saltelli 2002;Saltelli et al. 1999). We compute approximations to the sensitivity indices using a statistical surrogate model (Marrel et al. 2009;Neal 1997;O'Hagan 2002, 2004), or emulator, for the function Se(θ ). The emulator uses a Gaussian process regression (Higdon et al. 2008;Marrel et al. 2009;Neal 1997;Williams et al. 2006), which consists of fitting a Gaussian process Se g (θ ; η) to samples of Se(θ ) taken at different θ parameter sets specified by the rows of the M sample matrix Θ.
In the Bayesian regression approach, Se g is constructed from samples of the output Se i = Se(Θ i ), i = 1, 2, . . . , M. The mean and covariance of Se g are defined so that realizations of the simulated values have a maximized posterior probability given a prior distribution on the form of the covariance. The form for the covariance is specified so that when evaluating at a new parameter set, θ * , the variance of Se g (θ * ) increases for θ * further from the samples in the matrix Θ and goes to zero, if θ * lies in this sample set. The mean of Se g (θ * ) is related to the sampled values so that it is equal to Se i for θ * = Θ i . Thus, Se g (θ ; η) is an interpolant of the sample values.
Sensitivity indices of E η (Se g (θ ; η)) can be computed quickly once Se g is constructed from a sample set. We refer to Marrel et al. (2009), Oakley andO'Hagan (2004) for this computation. To construct the Gaussian process and to compute the sensitivity indices, we used the Los Alamos GPM/SA code 4 (Higdon et al. 2008;Williams et al. 2006).

Global optimization
Our goal is to find values for each of the parameters in θ for which SampEn, for the given activities, is either minimized (for increased regularity in scheduling) or maximized (for increased spontaneity). Optimizing over the complete 14-dimensional parameter space can be costly. Note that this 14-dimensional space is only for two activities; each additional activity adds 6 new parameters. Therefore, analyzing five activities would require optimization over a 32-dimensional space, which is computationally expensive for updating a schedule dynamically.
We use the global sensitivity indices to reduce the dimensionality of the optimization problem and identify parameters that contribute very little to the variance of SampEn. In an optimization step, these parameters are then fixed, and the remaining parameter space is searched using a global optimization procedure. If the number of parameters to which SampEn is sensitive is small, this can potentially result in a cheaper optimization procedure.
Schedules may be generated so that each activity has a desired level of regularity/irregularity by maximizing a single objective function, J (θ), involving the Sam-pEn statistics for each activity in the schedule. We define the objective function for a schedule of N activities, A 1 , A 2 , . . . , A N , by Here, the desired levels of SampEn for each activity are denoted by L i and weights, w i , are associated with each activity to control the importance of each term in the maximum of J (θ). It is important to note that we include the square of the absolute value in our objective function so that J (θ) is smooth. Maximization of these types of objective functions can attain specific goals allowing for more specificity in schedule design with regards to mixtures of regular and irregular activities. For instance, in a two-activity schedule we may choose w 1 = w 2 = 1 and L 1 = L 2 = 0 to obtain the objective function Maximization of (9) generates schedules where both activities have a high SampEn and, therefore, have irregular activity demand hour time series for both activities. Alternatively, taking w 1 = 1, w 2 = −1, and L 1 = L 2 = 0, we get Maximization of (10) will generate schedules in which A 1 has highly irregular demand and A 2 has very regular demand. More specific conditions can be met by specifying non-zero levels of SampEn for each activity. Setting w 1 = −1, L 1 = 0.9, w 2 = −0.5, and L 2 = 1.5 we get When maximizing (11) the contribution from the term involving SampEn 1 (θ ) has twice the effect of the contribution from the term involving SampEn 2 (θ ). Therefore schedules will be generated with SampEn 1 (θ ) ≈ 0.9, SampEn 2 (θ ) ≈ 1.5, and SampEn 2 (θ ) farther from 1.5 than SampEn 1 (θ ) is from 0.9. We use the harmony search (HS) global optimization algorithm ) to explore the parameter space. HS is a metaheuristic search algorithm, inspired from the improvisation process of jazz musicians, that optimizes (minimizes or maximizes) a certain objective function. Recently, HS has been successfully applied to a variety of problems including water distribution network design (Geem 2006b), parameter estimation , combined heat and power economic optimization (Vasebi et al. 2007), and even sudoku solving (Geem 2007). In many cases, it has been shown to outperform other commonly used search algorithms, such as simulated annealing (Kirkpatrick et al. 1983), tabu search (Glover 1989(Glover , 1990, and evolutionary algorithms (Bäck and Schwefel 1993).
In HS, sets of parameters (referred to as a harmonies) are randomly chosen (improvised) until the harmony memory is filled. A new harmony is improvised according to a set of rules: each parameter (note) may be chosen via random selection or memory consideration with an optional pitch adjustment (adjusting a parameter up or down slightly). The goodness of the new harmony is computed (in this case, the sum of the SampEn statistics for each activity), and if the harmony is better than the worst harmony stored in the harmony memory, the new harmony replaces the previously stored value.
HS features five main parameters: max_imp determines the maximum number of improvisations (iterations), hms is the harmony memory size (the number of best harmonies that should be remembered), hmcr is the harmony memory consideration rate (how often a note is chosen via memory consideration as opposed to random selection), par is the pitch adjusting rate (how often pitch adjustment is invoked), and mpap is the maximum pitch adjustment proportion (size of the perturbation).
A number of improvements and changes have been suggested since HS' first introduction. One change added the notion of ensemble consideration, an operation that considers relationships between decision variables (Geem 2006a). Another modification, dubbed improved harmony search, dynamically modifies the par and mpap parameters as the search progresses (Mahdavi et al. 2007). Global-best harmony search removes the mpap parameter altogether by altering the pitch-adjustment step so that values are drawn from the best harmony in the harmony memory (Omran and Mahdavi 2008). Most recently, a parameter setting-free variation was introduced that dynamically modifies both hmcr and mpap as the search progresses (Geem and Sim 2010).
For this study, we implemented the original HS algorithm in Python. The source code has been open-sourced and is available on GitHub. 5 At the start of our HS optimization for DASim, C, D, β u , and β p are allowed to vary in the range [0, 1], while α u and α p are allowed to vary in the range [0, 86400]. The parameters γ u and γ p are allowed to vary in the range [0, 10]. Notice, there is no need to normalize all inputs to a common range since the sensitivity indices rank the inputs relative to their ranges. HS is then combined with global sensitivity analysis to reduce the dimensionality of the search space, which is done iteratively as follows: Harmony Search with Global Sensitivity Analysis Algorithm 1. Provide allowable intervals for each parameter θ j ∈ [θ − j , θ + j ], j = 1, 2, . . . , 14. 2. M samples of SampEn are taken at different parameter sets in the sample matrix Θ M×14 . 3. Samples used to construct a Gaussian process emulator, Se g (θ ; η). 4. Sensitivity indices, {S j , S T j }, j = 1, 2, . . . , 14, are computed from Se g (θ ; η). 5. A subset of parameters, (θ k 1 , θ k 2 , . . . , θ k d ), with high sensitivity values (see Fig. 4) are chosen on which to perform HS. The remaining parameters are fixed (note that we arbitrarily fix them at the mean value of their interval). Here, we use notation for an arbitrary subset of distinct parameter subscripts of size d ≤ 14, {k 1 , k 2 , . . . , k d } ⊂ {1, 2, . . . , 14}. 6. HS is performed over the parameter subset to maximize a given functional of SampEn statistics for each activity. , we see that the α-parameters are much more sensitive than the other parameters. This implies that variations in α-parameters account for the majority of the variation in sample entropy values. One can also see some contribution to the variance by the γ -parameters. Moreover, since the first order sensitivity indices are near zero and the total effect indices are near one, we can tell that the interaction of the α-parameters is very nonlinear. In practice this means that a schedule's regularity is most sensitive to the location of the utility and priority function and less sensitive to its shape Each 12-week simulation of DASim for the 10-person test population takes approximately 5 seconds wall time to complete. We initialize max_imp to 2000, hms to 50, hmcr to 0.75, par to 0.5, and mpap to 0.25. HS consistently converged to solutions of approximately the same fitness over many test runs, each with initial harmonies selected uniformly at random. As a result, we determined that a parameter sweep of the HS parameters was unnecessary.

Results
The global sensitivity analysis (Fig. 4) shows that the offset parameters, (α u1 , α p1 , α u2 , α p2 ), have the largest effect on the variation of the sample entropy for both activities. Thus, these α parameters have the most impact on regularity. Recall that, as seen in Table 2, α u and α p control the activity duration in the utility and activity frequency in the priority, respectively.
We adjusted schedules to consist of two irregularly performed activities. This was done by maximizing (9), adjusting only four parameters, (α u1 , α p1 , α u2 , α p2 ), using the HS global optimization algorithm. The results of the four-parameter space were compared against tuning the entire 14-dimensional space. In Fig. 5, we show that, for small numbers of HS iterations (i.e., less than 350), the four-dimensional subspace search performs better on average; we can reach much closer to the maximum Sam-pEn in fewer iterations than a search over the whole parameter space. While running 500 iterations of HS over the whole parameter space will result in a better maximum SampEn, our results show that HS over the four-dimensional space will reach 90 % of the maximum SampEn with fewer than 100 iterations. Therefore, the search space should be chosen based on computational requirements.
We performed random sampling over the entire 14-dimensional parameter space and compared the variance in the SampEn for each activity against only varying the α-parameters. Our results show that the variation in SampEn caused by only varying the α-parameters was responsible for about 99 % of the variance in SampEn when the entire parameter set was allowed to vary. This result was consistent for each activity. This shows that our sensitivity analysis with the emulation gives realistic results and that optimization over the four-dimensional parameter space will suffice to approximate the minimum or maximum of the sample entropy or a functional thereof.
The maximization of SampEn over the α-parameters creates a schedule with a great deal of spontaneity. In addition to maximizing the sum of both SampEn statistics, we preferentially maximized and minimized each activity individually, ignoring the other activity. These SampEn-minimized and -maximized schedules, along with a schedule that has SampEn equal to the mean of the minimized and maximized Sam- Fig. 6 Demand hours resulting from SampEn optimization. These plots show demand hour output for the two test activities, denoted A 1 and A 2 , over the course of two weeks for the Minneapolis-Saint Paul region in Minnesota for the 10-person population. The top plot shows demand hours for A 1 , and the bottom plot shows demand hours for A 2 . The parameters for DASim were chosen using HS, which was directed to either maximize or minimize SampEn preferentially for one activity with m = 1. HS coupled with sensitivity analyses allow us to design order and disorder into DASim after first reducing the parameter search space. In both plots, as SampEn decreases, demand hour variability and range visibly reduces. Note that for activities considered on an hourly basis (i.e., when m = 1), the most regular time series (i.e., minimal SampEn) would be a constant line while the most chaotic time series (i.e., maximal SampEn) would be far from constant. We can see that, for both activities, as SampEn increases, regions of constant demand hours become more rare pEn schedules, are shown in Fig. 6. We see a visual difference: DASim output for a maximized SampEn schedule is more variable over a larger range than minimized SampEn schedules. Also, when SampEn is minimized, regions of constant demand hours are more prevalent, which is to be expected for activities considered on an hourly cycle (i.e., when m = 1).

Discussion
This study focuses on schedule realism in a human activity model, but the methods presented here are generic and can be applied to a variety of other problems where a specific property in a high-dimensional model is desired. These types of high-dimensional tuning/optimization problems are ubiquitous in modern complex computer simulations. Thus, there is a significant need for methods of automatic tuning that incorporate systematic dimension reduction. Our combination of global sensitivity analysis and a global optimization method is effective for the application presented here. Additionally, it is sufficiently general to warrant application in many other areas.
Dynamic scheduling for synthetic populations is necessary to make simulations of human behavior phenomena more realistic. The dynamic scheduling program DASim was designed to aid in large-scale agent-based infrastructure simulation (e.g., transportation and epidemic modeling). DASim can generate schedules that are different over demographics and change in response to events, such as disease outbreaks and non-pharmaceutical interventions.
To evaluate the realism of a dynamically-generated schedule, we must select metrics on which it should be evaluated. We presented a method for tuning a dynamic scheduling model for schedule regularity, which we quantify using the sample entropy (SampEn) statistic applied to population demand hours. Adjusting the SampEn statistic requires working with a high-dimensional optimization problem. We used global sensitivity analysis and statistical surrogate models to significantly lower the dimensionality of the search space. A global optimization algorithm, harmony search (HS), was used to efficiently tune the degree of regularity of a schedule.
Some of the major results of our study include: -Demand hour regularity of activities over a population can be controlled by tuning the SampEn statistic. -DASim parameters that most influence the SampEn statistic can be identified using global sensitivity analysis combined with a statistical surrogate model. We determined that the α parameters in the utility and priority functions have the largest effect on the variation of the sample entropy of an activity. -DASim parameters that result in close to optimal (i.e., minimized/maximized) SampEn values can be discovered using HS. Furthermore, this can be done efficiently with many fewer iterations by searching a parameter subspace determined by global sensitivity analysis first (just the α parameters in this study).
While we have shown how to reduce the search space and computation time when analyzing parameter importance under a particular metric, this process still takes a significant amount of compute time. We tuned our parameters in a reduced problem environment, using a 10-person population. Although this approach works for the measure of regularity discussed in this paper, it may not work for more complex measures of interest. Our initial search space of 14 dimensions is still relatively small; some simulations may have many tens, hundreds, or even thousands of dimensions. Understanding parameter importance and interactions in such highdimensional spaces may prove difficult or even impossible in some instances using our methods.
Our analysis is based on hourly regularity for demand hours of schedules. Many other granularities may be desirable. For example, work may be regular every 12 hours. Some studies suggest that the size of the dataset be at least 10 m and preferably at least 30 m in the approximate entropy (ApEn) algorithm (Pincus and Goldberger 1994). While this is certainly possible for small values of m (recall that m = 1 in this paper), larger values of m quickly become problematic (e.g., work would require at least 10 12 demand hour data points). Alternative measures of regularity may be considered for larger values of m.
We are considering other evaluative measures to quantify additional properties, beyond regularity, of a schedule's realism. Here, we analyzed measures of regularity of time usage and found that it is controlled by a small set of the defining parameters in the model. Another possibility would be to look at quantifying the efficiency of a dynamically-generated schedule in terms of location usage, whether an individual's schedule is geographically arranged in a sensible way given his or her current location. One could also look at the total percentage of time spent on an activity. Evaluation of a measure of each of these effects would lend testable realism to a generated schedule for a population. The use of statistical emulation, global sensitivity analysis, and optimization as demonstrated here would then allow for efficient tuning of these measures.
For models that rely on human activity patterns and movement, such as disease and infrastructure models, capturing realistic activity patterns is crucial for decision support. Therefore, new techniques such as the ones proposed here are needed for analyzing high-dimensional problems. However, more research still needs to be done related to efficiently solving these problems computationally and understanding human activity patterns and behavior. James M. Hyman is the Evelyn and John G. Phillips Distinguished Chair in Mathematics at Tulane University. He is the past president of the Society for Industrial and Applied Mathematics (SIAM) and former group leader of the Mathematical Modeling and Analysis Group at Los Alamos National Laboratory. His research interests include mathematical biology, the analysis of nonlinear dynamical systems, and the numerical solution of differential equations. He holds a Ph.D. in mathematics from the Courant Institute of Mathematical Sciences at New York University.