Infinite Horizon Simulation Experiments

Learning Objectives

• To be able to analyze infinite horizon simulations via the method of batch means and the method of replication-deletion

Infinite Horizon Simulation

• Infinite horizon: In an infinite horizon simulation, there is no well defined beginning time or initial conditions and no well defined ending time or condition. The planning period is over the life of the system, which from a conceptual standpoint lasts forever. Infinite horizon simulations are often called steady state simulations because in an infinite horizon simulation you are often interested in the long-term or steady state behavior of the system.

• Examples of infinite horizon simulations include:

  • A factory where you are interested in measuring the steady state throughput.
  • A hospital emergency room which is open 24 hours a day, 7 days of week.
  • A telecommunications system which is always operational.

• Infinite horizon simulations are often tied to systems that operate continuously and for which the long-run or steady state behavior needs to be estimated.

Statistical Issues for Infinite Horizon Simulation

• For an infinite horizon simulation, we may want to make one very long run; however, within replication observations are:

• Likely not independent

• Likely not identically distributed

• Likely not normally distributed

Thus, we cannot treat within replication observations as a random sample. This is an issue if we use only one replication in the analysis and for understanding the behavior within a single replication.

Cumulative Average Waiting Time of 1000 Customers

• Customer waiting time for M/M/1 queue. True waiting time is: \(W_q = \dfrac{L_q}{\lambda} = 1.6\bar{33} \; \text{minutes}\)
• As seen in the plot, the cumulative average starts out low and then eventually trends towards 1.2 minutes. \[\dfrac{1}{n} \sum_{i=1}^n W_i \; \text{for} \; n = 1,2,\ldots\]

Initialization Bias and Warm Up Period

• Because the system starts empty and idle, the performance measure is affected by the initial conditions.
• The idea is to delete the data affected by the initial conditions.

Welch Plot Analysis

• Make \(R\) replications. Typically, \(R \geq 5\) is recommended.

• Let \(Y_{rj}\) be the \(j^{th}\) observation on replication \(r\) for \(j = 1,2,\cdots,m_r\) where \(m_r\) is the number of observations in the \(r^{th}\) replication, and \(r = 1,2,\cdots,n\),

• Compute the averages across the replications for each \(j = 1, 2, \ldots, m\), where \(m = min(m_r)\) for \(r = 1,2,\cdots,n\).

  \[\bar{Y}_{\cdot j} = \dfrac{1}{n}\sum_{r=1}^n Y_{rj}\]

• Plot, \(\bar{Y}_{\cdot j}\) for each \(j = 1, 2, \ldots, m\)

• Apply smoothing techniques to \(\bar{Y}_{\cdot j}\) for \(j = 1, 2, \ldots, m\)

• Visually assess where the plot starts to converge

Using the KSL to Perform a Welch Plot Analysis

Within the KSL, Welch plotting is facilitated by using the classes within the ksl.observers.welch package. There are two classes of note:

• WelchDataFileCollector: This class captures data associated with a Response within the simulation model. Every observation from every replication is written to a binary data file. In addition, a text based data file is created that contains the meta data associated with the data collection process. The meta data file records the number of replications, the number of observations for each replication, and the time between observations for each replication.

• WelchDataFileAnalyzer: This class uses the files produced by the WelchDataFileCollector to produces the Welch data. That is, the average across each row of observations and the cumulative average over the observations. This class produces files from which the plots can be made.

Attach the WelchFileObserver to the Response

fun main(){
    val model = Model("Drive Through Pharmacy")
    val dtp = DriveThroughPharmacyWithQ(model, 1)
    dtp.arrivalRV.initialRandomSource = ExponentialRV(1.0, 1)
    dtp.serviceRV.initialRandomSource = ExponentialRV(0.7, 2)
    val rvWelch = WelchFileObserver(dtp.systemTime, 1.0)
    val twWelch = WelchFileObserver(dtp.numInSystem, 10.0)
    model.numberOfReplications = 5
    model.lengthOfReplication = 50000.0
    model.simulate()
    val rvFileAnalyzer = rvWelch.createWelchDataFileAnalyzer()
    val twFileAnalyzer = twWelch.createWelchDataFileAnalyzer()
    rvFileAnalyzer.createCSVWelchPlotDataFile()
    twFileAnalyzer.createCSVWelchPlotDataFile()
    val wp = WelchPlot(analyzer = rvFileAnalyzer)
    wp.defaultPlotDir = model.outputDirectory.plotDir
    wp.showInBrowser()
}

KSL Welch Plot Output

• 10000 observations looks like a good deletion point.

Setting the Warm Up Period

• Use the lengthOfReplicationWarmUp property to specify the length of the warmup period in terms of time.
• Since 1 observation equals 1 time unit, we need 10000 time units for the warmup period.
• This is called the method of replication-deletion.

fun main(){
    val model = Model("Drive Through Pharmacy")
    val dtp = DriveThroughPharmacyWithQ(model, 1)
    dtp.arrivalRV.initialRandomSource = ExponentialRV(1.0, 1)
    dtp.serviceRV.initialRandomSource = ExponentialRV(0.7, 2)
    model.numberOfReplications = 20
    model.lengthOfReplication = 50000.0
    model.lengthOfReplicationWarmUp = 10000.0
    model.simulate()
    model.print()
}

The Batch Means Method

• In the batch means method, only one simulation run is executed. After deleting the warm up period, the remainder of the run is divided into \(k\) batches, with each batch average representing a single observation.
• The idea is to make the batches appear to be IID normally distributed.

Bias Cases

\[\text{bias} = E\left[S^2/n\right] - var(\bar{X}) = \dfrac{-2 \gamma_0 R}{n - 1}\]

• Since \(\gamma_0 > 0\) and \(n > 1\) the sign of the bias depends on the quantity R and thus on the correlation.

• For positive correlation, \(0 \leq \rho_k \leq 1\), the bias will be negative, (\(- \gamma_0 \leq Bias \leq 0\)).

  • In the case of positive correlation, \(S^2/n\) underestimates the \(var(\bar{X})\).
  • Decisions based on positively correlated data will have a higher than planned risk of making an error based on the confidence interval.

• In the case of negatively correlated data, \(S^2/n\) over estimates the \(var(\bar{X})\).

  • A confidence interval based on \(S^2/n\) will be too wide and the true quality of the estimate will be better than indicated.
  • The true confidence coefficient will not be the desired \(1 - \alpha\); it will be greater than \(1 - \alpha\).

Defining the Batch Means

• The method of batch means method divides the data into sub-sequences of contiguous batches: \[\begin{gathered} \underbrace{X_1, X_2, \ldots, X_b}_{batch 1} \cdots \underbrace{X_{b+1}, X_{b+2}, \ldots, X_{2b}}_{batch 2} \cdots \\ \underbrace{X_{(j-1)b+1}, X_{(j-1)b+2}, \ldots, X_{jb}}_{batch j} \cdots \underbrace{X_{(k-1)b+1}, X_{(k-1)b+2}, \ldots, X_{kb}}_{batch k}\end{gathered}\]

• Let \(k\) be the number of batches each consisting of \(b\) observations, so that \(k = \lfloor n/b \rfloor\). If \(b\) is not a divisor of \(n\) then the last \((n - kb)\) data points will not be used.

• Define \(\bar{X}_j(b)\) as the \(j^{th}\) batch mean for \(j = 1, 2, \ldots, k\), where, \[\bar{X}_j(b) = \dfrac{1}{b} \sum_{i=1}^b X_{(j-1)b+i}\]

• Each of the batch means are treated like observations in the batch means series.

Batch Means Confidence Intervals

• Re-label the batch means as \(Y_j = \bar{X}_j(b)\). This produces series of data, (\(Y_1, Y_2, Y_3, \ldots, Y_k\))
• To form a \(1 - \alpha\)% confidence interval, compute the sample average and sample variance of the batch means series:

\[\bar{Y}(k) = \dfrac{1}{k} \sum_{j=1}^k Y_j\]

\[S_b^2 (k) = \dfrac{1}{k - 1} \sum_{j=1}^k (Y_j - \bar{Y}^2)\]

\[\bar{Y}(k) \pm t_{\alpha/2, k-1} \dfrac{S_b (k)}{\sqrt{k}}\]

Using the KSL to Capture Batch Means Statistics

• Add an instance of the StatisticalBatchingElement class to the model.

fun main(){
    val model = Model("Drive Through Pharmacy")
    // add DriveThroughPharmacy to the main model
    val dtp = DriveThroughPharmacyWithQ(model, 1)
    dtp.arrivalRV.initialRandomSource = ExponentialRV(1.0, 1)
    dtp.serviceRV.initialRandomSource = ExponentialRV(0.7, 2)
    model.numberOfReplications = 1
    model.lengthOfReplication = 1000000.0
    model.lengthOfReplicationWarmUp = 100000.0
    val batchingElement = model.statisticalBatching()
    model.simulate()
    val sr = batchingElement.statisticReporter
    println(sr.halfWidthSummaryReport())
}

Summary of Concepts

• Infinite horizon simulation situations attempt to model long term performance.
• Statistical issues related to within replication data characteristics necessitate specialized output analysis method if only 1 simulation run is used.
• The Welch plot is one method to assess the warm up period for an infinite horizon simulation.
• Replication-deletion is the recommended approach to analyzing infinite horizon simulation runs.
• The batch means method should be considered if using only 1 simulation run when performing an infinite horizon simulation analysis.
• Use the StatisticalBatchingElement class to implement the batch means method in the KSL.
• The KSL database functionality will also store the batching statistics.