Chapter 16 Simulation Algorithms for Regenerative Processes

doi:10.1016/S0927-0507(06)13016-9

Handbooks in Operations Research and Management Science

Volume 13, 2006, Pages 477-500

https://doi.org/10.1016/S0927-0507(06)13016-9 Get rights and content

Abstract

This chapter is concerned with reviewing the basic ideas and concepts underlying the use of regenerative structure in the development of efficient simulation algorithms. While it has long been known that discrete state space Markov chains exhibit regenerative structure, we argue that well-behaved discrete-event simulations typically also contain an embedded sequence of regeneration times. However, algorithmic identification of the corresponding regeneration times turns out to be a nontrivial problem. We discuss the theoretical and implementation issues involved in identifying the corresponding regeneration times, and describe how regenerative methodology supplies effective solutions to several difficult simulation problems. In particular, we discuss the use of regeneration in the context of steady-state simulation as a means of efficiently computing confidence intervals and correcting for initial bias. We also point out that regeneration has the potential to offer significant algorithmic efficiency improvements to the simulationist, and illustrate this idea via discussion of steady-state gradient estimation and computation of infinite horizon expected discounted reward.

Introduction

Let $V = (V (t) : t ⩾ 0)$ be a real-valued stochastic process in which $V (t)$ represents the simulation output collected at time t. Roughly speaking, the process V is said to be (classically) regenerative if there exist random times $T (0) < T (1) < \dots$ at which the process “regenerates” (in the sense that V probabilistically starts afresh at each time $T (i)$ , $i ⩾ 0$ , and evolves independently of the process prior to time $T (i)$ ). Such regenerative structure implies that V can be viewed conceptually as a sequence of independent “cycles” $(V (s) : T (i - 1) ⩽ s < T (i))$ that are “pasted together” (where we adopt the convention that $T (- 1) = 0$ ). Thus, the infinite time behavior of V over $[0, \infty)$ is implicitly captured in the behavior of V over a cycle. Hence, in principle, virtually any expectation of V over $[0, \infty)$ can be alternatively described as an expectation involving cycle-related quantities. This observation is the key insight that underlies regenerative simulation.

The use of regenerative structure as an algorithmic tool in the simulation setting has primarily focused on its use in the context of steady-state simulation. The first suggestion that regenerative cycles could play a useful role in steady-state simulation output analysis came from Cox and Smith (1961), and the idea was further developed in Kabak (1968). However, the first comprehensive development of the regenerative method for steady-state simulation output analysis came in a series of papers of Crane and Iglehart, 1974a, Crane and Iglehart, 1974b, Crane and Iglehart, 1975, as well as concurrent work by Fishman, 1973, Fishman, 1974. The great majority of subsequent work on algorithmic exploitation of regeneration has followed the historic tradition of focusing on its application to steady-state simulation output analysis.

In this chapter we focus our discussion on the key theoretical and algorithmic issues underlying the use of regeneration in the steady-state simulation context. We start, in Section 2, by describing the key challenges that confront a simulationist in analyzing steady-state simulation output, while Section 3 discusses the basic regenerative approach to forming an estimator for the so-called “time-average variance constant”. Section 4 offers some discussion of how the particular choice of regenerative structure influences the efficiency of the method, and Section 5 describes the regenerative solution to the initial transient problem and the construction of low-bias steady-state estimators. In Sections 6 When is a simulation regenerative?, 7 When is a GSMP regenerative?, 8 Algorithmic identification of regenerative structure we discuss the theoretical issue of when a simulation is regenerative, with a particular focus on when a discrete-event simulation contains algorithmically identifiable regenerative structure. Section 9 then discusses steady-state regenerative analysis from the perspective of martingale theory.

The last two sections of the chapter are intended to give the reader a hint of the role that regeneration can play in the development of computationally efficient algorithms for other simulation problems. In particular, we show that in computing either steady-state gradients or infinite-horizon expected discounted reward that regeneration offers the simulationist the opportunity to not only construct asymptotically valid confidence statements but to also improve computational efficiency. While regeneration is primarily understood within the simulation community as offering a vehicle for analysis of simulation output, our two examples are intended to argue that regeneration has the potential to also play a significant role in the variance reduction context.

Section snippets

The steady-state simulation problem

Let $V = (V (t) : t ⩾ 0)$ be a real-valued stochastic process in which $V (t)$ represents the value of the simulation output process at (simulated) time t. For example, $V (t)$ could represent the total work-in-process at time t in a production context or the inventory position at time t in a supply chain setting. Throughout this chapter, we use a continuous time formulation to describe the relevant theory. (Note that any discrete-time sequence $(V_{n} : n ⩾ 0)$ can be embedded into continuous time via the definition $V$

The regenerative estimator for the TAVC

To obtain a TAVC estimator that converges to $σ^{2}$ at rate $t^{- 1 / 2}$ , one needs to assume additional structure about the process V. To illustrate this idea, suppose that the simulation output process V is a (continuous-time) autoregressive process satisfying $d V (s) = - γ V (s) d s + d W (s),$ where $γ > 0$ and $W = (W (s) : s ⩾ 0)$ is a square integrable process with stationary independent increments for which $E W (s) = μ s$ and $var W (s) = η^{2} s$ for $s ⩾ 0$ . It is easily verified that $V (t) = e^{- γ t} V (0) + \int_{0}^{t} e^{- γ (t - s)} d W (s)$ and that V satisfies (2)

Choice of the optimal regeneration state

Given a simulation of V over the time interval $[0, t]$ , the natural point estimator for the steady-state mean α is, of course, the time-average $\bar{V} (t)$ . While it may be desirable to modify $\bar{V} (t)$ to deal with initial transient or initial bias effects, one would expect such modifications to be of small order asymptotically. Hence, any reasonable point estimator for α will either be exactly equal to $\bar{V} (t)$ or asymptotically equivalent to $\bar{V} (t)$ . Of course, the r.v. $\bar{V} (t)$ is not influenced in any way by

The regenerative approach to the initial transient and initial bias problems

As discussed in Section 2, one of the major challenges in steady-state simulation is the mitigation of effects due to the initial transient and initial bias. We deal first with the better understood issue of how to reduce biasing effects due to a nonstationary initialization.

It is usual, in the presence of (1), that there exists $ν > 0$ such that $E V (t) = α + O (e^{- ν t})$ as $t \to \infty$ , where $O (h (t))$ represents a function that is bounded by a constant multiple of $| h (t) |$ as $t \to \infty$ . For example, (13) is known to

When is a simulation regenerative?

As has been seen in preceding sections, regenerative structure turns out to be algorithmically useful in developing solutions to various aspects of the steady-state simulation problem. Furthermore, regenerative structure can be easily identified in the setting of discrete state space Markov chains, in either discrete or continuous time.

Of course, most real-world discrete-event simulations do not involve simulating a discrete state space Markov processes. Much more complicated models are

When is a GSMP regenerative?

Section 6 makes clear that regeneration is the rule rather than the exception for well-behaved steady-state simulations. This, however, leaves open the question of when a specific simulation model has the structure necessary to guarantee that the associated steady-state simulation is well behaved.

We shall focus exclusively, in this section, on conditions under which discrete-event simulations possess the required structure. We take the point of view here that a discrete-event simulation is

Algorithmic identification of regenerative structure

Our discussion of Sections 6 When is a simulation regenerative?, 7 When is a GSMP regenerative? makes clear that regenerative structure exists within the typical discrete-event steady-sate simulation. On the other hand, the TAVC estimator of Section 3, as well as the low bias estimators of Section 5, all depend upon the ability of the simulationist to identify the associated regeneration times. Of course, this identification is trivial in the setting of discrete state space Markov chains, where

A martingale perspective on regeneration

To illustrate the connection between martingales and regeneration, we focus here on the case in which $V (t) = g (X (t))$ , where $X = (X (t) : t ⩾ 0)$ is an irreducible finite state continuous time Markov chain with rate matrix $Q = (Q (x, y) : x, y \in S)$ . Given the performance measure $g : S \to R$ (where we choose to encode $g = (g (x) : x \in S)$ as a column vector), the linear system $Q h = - (g - α e)$ has a solution h. Here, e is the column vector in which all entries equal 1, and (19) is called Poisson's equation.

It is a standard fact in the

Efficiency improvement via regeneration: Computing steady state gradients

In many applications settings, it is of interest to compute the sensitivity of the system's performance to perturbations in an underlying parameter. For example, it may be that the arrival rate to a queue is only approximately known, so that computing the change in performance that corresponds to changing the arrival rate is relevant. In particular, computing the derivative (or, more generally, the gradient) of a steady-state performance measure with respect to the arrival rate is a

Efficiency improvement via regeneration: Computing infinite horizon discounted reward

We now offer a second illustration of the principle that the presence of regenerative structure can be usefully exploited to obtain efficiency improvements. Consider the infinite horizon expected discounted reward $α = E D$ , where $D = \int_{0}^{\infty} e^{- r t} V (t) d t$ for some $r > 0$ . From a computational standpoint, an algorithm based on simulating i.i.d. copies of the r.v. D cannot be operationalized, because it takes infinite time to generate the above r.v. As a consequence, one needs to consider computationally feasible

References (36)

P.W. Glynn et al.
Limit theorems for cumulative processes
Stochastic Processes and Their Applications
(1993)
P.W. Glynn et al.
Necessary conditions in limit theorems for cumulative processes
Stochastic Processes and Their Applications
(2002)
T.W. Anderson
The Statistical Analysis of Time Series
(1971)
S. Andradóttir et al.
Accelerated regeneration for Markov chain simulations
Probability in the Engineering and Informational Sciences
(1995)
P. Billingsley
Probability and Measure
(1995)
L. Breiman
Probability
(1968)
J. Calvin
Return state independent quantities in regenerative simulation
Operations Research
(1994)
Calvin, J.M., Glynn, P.W., Nakayama, M.K. (2006). The semi-regenerative method of simulation output analysis. ACM...
D. Cox et al.
Queues
(1961)
H. Cramér
Mathematical Methods of Statistics
(1946)

M. Crane et al.

Simulating stable stochastic systems, I: General multi-server queues

Journal of the Association for Computing Machinery

(1974)

M. Crane et al.

Simulating stable stochastic systems, II: Markov chains

Journal of the Association for Computing Machinery

(1974)

M. Crane et al.

Simulating stable stochastic systems, III: Regenerative processes and discrete-event simulations

Operations Research

(1975)

B. Efron et al.

An Introduction to the Bootstrap

(1993)

S.N. Ethier et al.

Markov Processes: Characterization and Convergence

(1986)

G. Fishman

Statistical analysis for queueing simulations

Management Science

(1973)

G. Fishman

Estimation in multiserver queueing simulations

Operations Research

(1974)

P.W. Glynn

Likelihood ratio gradient estimation for stochastic systems

Communications of the ACM

(1990)

Cited by (16)

On the numerical solution of functional equations with application to response time distributions
2024, Applied Mathematics and Computation
A unified approach is developed to solve functional equations defining generating functions. Such equations are often constructed as a means to solve recurrence relations, such as those arising from queue length probabilities and response time probability distributions in Markov models. Many such equations have been obtained over several decades. Some have been solved analytically, some numerically and for some no tractable solution has been found. Our unified approach is able to provide accurate numerical solutions to such equations, even when they include derivatives of the generating function. It solves the JSQ model with two queues for the first time, utilizing a novel “partial” generating function related to the one the functional equation defines. Numerical results, displayed in tables of moments and graphs of probability density functions, show good accuracy against simulations with 500 000 regenerative cycles.
A semi-product-form for a pair of queues with finite batches: Equilibrium state probabilities and response time densities
2020, Performance Evaluation
A Markovian network of two queues, with finite size batch Poisson arrivals and departures, is solved approximately, but to arbitrary accuracy, for its equilibrium state probabilities. Below a pair of thresholds on the queue lengths, a modification of the Spectral Expansion Method is used to construct a semi-product-form at all lengths of one queue in a finite lattice strip defined by the threshold of the other queue. No additional special arrival streams are required, for example at empty queues, from which it is already known that a product-form can be constructed. Hence the first exact closed form solution for the equilibrium probabilities in an unmodified Markovian queueing network with batches is obtained, the only constraint being finiteness of the batches. The method is illustrated numerically, first in a tandem network and then in a two-node network with feedback. Simulation results confirm convincingly the precision of the method, and partial batch forwarding and discarding are thereby compared quantitatively. Response time distributions are derived using the generating function method, which complements well the semi-product form of the equilibrium state probabilities. Again, agreement with simulation is excellent. The regenerative simulation method was used, so that no warm-up period was needed, and the statistical estimates for the 95% confidence bands are explained for the different cases of state occupancy and response time prediction at equilibrium.
Coupling from the past with randomized quasi-Monte Carlo
2010, Mathematics and Computers in Simulation
The coupling-from-the-past (CFTP) algorithm of Propp and Wilson permits one to sample exactly from the stationary distribution of an ergodic Markov chain. By using it n times independently, we obtain an independent sample from that distribution. A more representative sample can be obtained by creating negative dependence between these n replicates; other authors have already proposed to do this via antithetic variates, Latin hypercube sampling, and randomized quasi-Monte Carlo (RQMC). We study a new, often more effective, way of combining CFTP with RQMC, based on the array-RQMC algorithm. We provide numerical illustrations for Markov chains with both finite and continuous state spaces, and compare with the RQMC combinations proposed earlier.
Automatic regenerative simulation via non-reversible simulated tempering
2023, arXiv
Efficient Input Uncertainty Quantification for Regenerative Simulation
2023, Proceedings - Winter Simulation Conference
Response Time Distribution in a Tandem Pair of Queues with Batch Processing
2021, Journal of the ACM

View all citing articles on Scopus

View full text

Chapter 16 Simulation Algorithms for Regenerative Processes

Abstract

Introduction

Section snippets

The steady-state simulation problem

The regenerative estimator for the TAVC

Choice of the optimal regeneration state

The regenerative approach to the initial transient and initial bias problems

When is a simulation regenerative?

When is a GSMP regenerative?

Algorithmic identification of regenerative structure

A martingale perspective on regeneration

Efficiency improvement via regeneration: Computing steady state gradients

Efficiency improvement via regeneration: Computing infinite horizon discounted reward

Stochastic Processes and Their Applications

Stochastic Processes and Their Applications

The Statistical Analysis of Time Series

Accelerated regeneration for Markov chain simulations

Probability in the Engineering and Informational Sciences

Probability and Measure

Probability

Return state independent quantities in regenerative simulation

Operations Research

Queues

Mathematical Methods of Statistics

Simulating stable stochastic systems, I: General multi-server queues

Journal of the Association for Computing Machinery

Simulating stable stochastic systems, II: Markov chains

Journal of the Association for Computing Machinery

Simulating stable stochastic systems, III: Regenerative processes and discrete-event simulations

Operations Research

An Introduction to the Bootstrap

Markov Processes: Characterization and Convergence

Statistical analysis for queueing simulations

Management Science

Estimation in multiserver queueing simulations

Operations Research

Likelihood ratio gradient estimation for stochastic systems

Communications of the ACM