| t: | time (0 
t
T); |
| A(t): | the number of agents on the landscape at time t. |
Barry G. Lawson and Steve Park
{bglaws, park}@cs.wm.edu
Department of Computer Science
The College of William & Mary
Williamsburg, VA 23187
Epstein and Axtell (1996) define a representative artificial society model containing comprehensive rules that simulate many social processes. The authors' artificial society is a model rich with research topics of interest and provides an important basis for social science research via computer simulation. Our work was initially motivated by a careful examination of the results presented by Epstein and Axtell and our attempt to reproduce those results. We were particularly motivated by our inability to reproduce the results Epstein and Axtell obtained when agent movement and reproduction are combined, leading us to consider an alternative modeling approach. Consistent with the request for interdisciplinary coordination by Müller et. al. (1998), in this paper we offer an evaluation of, and suggestions for improvements to, the simulation techniques used in typical artificial society discrete-event simulation models.
Generally, artificial society models employ simple synchronous time evolution, i.e., all agent events on the landscape are synchronized to a simulation clock that evolves deterministically at a fixed rate. An alternative approach is to permit each type of event to occur randomly at its own characteristic rate; in this way, time evolves asynchronously. While we agree that synchronous behavior is correct for some applications, such as seasonally motivated movement, many other applications are modeled more accurately using asynchronous time evolution. For such asynchronous applications, forcing events to occur synchronously is unnatural and can produce simulation artifacts. In addition, forced synchronous time evolution produces an ambiguity in the order in which events should occur. Because all events occur randomly at their own characteristic rates, the use of asynchronous time evolution removes any such ambiguity.
In this paper, we demonstrate that very different behavior can be observed in a typical artificial society discrete-event simulation model if agent events occur asynchronously rather than synchronously. For the results presented, we use a variation of the artificial society model from Epstein and Axtell, with only resource growback, agent movement and agent reproduction. The states and attributes for the landscape cells and agents are real-valued. This variation of Epstein and Axtell's model provides functionality sufficient to observe behavior of interest yet remains simple enough to emphasize the difference between synchronous and asynchronous time evolution. We show that the use of synchronous time evolution can introduce simulation artifacts and produce output that is very different from the output obtained using asynchronous time evolution.
Asynchronous time evolution in the artificial society model is implemented using a next-event simulation approach1. For large models, the handling of the associated event list can incur a substantial loss in computational performance compared to a simple synchronous time evolution implementation. However, we discuss means of improving the performance using appropriate data structures and algorithms for the event list and we provide results demonstrating that acceptable performance can be achieved. These results confirm that an asynchronous time evolution implementation is an attractive approach for artificial society discrete-event simulation models.
The following is an outline of the remainder of this paper. In Section 2 we briefly describe the artificial society model used in this paper, specifically relating our model to the model defined by Epstein and Axtell. In Section 3 we discuss asynchronous time evolution and the corresponding next-event simulation implementation. In Section 4 we provide results demonstrating the different results possible when using asynchronous versus synchronous time evolution. In Section 5 we discuss means of improving the computational performance of a next-event artificial society simulation model and provide supporting results. In Section 6 we present conclusions and suggest areas of future study. In the Appendix we provide definitions of the random variable models referenced in this paper.
For brevity, we omit a complete description of our artificial society model. Except for the descriptions to follow, the rules used in our implementation are the same as those defined by Epstein and Axtell. A comprehensive description of the artificial society discrete-event simulation model used to obtain the results in this paper is located at
We denote the model attributes as
| T: | the maximum simulated time for the model; |
| X, Y: | the vertical and horizontal landscape dimensions respectively |
| t: | time (0
t
T); |
| A(t): | the number of agents on the landscape at time t. |
When necessary, we use the subscript a = 0, 1, 
, A(t) - 1 to distinguish agents.
Define 
= {0, 1, , X - 1} and 
=
{0, 1, , Y - 1}. Then each
cell in the X 
Y landscape is
distinguished by its (x, y) position with (x, y)
. Each (x, y) cell is
characterized by the attributes
(x, y): |
the resource capacity at cell (x, y); |
(x, y): |
the resource growback rate at cell (x, y) |
| r(x, y, t): | the level of resource at cell (x, y) at time t; |
| d(x, y): | the most recent time of resource depletion at cell (x, y); |
| o(x, y, t): | the occupancy status of cell (x, y) at time t. |
We define the real-valued resource capacity at an (x, y)
landscape cell by the equation
(x, y) = f(x - X/4,
y - Y/4) +
f(x - 3X/4,
y - 3Y/4)
exp(-(x /
x)2
- (y /
y)2
)
= 4.0,
x = 0.3X and
y = 0.3Y.
Although integer-valued but not otherwise explicitly defined, the resource
capacity used by Epstein and Axtell appears to be consistent with this
equation.
Each agent is characterized by the attributes
 : |
the field of view (FOV) attribute; |
 : |
the metabolic rate of resource consumption; |
 : |
sex (male or female); |
 : |
the time of birth; |
 : |
lifespan; |
 : |
the age when reproductive capability begins (puberty); |
 : |
the age when reproductive capability ends ( < ); |
 : |
the gestation period (for females) |
| w(t): | the amount of resource wealth (holdings) at time t; |
| (x, y): | the landscape cell occupied by the agent; |
| m(t): | pointer to the current mate (if any). |
, T. All events of interest must occur
precisely at these time steps. For instance, in the Epstein and Axtell model
all agents attempt to move simultaneously, once each time step. All agents
attempting to reproduce must also do so once each time step. For serial
program execution, this kind of parallel activity is not possible and so it is
necessary to randomize the order in which agents act at each time step.
Without this randomization, the agents will always act in the same order,
introducing a bias that can produce artifacts (Epstein and
Axtell, 1996).
Given the model attribute T, which defines the number of simulated time steps, the artificial society model evolves synchronously in time according to the following pseudocode algorithm.
/* initialize the landscape */
/* initialize the agents */
t = 0;
while (t <= T) {
/* move and reproduce the agents */
/* update the landscape */
/* update the agent list */
/* randomize (shuffle) the agent list */
t++;
}
Note the ambiguity in the order of events inside the while loop. For example, should the landscape be updated before or after the agents move?
To facilitate asynchronous time evolution, we construct our artificial society model using the next-event simulation approach. For a next-event simulation model, we must clearly define the system state, the event types, the simulation clock, and a set of algorithms for each event type that define the state changes that occur when each type of event occurs (Park and Leemis, 1999).
The simulation clock is the real-valued global time parameter t. As
defined in Section 2.1, at any time 0 t T, the state of the
system is defined by the following set of variables4.
; , A(t) - 1.
Accordingly, there are four types of events that can change the state of the system.
Because there is a gestation period, the state of the system can change when agents mate and when a new agent is born. For this reason, the reproduction rule, as defined in Section 2.2, is modeled by two distinct types of events. A mating event type initiates the gestation period if the female agent becomes pregnant. Since a mating event is an attempt to reproduce, the mating event does not guarantee pregnancy. A birth event type signals the end of the gestation period when the new agent is born (if possible). Both of these event types can change the state of the system as described above.
) where
is the rate of occurrence of that event type. Equivalently, each event type is
modeled as a stationary Poisson process with rate . The Poisson rate for the agent movement process is = 1.0
5. The mating event type can either be coupled with the
movement/consumption event type or modeled as a separate stationary Poisson
process. If the two event types are coupled, a mating event occurs for an agent
each time a movement/consumption event occurs for that agent. For all the
results in this paper, if the two event types are uncoupled, the Poisson rate
for agent mating is also = 1.0.
For each (x, y) landscape cell, as time evolves d(x, y) is defined by the last time an agent departed that cell. At time t > d(x, y), we can compute the current level of resource at cell (x, y) by the equation
(x, y), (x, y)(t - d(x, y))
}.
/* initialize the landscape */
/* initialize the agents */
t = 0.0;
while (t <= T && A > 0) {
event = GetNextEvent();
t = event.time;
switch (event.type) {
case movement:
/* unoccupy the current cell */
/* select the (x,y) cell within the FOV with maximum resource */
/* occupy the selected (x,y) cell */
/* update the wealth of the agent */
/* schedule the next movement event */
case death:
A(t)--;
/* unoccupy the current cell */
case mating:
/* execute the mating algorithm */
if ( /* the mating is successful */ ) {
/* compute endowments and update the wealths of the parents */
/* update the mate pointers of each parent */
/* schedule a birth event */
/* cancel any parent agent events to occur before the birth */
}
case birth:
/* update the wealths of the parent agents */
/* schedule the next event occurrences for each parent agent */
if ( /* there is an (x,y) cell for the new agent */ ) {
A(t)++;
/* initialize the new agent */
/* occupy the selected (x,y) cell with the new agent */
/* update the wealth of the new agent */
/* schedule the next event occurrences for the new agent */
}
}
}
window at times t, t + 1, as shown in Figure 3-1 where bullets denote typical agent
events.
The movement/consumption and mating events for an agent occur at the same time. Birth and/or death events, if they are to occur, also take place at this time.
The following figures were motivated by an attempt to replicate the results
Epstein and Axtell obtained when combining agent movement with reproduction.
Despite considerable experimentation, we were never able to reproduce the
results in Figure III-4 (Epstein and Axtell,
1996). Since Epstein and Axtell never reveal their model for the
landscape capacity (perhaps for brevity), there is some uncertainty in our
choice of (x, y). However,
we conjecture that the oscillatory behavior shown in Figure III-4 is a
simulation artifact caused primarily by synchronous time evolution. Moreover,
this oscillatory behavior is very much dependent upon initial conditions which
to date we have been unable to determine.
Y = 50 50
landscape. The resource capacity at each cell is given by the equation in Section 2.1, and the resource growback rate for each
(x, y) cell is (x,
y) = 1. The landscape is initially occupied with A(0) = 400
randomly placed agents. The initial wealth of each agent is w(0) = 10.
The attributes for each agent are initialized as random variates according
to6
| : |
Equilikely(1, 6); |
| : |
Uniform(1, 4); |
| : |
Uniform(60, 100); |
| : |
Uniform(12, 15); |
| : |
Uniform(40, 50) for females, Uniform(50, 60) for males. |
= 0. Reproduction and movement are coupled, i.e.,
the mating event type is not modeled as a separate process.
Figure 4-1 depicts the time evolution of the agent population for five different initial random number generator seeds with T = 2500. In general, the output produced using synchronous time evolution is the more oscillatory of the two curves, while asynchronous time evolution produces the more stable darker curves. Epstein and Axtell propose explanations of the large oscillatory behavior in terms of population dynamics. Although in a few atypical cases asynchronous time evolution produced large oscillations in the population7, we believe that, in general, the large oscillations produced by Epstein and Axtell are simulation artifacts produced by the use of synchronous time evolution.
Figure 4-2 depicts the time evolution of the agent
population for an X Y = 100
100 landscape with an initial A(0) =
1600 randomly placed agents. The remaining state and attribute parameters are
the same as in Figure 4-1. Again, undamped
oscillatory behavior is present in the synchronous time evolution output but
not in the asynchronous time evolution output8. While some initial oscillation is present in
the asynchronous time evolution results, the amplitude is damped
considerably. Moreover, for all initial seeds, changes in the output are far
less dramatic in the asynchronous case. That is, unlike the synchronous
results, the asynchronous results are not sensitive to the random variate
sequence of movement and reproduction, as manifested by the choice of initial
seed.
We felt that T = 2500, as suggested by Epstein and Axtell, may not yield a sufficient measure of the long-term time evolution of the agent population. For this reason, Figure 4-3 depicts the time evolution for the same initial conditions as Figure 4-1, but with T = 25,000. In these figures, the highly oscillatory behavior of the synchronous time evolution output persists, while the asynchronous time evolution produces a stable agent population in virtually all cases.
Similarly, Figure 4-4 depicts the time evolution for the same initial conditions as Figure 4-2, but with T = 25,000. In this figure, we see the oscillatory behavior persisting in a majority of the synchronous time evolution output. In contrast, the asynchronous time evolution quickly produces a stable agent population. We present the results in Figures 4-1, 4-2, 4-3 and 4-4 as strong evidence that very different behavior can result from the choice of modeling time evolution as synchronous or asynchronous.
= 1. Again, the darker curves represent
output from the asynchronous time evolution model. We see that not only is the
amplitude of oscillations diminished dramatically using asynchronous time
evolution, but when compared with Figure 4-1 the mean
of the agent population is lowered as well. We reemphasize that we are not
promoting asynchronous time evolution as the best approach for all
artificial society models, but these results show that very different output
can result from modeling time evolution as synchronous or asynchronous.
Given that each agent can execute only one type of event at a time, we can simplify the event list by including only the most imminent event for each agent. In this way, there is one event per agent in the event list. Since the type of the most imminent event for each agent can be determined from the agent data structure, we no longer need the type of the event in the event list. If we use a singly-linked list data structure to represent the event list and sort the events by increasing time of occurrence, we obtain a sample event list like the one shown in Figure 5-1.
Using this approach, the next event to occur in simulated time is always the first element in the linked list. In this way, the GetNextEvent() procedure from Algorithm 3-2 is implemented by simply returning the head of the linked list. However, when a new event is inserted into the event list, a linear search must be executed from the head of the list to determine the appropriate insertion point. Similarly, any search for a canceled event also requires a linear search. For a large number of events, i.e., as the number of agents grows, such linear searches seriously degrade the computational performance of the simulation.
Figure 5-2 depicts sample execution times9 for a standard synchronous implementation of the artificial society model and an asynchronous implementation using the previously described singly-linked list event list approach. Each value displayed is the average time for five replications along with the corresponding 95% confidence interval. The simulation was executed with T = 500 for the increasing landscape sizes shown, each initially populated with 16% agents10. The initial values for agent states and attributes were the same as described in Section 4.1, except initial agent wealth was w(0) = 25. Agents were not permitted to reproduce so that the number of agents monotonically decreased for the duration of the simulation. As shown, as the initial number of agents increases (with landscape size), the performance of the sorted linked list implementation degrades quickly.
As depicted in Figure 5-3, a linear search is avoided through the use of an auxiliary search array. Each element in the search array contains an event time and a pointer to an event in the event list. The search array is sorted by decreasing event time. A typical binary search of the search array, a variation on Henriksen's approach, is used to find the smallest time in the search array greater than the event time sought. Beginning with the event pointed to by the selected search array element, a linear search of the event list in descending order ensues until the desired event time is found. During the binary search, Henriksen's ``pull'' technique (1983) is used to keep the search array pointers evenly distributed throughout the event list. The end of the search array is maintained, and the size of the array is increased if necessary11.
Figure 5-4 depicts sample execution times for an asynchronous implementation of the artificial society model using the sorted linked list approach and the modified Henriksen's approach. Again, each value shown is the average time for five replications along with the corresponding 95% confidence interval. The simulation was executed with T = 500 for the increasing landscape sizes shown, each initially populated with 16% agents. As shown, the data structure and algorithm modifications used in the modified Henriksen's approach yield execution times comparable to those of the standard synchronous implementation given in Figure 5-2. Also shown are the times obtained using the modified Henriksen's approach with agents able to reproduce. Although the number of agents is no longer monotonically decreasing (in time) in this case, acceptable performance is still achieved. These results show that when implemented properly, the loss in computational performance resulting from the use of asynchronous time evolution is minimal in a typical artificial society discrete-event simulation.
Asynchronous time evolution is implemented using a next-event simulation approach. If implemented poorly, a next-event artificial society simulation model can suffer dramatically in terms of execution time. However, through the use of proper event list data structures and algorithms, execution times are comparable to those of a typical synchronous time implementation.
The results in this paper suggest that when judiciously implemented, asynchronous time evolution in an artificial society discrete-event simulation model is an attractive alternative to synchronous time evolution in terms of both the quality of the resulting output and computational performance. Future work involves the systematic study of both standard and novel event list data structures and algorithms on a more comprehensive artificial society model. We also plan to investigate various approaches for executing the artificial society simulation in parallel.
= { x
| a < x < b } 
) if and only if
satisfies > 0 = { x
| x > 0} 
= {a,
a + 1, , b} 
= 0, the
behavior of agents under our reproduction rule is the same as if the gestation
period was omitted from the rule.
> 0, during the gestation period other
agents may occupy all the cells adjacent to the parents. Consequently,
at the time of birth a cell may not be available for the child to occupy.
. Accordingly, the state
description is not minimal.
= 1.0 corresponds to the fixed-increment time step rate of the
synchronous approach. Systematically increasing produces time evolutions similar to those presented in Section 4 but with lower population means. We also
experimented with inter-event times as Uniform, Erlang and
Normal random variables, each with a mean of 1.0 (the Normal was
truncated to positive values). These other distributions for the inter-event
times produce results similar to the stationary Poisson process results.
50 landscape results,
asynchronous time evolution produced no atypical behavior in any of the
results using 100 different initial seeds on a 100 100 landscape. To date we have not performed a more
systematic study of this atypical behavior in either the 50 50 or 100 100 landscape
case.
to
represent the mean and to represent the standard deviation
is consistent with conventional statistical notation and is not to be confused
with the and notations
from Section 2.1.
Author's present address:
Steve Park has been a professor of computer science at the College of William & Mary since 1986 and the chair of that department since 1991. Prior to 1986, he was a member of the research staff at NASA Langley Research Center. He received a doctorate of applied mathematics from North Carolina State University in 1970.
Dr. Park is a member of several professional societies including the Association of Computing Machinery, the IEEE Computer Society and SPIE. His research interests include: the modeling, simulation and performance analysis of sampled imaging systems; algorithms for the reconstruction and restoration of sampled data; and stochastic methods in discrete-event systems simulation. He has published more than 85 research articles in these and related areas.
Author's present address: