1 CHAPTER 2 II. DISCRETE VENT SIMULATION 2.1. Introduction Simulation modelling has been used in a wide range of physical and social sciences and engineering fields, ranging from nuclear fusion to economic forecast to space shuttle design. For different types of situations and systems, different types of models are used. In classifying simulations, there are important distinctions among the types of models that are being simulated, and among the types of program structures that are used to carry out the simulation. 2.2. Deterministic Simulation Models Simulation models may be either deterministic or stochastic (meaning probabilistic). In a deterministic simulation, all of the events and relationships among the variables are governed entirely by a combination of known, but possibly complicated, rules. Deterministic simulations are often used to study the behaviour of physical systems. As a simple example, consider the problem of determining the time required to pay off a $100 loan if you pay $10 per month and the interest rate is 10 percent per year, compounded monthly. In this case, the answer (namely 11 months) can be found very easily, by writing a simple program to ‘‘simulate’’ the history of the loan: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 const rate := .10 / 12 % interest rate per month, expressed as a fraction const payment := 10 % payment per month var loan : real := 100 % current value of the loan (initially $100) var month :=0 % current month put "end of month:" : 15, "loan balance:" : 15, "payment due:" : 10 put month : 10, loan : 15 : 2 loop month += 1 loan *= (1 + rate) const thisPayment := min (payment, loan) loan -= thisPayment put month : 10, loan : 15 : 2, thisPayment : 15 : 2 exit when loan <= 0.0 end loop By running the program, we see that the answer is 11 months: end of month: 0 1 2 3 4 5 6 7 8 9 10 11 loan balance: 100.00 90.83 81.59 72.27 62.87 53.40 43.84 34.21 24.49 14.70 4.82 0.00 payment due: 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 4.86 In this case, it was not really necessary to resort to simulation, since standard mathematical techniques for dealing with series can be used to obtain a general solution to the number of loan periods, P, required assuming the periodic payment is a fraction D of the initial loan, and the interest per period is R:
MODELLING AND SIMULATION 2 ln (D) − ln (D−R) P = , ln (1 + R) which in this case tells us the answer is 10.48 months. The advantage of simulation is that you can still answer the question even if the model is too complicated to solve analytically. For example, consider the trajectory of a baseball under the combined influence of gravity and aerodynamic drag. Once the ball leaves the pitcher’s fingertips, with a given orientation (i.e., the direction that the seams are facing), speed, direction and spin (i.e., the axis of rotation and rotational speed), we can in principle determine its flight path to the catcher by assuming that gravity causes it to fall towards the ground with a constant rate of acceleration, and that at the same time aerodynamic drag causes it to slow down its rate of spin, reduce its speed, and change its direction of flight, at rates that depend on its current orientation, speed and spin. We could use a deterministic simulation to estimate the flight path of the ball by, in effect, ‘‘flying’’ the ball. Periodically, we take a snapshot of the current state of the ball (i.e., its position, speed, spin, etc.), and then calculate all the forces acting on the ball in that state. We then let the ball move ahead a short distance under the combined influence of those forces, assuming that they remain unchanged until we take the next snapshot, recalculate all the forces, and so on. In general, the assumption of unchanging forces causes an error in the resulting solution, so we must limit the distance between successive recalculations of the forces to be small enough for the error to be acceptable. Nevertheless, when the forces acting on the ball are complicated, it can be much simpler to find the flight path from such a deterministic simulation than it would be to solve all the equations to obtain an analytic solution! 2.3. Stochastic Simulation Models In a stochastic simulation, ‘‘random variables’’ are included in the model to represent the influence of factors that are unpredictable, unknown, or beyond the scope of the model we use in the simulation. Throughout the rest of this chapter, we will be discussing stochastic simulation models. In many applications, such as a model of the tellers in a bank, it makes sense to incorporate random variables into the model. In the case of a bank, we might wish to assume that there is a stream of anonymous customers coming in the door at unpredictable times, rather than explicitly modelling the behaviour of each of their actual customers to determine when they plan to go to the bank. It is worth noting here that it is well known in statistics that when we combine the actions of a large population of more-or-less independently operating entities (customers, etc.) the resulting behaviour appears to have been randomly produced, and that the patterns of activity followed by the individual entities within that population are unimportant. For example, all of the telephone systems designed in the last 60 years are based on the (quite justified) assumption that the number of calls made in fixed length time intervals obeys the Poisson distribution. Thus the generation and use of random variables is an important topic in simulation, and we will discuss it in a separate section below. 2.4. Static Versus Dynamic Simulation Models Another dimension along which simulation models can be classified is that of time. If the model is used to simulate the operation of a system over a period of time, it is dynamic. The baseball example above uses dynamic simulation. On the other hand, if no time is involved in a model, it is static. Many gambling situations (e.g., dice, roulette) can be simulated to determine the odds of winning or losing. Since only the number of bets made, rather than the duration of gambling, matters, static simulation models are appropriate for them. Monte Carlo Simulation (named after a famous casino town1 in Europe) refers to the type of simulation in which a static, approximate, and stochastic model is used for a deterministic system. This is in contrast to the baseball example at the beginning of this chapter, where both the system and the model are dynamic and deterministic. Let us now look at an example of Monte Carlo simulation. Consider estimating the value of π by finding the approximate area of a circle with a unit radius. The first quadrant of the circle is enclosed by a unit square (See Figure 1). If pairs of uniformly distributed pseudo-random values for the x and y coordinates in [0, 1] are generated, the probability of the corresponding points falling in the quarter circle is simply the ratio of the area of the quarter circle to that of the unit square: 1 Interestingly, John Von Neumann, a pioneer in computer design and use, used Monte Carlo as the code name for his military project during the Second World War.
3 CHAPTER 2 y (1,1) (0,1) .. .. . . .. .. . . .. .. . . . . . . . . .. .. .. .. .. . . . .. .. . .. .. . . .. . . . . . . .. .. . . . . . . . . . . . . . . . . . . . . . . .. . .. . . . . . . . . . . . .(x,y) . .. . .. . . .. .. .. .. .. .. . .. . . . . . . . . . . . . . . . . . . . . . . .. .. .. .. .. .. .. .. .. . . . . . . . . . . . .. .. .. .. .. .. .. .. .. . . .. .. .. .. . .. . .. .. . . . . . . . . .. . . . . . . . . . .. . . . . . . . . . . . .. (0,0) x (1,0) Figure 1. Monte Carlo Integration of a Quarter Circle. number of points falling in the quarter circle total number of points generated π = 4 The simple Turing program below simulates this situation. Some sample output follows. A good approximation can be achieved by using a large number of points. The convergence is quite slow, though. In the next chapter, we will discuss a faster way through numerical integration. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 randseed(1, 1)% could have used interactive input to initialize, but this is simpler randseed(2, 2) var x, y: real var repeatCount, inCount := 0 put "repetition count: " .. get repeatCount for i: 1 .. repeatCount randnext(x, 1) randnext(y, 1) if x*x + y*y < 1 then inCount += 1 end if end for put "Approximate value of pi (based on ", repeatCount, " samples): ", 4 * (inCount / repeatCount) Repetition count Approximate value of π 100 3.24 1,000 3.16 10,000 3.1608 100,000 3.14508 1,000,000 3.14098 The above program is an example of Monte Carlo integration, by which definite integral of arbitrary but finite functions can be estimated. Consider the following integral: b ∫a f (x)dx
MODELLING AND SIMULATION 4 Such an integral can be partitioned into segments above or below the horizontal axis. Without loss of generality, let us assume f (x)≥0, a ≤ x ≤ b. We can then bound the curve with a rectangle with borders of x = a, x = b, y = 0, and y = y max , where y max is the maximum value of f (x) in the interval [a, b]. By generating random points uniformly distributed over this rectangular area, and deciding whether they fall above or below the curve f (x), we can estimate the integral. 2.5. Program Structures for Dynamic Simulation In the remainder of chapter, dynamic stochastic simulation models will be emphasized. A dynamic simulation program that is written in a general purpose programming language (such as Turing) must: i) keep track of ‘‘simulated’’ time, ii) schedule events at ‘‘simulated’’ times in the future, and iii) cause appropriate changes in state to occur when ‘‘simulated’’ time reaches the time at which one or more events take place. The structure of a simulation program may be either time-driven or event-driven, depending on the nature of the basic loop of the program. In time-driven models (see Figure 2a), each time through the basic loop, simulated time is advanced by some ‘‘small’’ (in relation to the rates of change of the variables in the program) fixed amount, and then each possible event type is checked to see which, if any, events have occurred at that point in simulated time. In event-driven models (see Figure 2b), events of various types are scheduled at future points in simulated time. The basic program loop determines when the next scheduled event should occur, as the minimum of the scheduled times of each possible event. Simulated time is then advanced to exactly that event time, and the corresponding event handling routine is executed to reflect the change of state that occurs as a result of that event. initialize initialize time := startingTime TimeToStop := false loop loop exit when time >= endingTime exit when TimeToStop % Collect statistics and record measurements. % remove entry for next event from event list % Call event routines for each type of event time := NextEventTime % that has occurred by this time. % Collect statistics and record measurements. time += timeIncrement % Call event routine for next event end loop % (end of simulation event sets TimeToStop). % Insert next occurrence of this event type % in event list. end loop a. Time Driven Structure b. Event Driven Structure Figure 2. Structure of the basic program loop. The choice of whether a time- or event-driven program structure should be used depends on the nature of the model. In general, an event driven structure means that the basic program loop gets executed fewer times (since an event is guaranteed to take place every time through the loop). However, there is overhead involved in managing the event list, so that the ‘‘cost’’ (i.e., execution time) of going around the loop once is, in general, substantially higher than it would be for a time-driven structure. Thus, a time-driven structure is generally better for models in which it is possible to find a time increment that is both small enough to distinguish among the times at which events occur that cannot be treated as simultaneous, and also large enough for at least one event to occur at a significant proportion of the times through the basic program loop. If such a compromise is not possible (because events occur at irregular intervals), an event-driven structure is preferable.