A terminating system may or may not have a transient period following the start of the simulation. For example the queues which form in a bank may build for half an hour after the doors open and then remain stable (in the sense discussed above) for the remainder of the day. If you are interested in the time an average customer spends waiting in the queue you may wish to exclude this transient period. However in the more usual case the arrival rate of customers varies throughout the day and no initial transient can be recognised and averages are most usefully defined over the whole days trading.

Whatever period is chosen for analysis the basic procedure is the same. The simulation run is replicated n times using a different random number stream in each case to ensure independent trials. If we assume the variation in any output parameter Y is approximately normally distributed we can compute an estimate of the mean using the standard formulae

where the Y_i are the means of the individual replications. The sample variance is given by

This gives us an approximate 100(1 - a ) percent (0<a <1) confidence interval for the expected value of

Notice that to halve the uncertainty in the mean you need to average four times as many runs.

For the data in table 5.1 we can estimate the expected average delay in the queue to be given by Y_m(10) = 10.2 and the variance S²(10) = 54.48 This gives a variance of the expected value of S²(10)/N = 5.44.

If we accept a confidence level of 90% and apply a 'two tailed' Student's t test to this data we find :

In other words we have 90% confidence that the expected waiting time in the queue will be 10.2 ± 4.2 minutes.

n

In example 11.4 we were able to estimate the precision obtained from replicating the basic experiment 10 times. The value, in that example, turns out to be quite large. The question naturally arises as to how many replications we have to perform to achieve a specified precision ?. Our estimated mean Y_m is only an approximation to the true mean (which is unknown) of m . Let us define the absolute error of our estimate as b where

b = |Y_m - m |

If we now make replications of our experiment until the half length L_1/2 of the 100(1 - a ) percent confidence interval is less than or equal to b then :

In other words Ym has an absolute error of at most b with a probability of 1 - a . Suppose we have constructed a confidence interval for m based on n replications then if we assume that the estimate S²(n) of the population variance will not change appreciably as the number of replications increase an approximate estimate of the number of replications n(b ) required to reduce the absolute error to b or less is given by :

The colon is read as 'such that'.

To find n we can use iteration to increase I by 1 from n until the condition is met.

In practice no simulation results should be quoted or inferences drawn from them unless at least 5 to 9 independent replications have been performed and the variance in the results estimated. Depending on decisions to be made on the basis of the simulation (for example a comparison of two queuing policies) this variance estimate can then be used to decide if further replications are required and if so, how many should be taken.

The analysis outlined above assumes that the distributions are normal or near normal in form. If the distributions are seriously skewed (e.g. some Weibull distributions) then calculated confidence intervals may deviate significantly from the calculated (1 - a ) values.

In a non terminating system we are usually interested in estimating the steady state values of the output parameters. This assumes that we have successfully identified the transient zone and removed its influence on the analysis. This is difficult to do because the onset of steady state has no precise definition and in general there are no simple rules that can be universally applied to confirm the system is operating in the steady state condition.

The following data was obtained from a simple M/M/1 ARENA simulation with an exponential inter-arrival times with a mean of 22 and an exponential process time with a mean of 20.

Over the first 2000 time intervals the data appears to show a classic approach to the steady state. However if we allow the simulation to run for 10,000 time units (Fig. 5.4) we observe that we can no longer make such a simple assumption

Finally if we extend the experimentation out to 500,000 time units we appear to have a very chaotic situation (Fig. 11.5)

n

At first site it appears that there is no evidence for the onset of steady state conditions in Fig. 11.5 but if we construct a moving average for the results a different picture begins to emerges.

We define the moving average over a window of observations w as follows :

As we increase the size of the window 'w' we average out over a larger range of observations. This has the effect of progressively removing higher frequency components in the output parameter. Fig. 11.6 shows this applied to the data of example 11.6. You can see that if we use a window of w = 101 we are left only with the lowest frequencies.

In practice because we are interested in defining the onset of the steady state condition we often do not include the calculation of the moving average over the first w observations. This is the approach taken by the MOVAVERAGE command in ARENA (see section 11.6).

In many cases an even clearer demonstration of the onset of the steady state can be obtained if we compute the cumulative moving average defined as follows :

Fig. 11.7 shows this average applied to the data of example 11.6 and we can now clearly see the onset of steady state behaviour and can predict a mean of 9.6. This agrees very well with the theoretical mean queue size for this M/M/1 queue of

Having identified the steady state region we can now apply the same analysis as used in section 11.4 for terminating systems. To find the confidence interval requires us to replicate the experiment a number of times using different random number streams to ensure independent trials. This procedure is very easy to understand and will give the most reliable results and if computing time and facilities are not limited is probably the best approach to take.

From the data graphed in Fig. 11.7 we can conclude that we are in the steady state region beyond a time of 300,000 units. Table 11.2 shows the results of 9 independent replications of the experiment recording the queue size for the data in the steady state region. From this data we can conclude using the equations of section 11.4 that with a confidence level of 90%:

Average queue length : 9.5 ± 1.0

Maximum queue length : 56 ± 6

n

Experiment replication does however have the drawback that for each experiment the transient region has to be passed and its data rejected before the analysis is performed. There have been many proposals for using the data from a single run to estimate the mean and confidence interval to avoid this problem. We will describe two such procedures chosen because facilities exist within ARENA for implementing them.

We could avoid the inefficiency of many replications if we could form an estimate of the variance of the mean from a single replication. The problem to be overcome is that near adjacent observations of Y are clearly correlated (e.g. large observations tend to be adjacent to other large observations etc.). We can solve this problem if we can assume that observations are covariance stationary, that is that the covariance between any two observations that are 'k' apart is the same where ever these observations are taken in the sequence. Assuming 'n' observations are left after removing the transient data we compute the sample mean, variance and covariances as follows :

The mean value

The sample variance

The covariance for a lag of k

We now estimate the variance in the mean from a weighted sum of all the covariances up to k

To apply the method we need to check the assumption that the observations are covariance stationary. We can achieve this by calculating what is known as the correlation at a lag of k :

To select a suitable value of k we plot a correlogram which is a plot of m_{k ~} k. Fig. 11.8 shows an idealised correlogram.

As the correlation values fall to zero we can assume that stationary covariance has been achieved (at approximately 10 in Fig. 11.8).

The main problem with the method is that sample covariance is a biased estimator of the covariance of the mean and thus a large 'n' and a 'k'<<'n' is required to minimise the bias. Generally speaking if k exceeds n/4 the method will probably produce unreliable results.

Again we take a single replication and discard the observations from the transient region. Now divide the remaining n observations into m batches each of size k (n = mk). If we can assume each batch is equivalent to an independent replication then we can apply the formulae of section 11.4 directly. Now it is clear that the last few observations of batch (k-1) will be correlated with the first few observations of batch k but if the batches are large enough our assumption of independence should be valid. As an estimate of the minimum batch size required we can use the correlogram of the observations and set the minimum batch size to the point where the graph approaches zero. There is no reason why we cannot use batches larger than this minimum provided we retain a total of at least 10 batches to ensure that the Student 't' statistic does not become too large and degrade the confidence half width.

Figure 11.9 shows the correlogram for the data from example 11.6 with a lag of 2,500. Provided that we chose a batch size in excess of 1,200 we can expect a reasonable approximation to independent replications. We chose a batch size of 12,000 which gives us 16 batches between a time of 300,000 and 500,000 time units. We can summarise the results as follows :