- Email: [email protected]
Infinitesimal perturbation analysis of a birth and death process
Volume 7, Number I
OPERATIONS RESEARCH LETTERS
February 1988
INFINITESIMAL PERTURBATION ANALYSIS OF A BIRTH AND D E A T H P R O C E S S *
Paul GLASSERMAN Division of Applied Sciences, Harvard University, (79 Pierce Hall, Cambridge, MA.02138, USA Received September 1987 Revised December 1987
Using a birth and death process as an illustrative example, we introduce the notion of alternative representations of stochastic processes and discuss its importance for infinitesimal perturbation analysis derivative estimation. Through a different choice of representation, we are led to an IPA algorithm for a birth and death process better than one discussed by other authors. simulation, birth-death queues * networks of queues
1. Introduction
Infinitesimal perturbation analysis (IPA) is a method of obtaining estimates of derivatives of performance with respect to parameters from a single sample path or simulation of a discrete-event system (Ho and Cao [5], Suri [8]). The consistency of IPA estimates generally depends on the vah'dity of assuming that small parameter changes will effect small changes in state holding times along a sample path but no changes in the sequence Of states. In fact, changing a parameter,'*say 8, to 0 + AO often will cause some change in the state sequence of many sample paths. However, if the probability of a sequence change is 0(210) and the expected change in performance given such a change is also O(AO) then the combined effect is o(AO) and is therefore negligible so far as the derivative of performance is concerned (Cao [1], Zazanis and Suri [9]). In deriving an IPA algorithm, we view a stochastic p~'ocess as one of a parametric family of
* This work was supported in part by Office of Naval Research Contracts N00014-84-K-0465 and N00014-85-K-0075, by National Science Foundation Grants ECS-85-15449 and CDR-85-001-08 and by Army Contract DAAL 03-86-K0171.
processes defined on a common probability space. From the point of view of. applications, the most natural way to achieve this is by taking a simulation algorithm as a construction of a stochastic process. Frequently, the underlying random component in an idealized simulation is a sequence of uniform random variables - the output of the random number generator - and this can be used to construct an ~ underlying probability space (fL F, P) on which are defined all processes to be simulated. Suppose that a family of processes ~ is parameterized by 0; then the idea of IPA is to view a sample path ~(0, ¢0) as a function of 0 for fixed co. Sample derivatives obtained by considering changes in 0 are then used to estimate derivatives of ensemble averages of performance indices for ~. • In general, there may be different ways to effect the construction of the processes on a common probability space, even once the space itself is fixed. Just as there are different ways to simulate random variables, there are different ways to simulate sample paths and each choice of simulation implicitly defines a construction of a process. On a single probability space there may be different functions ~ and ~ such that for each 0, ~(0, .) and ~(0, .) are the same process (in the sense that they induce the same measure on a space of sample paths) though ~(., co) and ~(., ~o) are
0167-6377/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
43
Volume 7, Number I
OPERATIONS RESEARCH LETrERS
different functions of 0 for all to. In analogy with the terminology in Olynn [2], we refer to 4~and as alternative representations of a family of processes. Notice that if ~(00, to) ffi ~(0o, to) but 4~(0o÷ e, to) ~ ( 0 o - e, to) for e ¢ 0, then sample obtained from the same sample path under these two different representations ~ typicany be different. We ~ suggest that, for IPA, a queueing representation - whenever possible - is preferable to an abstract Markov chain representation. The preceding discussion should not obscure the important fact that in IPA no actual parameter change is ever introduced: derivatives evaluated at 0 are estimated from a process simulated at 0 only. Thus, the representation of a family of processes is only important in the derivation of an IPA algorithm. The kind of construction described above could be viewed merely as a 'thought experiment' by which one formulates an IPA estimator. Heidelberger, Cao, Zazanis and Suri [4] develop an 'IPA-like' algorithm for a finite birth and death process with state-dependent rates. Specifically the process has state space {0, 1,...,N} transitions from n to n + 1 occur at rate ?,n (n < N ) and transitions from n to n - 1 occur at rate/t= (n > 0). The IPA-like algorithm estimates d~o/d~ k, where ~r0 is the stationary probability of state 0 and 1 < k < N. Heidelberger et al. show that this algorithm yields estimates which converge to the wrong value and actually predict the wrong sign for the derivative. This is included as an example of the danger in assuming that parameter changes do not lead to changes in the state sequence. In this paper, we will show how a different IPA algorithm circumvents the problems faced by the one in [4]. No new techniques are introduced; the IPA reasoning behind our algorithm is essentially the same as that in [4]. The main difference is in how we view the generation of sample paths. Specifically, we take the birth and death process as a model of a simple cyclic network, while no such interpretation is used in [4]. This point of view leads to a different representation and to different sample path derivatives. Birth and death processes are analytically tractable and so, as a practical matter, there is no reason to simulate them. However, they serve as an interesting illustration of the importance of the choice of representation. 44
February 1988
2. Two models of the same process
Since the essential difference between the algorithm in [4] and ours is in the choice of representation, we will describe the two different approaches. We ~ not present the detailed construction of the tw~ representations; rather, we present them informally as methods of simulation. From this discussion (and the algorithm to be given below), it is hoped that the interested reader could produce the impficit representation involved in each case. This could be accomplished by applying a general construction in [8]. Were we to write out the representations explicitly, we could show that all the sample derivatives to be used below actually exist (almost everywhere) and are genuine random variables. For our present objective, we adopt these as assumptions. (The existence of sample derivatives for a broad class of systems has also been established in [8].) Heidelberger et al. model the birth and death process as follows: every time the process enters state n, for 1 < n < N, samples X1 and X2 are drawn from exponential distributions with rates ~n and ttn, respectively. If X1 < X2 then the holding time in the state is X1 and the next state is n + 1, otherwise the opposite occurs. At the boundary states, only a single sample need be generated. A regenerative cycle begins in state n - 0 and is terminated upon the first return to that state. Based on this representation, [4] uses the following IPA-fike rule to compute d~o/d~k: let T be the length of a regenerative cycle and To the time spent in 0 during a cycle. At the i th transition from k to k + 1 on the original (nominal) path compute -X~O/~k where X~° is the holding time in k immediately preceding the transition. This is the derivative of the holding time in k (corresponding to the fact that ~k is a scale parameter ([8])). The sum of such terms over all transitions from k to k + 1 in a regenerative cycle is T', the derivative of the length of the cycle. Clearly, ~r0 --E[To]/E[T]. Moreover, under this representation, in almost any cycle, a sufficiently small change in ?~k results in no change in To; only the durations of those visits to k that are followed by visits to k + 1 will change. Thus, the ~ample derivative of ~r0 under this representation is (T~T- TO_T')/T 2 - - T'(TO/T2). If T'(m), To(m) and T(m) denote the averages of T', To
o
VoluL:e 7, Number 1
OPERATIONS RESEARCH LETrERS
and Tover m regenerative cycles, then the IPA-like estimator o b t a i n e d from m cycles is - T ' ( m ) { o ( m ) / T ( m ) 2 ) . Over infinitely many regenerative cycles the IPA-like estimator converges to -E[T'I(E[To]/E[T]2), which is positive since T ' < 0. Direct computation shows that d~ro/d?~k is negative. As is pointed out in [4], this algorithm ignores the fact that if any actual perturbation A~,k were introduced there would be a positive probability that a transition from k to k - 1 on the nominal path would become a transition from k to k + 1 on the perturbed path. This would appear to be why this algorithm yields incorrect estimates. However, as we shall see, the real reason for the failure of this algorithm - of this choice of representation - has to do with what happens given such a change rathe~ than with the possibility of such a change. Consider, now, a cyclic queueing network in which N customers alternate service at two exponential servers with load-dependent rates. Server 1 works at "rate it, when its queue contains n customers and server 2 works at rate ~,~_, when its queue contains n customers. (The customer in service is counted among those in the queue.) Then the number of customers at server 1 is a birth and death process like the one described above and we shall take it as our model of the process. Our regenerative cycles begin and end with the termination of an idle period at server 1: if n is the number of customers at server 1 then the process begins just after a jump from ~ n = 0 to n = 1 and ends the next time such a jump occurs. Thus, each cycle consists of a single busy period and idle period at server 1. Using tKe notation defined above (but interpreted according to the new representation) the estimator obtained from m cycles is ~'(~( m ) = ((To'(m) T ( m ) ~ o ( m ) ~ ' ( m ) ) / T ( m ) 2. It is worth noting that the sample derivatives have different values under this representation than they do in [4]; in particular, it will become clear that TO' is no longer zero.
February 1988
Ho and Yang [7]. Suppose the scheduled service completion time of a customer in service at a load-dependent server with n customers in queue and rate It, is t}"). Then if another customer arrives at time to, before the service completion, the new service completion time t} "+1) can be taken to be (1)
Itn+l It follows that"
d t ( n + l ) ( 1 - I t - - - L ) d - ~ - t + It, d d it k ! = It.+1 d it l, ° It,,+1 d lt--'-kkt } " )
+('I ")- t -,~ d { I t .
)"
(2)
The derivatives of the original scheduled departure time, of the arrival time and of the ratio of rates all contribute to the derivative of the new scheduled completion time. Suppose the i th customer begins service with n customers present, n > 0, and the previous customer completed service with (n + 1) customers present; then d
dItk .f d dpktl")(i)=
~
n
tIn+')(t - 1)
[ -(tln)(i)-tln+')(i - 1 ) ) / / z k, n-k,
(3) where t~')(j) is the service completion time of the j t h customer, if this customer completed service in state m. The term common to both cases reflects the fact that a small delay in the previous service completion time is propagated to the next service completion time, if there is no idle period in between. The addition~ term in the case n = k is the derivative of the sample service time when service is provided at rate It k- (As before, the form of this term is due to the fact that Itk is a scale parameter.) On the other hand, if n - 0 and k > 1, and we take t~°) = 0, (2) reduces to
3. The algorithm
d t},)(i)= dp~
d t.(i),
Before presenting the algorithm we briefly describe the infinitesimal perturbation propagation rules for load-dependent servers, first derived in
where to(i) is thearrival time of the ith customer. That is, a small delay in a customer's arrival time is propagated to a delay in that customer's sched-
(4)
45
Volume 7, Number 1
OPERATIONS RESEARCH LETI'ERS
uled d e p ~ time if the customer finds the server idle. Eq. (2), (3) and (4) describe the perturbation propagation rules for load-dependent
Table 1 d % / d g k for N ffi 4 Parameters
servers.
To calculate ~Y' and To' we apply these rules at ~ transition along a regenerative cycle of the p~, ~ t T - T o - So so that So is the ti_me server 1 becomes idle. Since So and T are service completion times at servers 1 and 2, respectively, T ' and ~ (hence, also To') can be calculated as the p r o o ~ is simulated, using rules (2)-(4). It will be useful to define some additional notation for the algorithm. Let rl(.0ffi/zJ~t~+l, ra(i) - ~ , J ~ i - l , i - 1,..., N - 1. Define dr2(i ) to be 1/~k-1 if i - k , --Xk+l/X2k if i ffik+ 1 and zero otherwise. Denote by Y1 and Y2 generic independent exponential random variables with mean 1. At any time t, t / , - tl,(t ) is the next service completion time at server i. As before, n denotes the number of customers at server 1. With this notation, we have the following algorithm to compute T ' and To' from a regenerative cycle of the system.
Algorithm Step 1. (Initialization) t "ffi0; n "ffi1, generate Y1, tl, "ffi Y1/Izl, gene;~ate Y2, if, "= Y2/~, I, endcycle ,= false.
Step 2.
if t/,
then t ' - t l , , 80 to Step 3, else t ,= t/,, 80 to Step 4.
Step 3. (Departure from server 1) n*== n - - 1;
if (n>0), then generate Y1, t/, "ffi t + Y1/Itn,
else t/, "ffi oo, if (n < N - 1) then dt 2 "-- r(n + 1) × d t 2 + (1 - r2(n + 1)) × dt 1 + dr2(n + 1) × (t/2 - t)(rule (2)), tlz "= t + rz(n + 1) × (t/2 - t), if (n = N - 1) then dt2 :ffi dr1 (rule (4)), generate Y2, tl 2 "= t + Y 2 / ~ s - 1, if ( n = 0 ) then SO "ffi t,
dSo "ffidtl, go to Step 5.
Step4. (Arrival to server 1) n :-- n + 1,
46
February 1988
k
Theoretical
Experimental -0.0407
0
-0.0408
ffi(4.0, 3.0, 2.0,1.0)
1
-0.0357
-0.0357
~, = (3.0, 3.0, 3.0, 3.0)
2 3
-0.0306 -0.0230
-0.0304 -0.0229
0
-0.1126 -0.0447 -0.0288
= (5.0,3.0, 3.0, 5.0)
1
ffi(2.0,4.0, 4.0, 2.0)
2
-0.u24 -0.0446 -0.0290
3 0
-0.0166 -0.2153
-0.0166 -0.2151
1 -0.~554 2 " -0.0290
-0.0542 -0.0288
3
-0.0166
= (2.0,2.0, 2.o, 2.0) = (1.0, 3.0,1.0, 3.0)
-0.0166
if (n = 1) then generate Y1, t/a "ffi t + Y1/!~1, dtl :-- dt 2 (rule (4)), endcycle , - true, if (n > 1) then dtl "ffirl(n - 1) × dt 1 + (1 - rl(n - 1)) × dt 2 (rule (2)), tf, "ffit + r l ( n - 1) × (tn - t), if (n < N ) then generate Y2, if2 "ffit + Y2/Xn, else '~2 "= ¢~, if (n - k) then dt 2 .ffidt2)(tl2- t)/)t 2 (rule (3)).
Step 5. if (not endcycle) then go to Step 2, T:ffi t,
To'ffir-So, dT'ffi dtl, dTo "ffidtl - d S o , STOP. At the end of execution, d T and d T o contain the exact values of the sample derivatives of T
Table 2 d%/d~t k for N - 6 Parameters
k Theoretical Experimental
0 /s - (1.0, 2.0, 3.0, 4.0, 3.0, 2.0) 2 ffi(3.0, 3.0, 3.0, 2.0, 2.0, 2.0) 4 1 /t - (1.0, 2.0,1.0, 2.0,1.0, 2.0) 3 ~, ffi (2.0,1.0, 2.0,1.0, 2.0,1.0) 5 0 /t = (4.0, 3.0, 2.0,1.0, 2.0, 4.0) 2 = (2.0, 2.0, 2.0, 2.0, 2.0, 2.0) 4
-0.0173 -0.0098 -0.OO45 -0.0700 -0.0400 - 0.0100 -0.0964 -0.0681 -0.0340
-0.0170 -0.0093 -O.0O45 -0.0695 -0.0401 -0.0096 -0.0956 -0.0642 -0.0334
Volume7, Number 1
OPERATIONSRESEARCHLETTERS
and To for a regenerative cycle. We saw in Section 2 how these derivatives combine with estimates of T and To (obtained simultaneously) to produce estimates of ~r~. (Essentially the same idea could be used to define estimators of derivatives of the other state probabilities.) Numerical results displayed in Tables 1 and 2 indicate that our estimator converges to the correct value. The displayed results represent averages over 100,000 regenerative cycles.
4. Queueing networks ver~ns Markov chains The algorithm presented here succeeds where another failed because it utilizes a different, more suitable, representation of the processes. The dependence of each path on ~k in the simple network described is different from that in the abstract birth and death process, though the statistical dependence (over the ensemble of paths) is clearly the same. There are usually many ways to implement a simulation. Traditional concerns in the comparison of implementations include, for example, variance reduction and computational efficiency. IPA introduces a new consideration; namely, the properties of sample derivatives computed from a simulation. When using IPA, one should - as always - choose the implementation to advantage. We will now try to give some insight into why the particular representation used here performs well and why, more generally, queueing network representations (whenever possible) will do better than abstract Markov chains. We will make the comparison for the somewhat simpler case of a birth and death process all of whose rates are initially the same, say, / L ~ - ~ _ 1 - 1 . (Note that algorithm in [4] fails in this case as well.) We consider the effect of a small change in k k, 0 < k
February1988
first transition is determined as in AMC but the
next transition would be determined by min(X2, ]'1 - X1) if we assume that Y1 > X1. Let us further suppose that on some nominal path ]'1 - X1 < X2 so that after the transition from k to k - 1 we have one from k - 1 to k. Consider the effect of a small increase in ~kFor both simulations, this means that the transition out of k will now be determined by min(Xl, Y l - &l) for some A~ > 0. Suppose that, in fact, ]'1 - A1 < X1. "l~henwe would have a change in the state sequence: the transition from k to k - 1 on the nominal path becomes a transition from k to k + 1 on the perturbed path for both simulations. But now consider the next transition. In the AMC simulation, we look at min(X2, Y2) and make a transition up or down in exactly the same way as on the nominal path. Thus, if the nominal state sequence had been (k, k - 1,k,...) the perturbed state sequence would be (k, k + 1, k + 2,...). On the' other hand, if we look at the QN simulation we draw a different conclusion. The second transition is determined by min(X~ - Y1 + A1, Y2)- Then, since X t - ]'1 < 0, for A1 small enough, we have min(X1 - ]'1 + A1, Y2)ffi X1 - Yl + A~ so that, for ,~ in some interval, the perturbed state sequence will be (k, k + 1, k , . . . ) where the nominal state sequence was (k, k - 1, k,...). Furthermore, on the perturbed path the next transition will be determined by min(X 2, Y2 - (X1 - ]'1 + AI) - &2)(where A2 is the new perturbation generated in state k) and on the nominal path it is determined by rain(X2- (Y1 XO, ]'2). Thus, not only do the paths converge at k.after a brief separation, their transition mechanisms also stay close together: if A~ + A2 is sufficiently small then the next state will be the same on both paths. As more and more perturbations are accumulated, it becomes more likely that the nominal and perturbed paths will be quite different. Nevertheless, it should be clear from the example above that the QN simulation is much more tolerant of perturbations, tending to keep the nominal and perturbed paths together while the AMC simulation actually tends to keep them apart. The key observation is this: using residual lifetimes rather than generating new lifetimes at eacb transition restricts the possible ch~ges in event sequences for small parameter changes. 47
Volume 7, NumberI
OPERATIONSRESEARCHLETTERS
We have only analyzed a specific case but it is typical of the situation in general. Similar arguments can be made for the case where the transition rates are initially state-dependent. Indeed, similar insight can be obtained for larger networks where difference between the two methods is even more pronounced sin~ the number of possible ~ t i o u s out of any state is greater. If, in a network, a transition is determined by
mm(Xl",...,
(4)
and the next transition is determined by
(
- x,,,) )),
(5)
then, for A in some interval, the only effect of replacing Xa(1) with X~(1) + A is to reverse the order of these two transitions- that is, to change a single state in the state sequence. In light of the discussion above, let us reiterate a point made earlier. It might appear that the IPA-like algorithm in [4] fails because it ignores changes in the state sequence. Our results (as well as those in [1] and [9]) should indicate that the real issue is the effect of changes in the state sequence. Our algorithm also ignores changes in the state sequence but it works because it relies on a representation of the sample paths that limits the 'damage' caused by such changes. For some systems, order changes cannot be ignored regardless of the representation. New perturbation analysis techniques have been proposed to deal with such cases (Gong and Ho [3], Ho anti Li [6]).
February1988
parameter changes in the representation used is irrelevant. From this perspective, an IPA estimator is simply a sample path functional. Once its definition is fixed (using a particular representation) how the sample paths on which it acts are generated is immaterial. A similar point is made in [2] and [8], where it is shown that derivatives of real-valued random variables can often be generated using an expression derived in [8] - based on viewing the random variables as having been generated by the inversion method - even if the random variables are actually generated by a different method. In our birth and death example, the IPA algorithm obtained using the queueing network representation could be used on a simulation that follows the representation in [4]: given a sample path generated by any means we could apply our IPA estimator by viewing the path as though it had been generatedusing our representation. It is not hard to see from the algorithm in Section 3 that the sample derivatives calculated from a regenerative cycle depend only on the state sequence and the sojourn times. Provided these have the correct joint distribution, IPA estimates obtained from them are independent of how they were generated. t
Acknowledgments I am grateful to P~'~fessor Y.C. Ho, W.B. Gong, P. Vakili and Professor R. Suri for helpful discussions; and to an anonymous Associate Editor for a careful reading of the manuscript.
$. A note on implementation Thus far, we have emphasized that different representations lead to different IPA algorithms and different derivative estimators. As a final remark, we point out that although the choice of representation is important in the derivation of the algorithm, for the actual implementation it is usually the case that any simulation based on some valid representation will do - provided the IPA a l g o r i ~ is incorporated into the simulation as though the 'right' representation were being used, The crucial point is ~ a t in IPA no parameter change is ever really ~trodUced; hence, the actual sample path by sample path effect of 48
References
[1] X.R. Cao, "Convergenceof parametersensitivityestimates in a stochastic experiment", IEEE Transactions on Automatic Control AC-30,845-853(1985). [2] P.W. Glynn, "Process-differentiable representations of parametric families of random variables", Technical Report, Department of Industrial Engineering, Univet.~ityof Wisconsin-Madison,Wl, 1987. [3] W.B.Gongand Y.C. Ho, "Smoothedperturbation analysis of discrete-eventdynamicsystems", IEEE Transactions on Automatic Control AC-3Z856-868 (1987). [4] P. Heidelberger,X.R. Cao, M. Zazanis and R. Suri,"Convergence properties of infinitesimalperturbation analysis estimates", subn~tt=d for pubEca~n, 1987(revisedversion of IBM Research Report RCl1891 (53575)).
Volume 7, Number I
OPEI~ATIONS RESEARCH LETrERS
[5] Y.C. Ho and X.R. Cao, "Perturbation analysis and optimization of queueing networks", Journal of Optimization Theory and Applications 40 559-582 (1983). [6] Y.C. Ho and S. Li, "Extensions of perturbation analysis of discrete-event dynamic systems", IEEE Tra~sactions on Automatic Control, forthcoming May, 1988. [7] Y.C. Ho and P.Q. Yang, "Equivalent network, load-dependent servers and perturbation analysis - An experimental
February 1988
study", Proc. Int. Conf. Teletraffic Anal. and Comp. Perf. Ev., North-Holland, New York, 1986. [8] R. Suri, "Infinitesimal perturbation analysis for general discrete event systems", Journal of the Association for Computing Machinery 34 686-717 (1987). [9] M. Zazanis and R. Suri, "Estimating second derivatives of performance measures for G / G / 1 queues from a single sample path", submitted for publication, 1985.
49
OPERATIONS RESEARCH LETTERS
February 1988
INFINITESIMAL PERTURBATION ANALYSIS OF A BIRTH AND D E A T H P R O C E S S *
Paul GLASSERMAN Division of Applied Sciences, Harvard University, (79 Pierce Hall, Cambridge, MA.02138, USA Received September 1987 Revised December 1987
Using a birth and death process as an illustrative example, we introduce the notion of alternative representations of stochastic processes and discuss its importance for infinitesimal perturbation analysis derivative estimation. Through a different choice of representation, we are led to an IPA algorithm for a birth and death process better than one discussed by other authors. simulation, birth-death queues * networks of queues
1. Introduction
Infinitesimal perturbation analysis (IPA) is a method of obtaining estimates of derivatives of performance with respect to parameters from a single sample path or simulation of a discrete-event system (Ho and Cao [5], Suri [8]). The consistency of IPA estimates generally depends on the vah'dity of assuming that small parameter changes will effect small changes in state holding times along a sample path but no changes in the sequence Of states. In fact, changing a parameter,'*say 8, to 0 + AO often will cause some change in the state sequence of many sample paths. However, if the probability of a sequence change is 0(210) and the expected change in performance given such a change is also O(AO) then the combined effect is o(AO) and is therefore negligible so far as the derivative of performance is concerned (Cao [1], Zazanis and Suri [9]). In deriving an IPA algorithm, we view a stochastic p~'ocess as one of a parametric family of
* This work was supported in part by Office of Naval Research Contracts N00014-84-K-0465 and N00014-85-K-0075, by National Science Foundation Grants ECS-85-15449 and CDR-85-001-08 and by Army Contract DAAL 03-86-K0171.
processes defined on a common probability space. From the point of view of. applications, the most natural way to achieve this is by taking a simulation algorithm as a construction of a stochastic process. Frequently, the underlying random component in an idealized simulation is a sequence of uniform random variables - the output of the random number generator - and this can be used to construct an ~ underlying probability space (fL F, P) on which are defined all processes to be simulated. Suppose that a family of processes ~ is parameterized by 0; then the idea of IPA is to view a sample path ~(0, ¢0) as a function of 0 for fixed co. Sample derivatives obtained by considering changes in 0 are then used to estimate derivatives of ensemble averages of performance indices for ~. • In general, there may be different ways to effect the construction of the processes on a common probability space, even once the space itself is fixed. Just as there are different ways to simulate random variables, there are different ways to simulate sample paths and each choice of simulation implicitly defines a construction of a process. On a single probability space there may be different functions ~ and ~ such that for each 0, ~(0, .) and ~(0, .) are the same process (in the sense that they induce the same measure on a space of sample paths) though ~(., co) and ~(., ~o) are
0167-6377/88/$3.50 © 1988, Elsevier Science Publishers B.V. (North-Holland)
43
Volume 7, Number I
OPERATIONS RESEARCH LETrERS
different functions of 0 for all to. In analogy with the terminology in Olynn [2], we refer to 4~and as alternative representations of a family of processes. Notice that if ~(00, to) ffi ~(0o, to) but 4~(0o÷ e, to) ~ ( 0 o - e, to) for e ¢ 0, then sample obtained from the same sample path under these two different representations ~ typicany be different. We ~ suggest that, for IPA, a queueing representation - whenever possible - is preferable to an abstract Markov chain representation. The preceding discussion should not obscure the important fact that in IPA no actual parameter change is ever introduced: derivatives evaluated at 0 are estimated from a process simulated at 0 only. Thus, the representation of a family of processes is only important in the derivation of an IPA algorithm. The kind of construction described above could be viewed merely as a 'thought experiment' by which one formulates an IPA estimator. Heidelberger, Cao, Zazanis and Suri [4] develop an 'IPA-like' algorithm for a finite birth and death process with state-dependent rates. Specifically the process has state space {0, 1,...,N} transitions from n to n + 1 occur at rate ?,n (n < N ) and transitions from n to n - 1 occur at rate/t= (n > 0). The IPA-like algorithm estimates d~o/d~ k, where ~r0 is the stationary probability of state 0 and 1 < k < N. Heidelberger et al. show that this algorithm yields estimates which converge to the wrong value and actually predict the wrong sign for the derivative. This is included as an example of the danger in assuming that parameter changes do not lead to changes in the state sequence. In this paper, we will show how a different IPA algorithm circumvents the problems faced by the one in [4]. No new techniques are introduced; the IPA reasoning behind our algorithm is essentially the same as that in [4]. The main difference is in how we view the generation of sample paths. Specifically, we take the birth and death process as a model of a simple cyclic network, while no such interpretation is used in [4]. This point of view leads to a different representation and to different sample path derivatives. Birth and death processes are analytically tractable and so, as a practical matter, there is no reason to simulate them. However, they serve as an interesting illustration of the importance of the choice of representation. 44
February 1988
2. Two models of the same process
Since the essential difference between the algorithm in [4] and ours is in the choice of representation, we will describe the two different approaches. We ~ not present the detailed construction of the tw~ representations; rather, we present them informally as methods of simulation. From this discussion (and the algorithm to be given below), it is hoped that the interested reader could produce the impficit representation involved in each case. This could be accomplished by applying a general construction in [8]. Were we to write out the representations explicitly, we could show that all the sample derivatives to be used below actually exist (almost everywhere) and are genuine random variables. For our present objective, we adopt these as assumptions. (The existence of sample derivatives for a broad class of systems has also been established in [8].) Heidelberger et al. model the birth and death process as follows: every time the process enters state n, for 1 < n < N, samples X1 and X2 are drawn from exponential distributions with rates ~n and ttn, respectively. If X1 < X2 then the holding time in the state is X1 and the next state is n + 1, otherwise the opposite occurs. At the boundary states, only a single sample need be generated. A regenerative cycle begins in state n - 0 and is terminated upon the first return to that state. Based on this representation, [4] uses the following IPA-fike rule to compute d~o/d~k: let T be the length of a regenerative cycle and To the time spent in 0 during a cycle. At the i th transition from k to k + 1 on the original (nominal) path compute -X~O/~k where X~° is the holding time in k immediately preceding the transition. This is the derivative of the holding time in k (corresponding to the fact that ~k is a scale parameter ([8])). The sum of such terms over all transitions from k to k + 1 in a regenerative cycle is T', the derivative of the length of the cycle. Clearly, ~r0 --E[To]/E[T]. Moreover, under this representation, in almost any cycle, a sufficiently small change in ?~k results in no change in To; only the durations of those visits to k that are followed by visits to k + 1 will change. Thus, the ~ample derivative of ~r0 under this representation is (T~T- TO_T')/T 2 - - T'(TO/T2). If T'(m), To(m) and T(m) denote the averages of T', To
o
VoluL:e 7, Number 1
OPERATIONS RESEARCH LETrERS
and Tover m regenerative cycles, then the IPA-like estimator o b t a i n e d from m cycles is - T ' ( m ) { o ( m ) / T ( m ) 2 ) . Over infinitely many regenerative cycles the IPA-like estimator converges to -E[T'I(E[To]/E[T]2), which is positive since T ' < 0. Direct computation shows that d~ro/d?~k is negative. As is pointed out in [4], this algorithm ignores the fact that if any actual perturbation A~,k were introduced there would be a positive probability that a transition from k to k - 1 on the nominal path would become a transition from k to k + 1 on the perturbed path. This would appear to be why this algorithm yields incorrect estimates. However, as we shall see, the real reason for the failure of this algorithm - of this choice of representation - has to do with what happens given such a change rathe~ than with the possibility of such a change. Consider, now, a cyclic queueing network in which N customers alternate service at two exponential servers with load-dependent rates. Server 1 works at "rate it, when its queue contains n customers and server 2 works at rate ~,~_, when its queue contains n customers. (The customer in service is counted among those in the queue.) Then the number of customers at server 1 is a birth and death process like the one described above and we shall take it as our model of the process. Our regenerative cycles begin and end with the termination of an idle period at server 1: if n is the number of customers at server 1 then the process begins just after a jump from ~ n = 0 to n = 1 and ends the next time such a jump occurs. Thus, each cycle consists of a single busy period and idle period at server 1. Using tKe notation defined above (but interpreted according to the new representation) the estimator obtained from m cycles is ~'(~( m ) = ((To'(m) T ( m ) ~ o ( m ) ~ ' ( m ) ) / T ( m ) 2. It is worth noting that the sample derivatives have different values under this representation than they do in [4]; in particular, it will become clear that TO' is no longer zero.
February 1988
Ho and Yang [7]. Suppose the scheduled service completion time of a customer in service at a load-dependent server with n customers in queue and rate It, is t}"). Then if another customer arrives at time to, before the service completion, the new service completion time t} "+1) can be taken to be (1)
Itn+l It follows that"
d t ( n + l ) ( 1 - I t - - - L ) d - ~ - t + It, d d it k ! = It.+1 d it l, ° It,,+1 d lt--'-kkt } " )
+('I ")- t -,~ d { I t .
)"
(2)
The derivatives of the original scheduled departure time, of the arrival time and of the ratio of rates all contribute to the derivative of the new scheduled completion time. Suppose the i th customer begins service with n customers present, n > 0, and the previous customer completed service with (n + 1) customers present; then d
dItk .f d dpktl")(i)=
~
n
tIn+')(t - 1)
[ -(tln)(i)-tln+')(i - 1 ) ) / / z k, n-k,
(3) where t~')(j) is the service completion time of the j t h customer, if this customer completed service in state m. The term common to both cases reflects the fact that a small delay in the previous service completion time is propagated to the next service completion time, if there is no idle period in between. The addition~ term in the case n = k is the derivative of the sample service time when service is provided at rate It k- (As before, the form of this term is due to the fact that Itk is a scale parameter.) On the other hand, if n - 0 and k > 1, and we take t~°) = 0, (2) reduces to
3. The algorithm
d t},)(i)= dp~
d t.(i),
Before presenting the algorithm we briefly describe the infinitesimal perturbation propagation rules for load-dependent servers, first derived in
where to(i) is thearrival time of the ith customer. That is, a small delay in a customer's arrival time is propagated to a delay in that customer's sched-
(4)
45
Volume 7, Number 1
OPERATIONS RESEARCH LETI'ERS
uled d e p ~ time if the customer finds the server idle. Eq. (2), (3) and (4) describe the perturbation propagation rules for load-dependent
Table 1 d % / d g k for N ffi 4 Parameters
servers.
To calculate ~Y' and To' we apply these rules at ~ transition along a regenerative cycle of the p~, ~ t T - T o - So so that So is the ti_me server 1 becomes idle. Since So and T are service completion times at servers 1 and 2, respectively, T ' and ~ (hence, also To') can be calculated as the p r o o ~ is simulated, using rules (2)-(4). It will be useful to define some additional notation for the algorithm. Let rl(.0ffi/zJ~t~+l, ra(i) - ~ , J ~ i - l , i - 1,..., N - 1. Define dr2(i ) to be 1/~k-1 if i - k , --Xk+l/X2k if i ffik+ 1 and zero otherwise. Denote by Y1 and Y2 generic independent exponential random variables with mean 1. At any time t, t / , - tl,(t ) is the next service completion time at server i. As before, n denotes the number of customers at server 1. With this notation, we have the following algorithm to compute T ' and To' from a regenerative cycle of the system.
Algorithm Step 1. (Initialization) t "ffi0; n "ffi1, generate Y1, tl, "ffi Y1/Izl, gene;~ate Y2, if, "= Y2/~, I, endcycle ,= false.
Step 2.
if t/,
then t ' - t l , , 80 to Step 3, else t ,= t/,, 80 to Step 4.
Step 3. (Departure from server 1) n*== n - - 1;
if (n>0), then generate Y1, t/, "ffi t + Y1/Itn,
else t/, "ffi oo, if (n < N - 1) then dt 2 "-- r(n + 1) × d t 2 + (1 - r2(n + 1)) × dt 1 + dr2(n + 1) × (t/2 - t)(rule (2)), tlz "= t + rz(n + 1) × (t/2 - t), if (n = N - 1) then dt2 :ffi dr1 (rule (4)), generate Y2, tl 2 "= t + Y 2 / ~ s - 1, if ( n = 0 ) then SO "ffi t,
dSo "ffidtl, go to Step 5.
Step4. (Arrival to server 1) n :-- n + 1,
46
February 1988
k
Theoretical
Experimental -0.0407
0
-0.0408
ffi(4.0, 3.0, 2.0,1.0)
1
-0.0357
-0.0357
~, = (3.0, 3.0, 3.0, 3.0)
2 3
-0.0306 -0.0230
-0.0304 -0.0229
0
-0.1126 -0.0447 -0.0288
= (5.0,3.0, 3.0, 5.0)
1
ffi(2.0,4.0, 4.0, 2.0)
2
-0.u24 -0.0446 -0.0290
3 0
-0.0166 -0.2153
-0.0166 -0.2151
1 -0.~554 2 " -0.0290
-0.0542 -0.0288
3
-0.0166
= (2.0,2.0, 2.o, 2.0) = (1.0, 3.0,1.0, 3.0)
-0.0166
if (n = 1) then generate Y1, t/a "ffi t + Y1/!~1, dtl :-- dt 2 (rule (4)), endcycle , - true, if (n > 1) then dtl "ffirl(n - 1) × dt 1 + (1 - rl(n - 1)) × dt 2 (rule (2)), tf, "ffit + r l ( n - 1) × (tn - t), if (n < N ) then generate Y2, if2 "ffit + Y2/Xn, else '~2 "= ¢~, if (n - k) then dt 2 .ffidt2)(tl2- t)/)t 2 (rule (3)).
Step 5. if (not endcycle) then go to Step 2, T:ffi t,
To'ffir-So, dT'ffi dtl, dTo "ffidtl - d S o , STOP. At the end of execution, d T and d T o contain the exact values of the sample derivatives of T
Table 2 d%/d~t k for N - 6 Parameters
k Theoretical Experimental
0 /s - (1.0, 2.0, 3.0, 4.0, 3.0, 2.0) 2 ffi(3.0, 3.0, 3.0, 2.0, 2.0, 2.0) 4 1 /t - (1.0, 2.0,1.0, 2.0,1.0, 2.0) 3 ~, ffi (2.0,1.0, 2.0,1.0, 2.0,1.0) 5 0 /t = (4.0, 3.0, 2.0,1.0, 2.0, 4.0) 2 = (2.0, 2.0, 2.0, 2.0, 2.0, 2.0) 4
-0.0173 -0.0098 -0.OO45 -0.0700 -0.0400 - 0.0100 -0.0964 -0.0681 -0.0340
-0.0170 -0.0093 -O.0O45 -0.0695 -0.0401 -0.0096 -0.0956 -0.0642 -0.0334
Volume7, Number 1
OPERATIONSRESEARCHLETTERS
and To for a regenerative cycle. We saw in Section 2 how these derivatives combine with estimates of T and To (obtained simultaneously) to produce estimates of ~r~. (Essentially the same idea could be used to define estimators of derivatives of the other state probabilities.) Numerical results displayed in Tables 1 and 2 indicate that our estimator converges to the correct value. The displayed results represent averages over 100,000 regenerative cycles.
4. Queueing networks ver~ns Markov chains The algorithm presented here succeeds where another failed because it utilizes a different, more suitable, representation of the processes. The dependence of each path on ~k in the simple network described is different from that in the abstract birth and death process, though the statistical dependence (over the ensemble of paths) is clearly the same. There are usually many ways to implement a simulation. Traditional concerns in the comparison of implementations include, for example, variance reduction and computational efficiency. IPA introduces a new consideration; namely, the properties of sample derivatives computed from a simulation. When using IPA, one should - as always - choose the implementation to advantage. We will now try to give some insight into why the particular representation used here performs well and why, more generally, queueing network representations (whenever possible) will do better than abstract Markov chains. We will make the comparison for the somewhat simpler case of a birth and death process all of whose rates are initially the same, say, / L ~ - ~ _ 1 - 1 . (Note that algorithm in [4] fails in this case as well.) We consider the effect of a small change in k k, 0 < k
February1988
first transition is determined as in AMC but the
next transition would be determined by min(X2, ]'1 - X1) if we assume that Y1 > X1. Let us further suppose that on some nominal path ]'1 - X1 < X2 so that after the transition from k to k - 1 we have one from k - 1 to k. Consider the effect of a small increase in ~kFor both simulations, this means that the transition out of k will now be determined by min(Xl, Y l - &l) for some A~ > 0. Suppose that, in fact, ]'1 - A1 < X1. "l~henwe would have a change in the state sequence: the transition from k to k - 1 on the nominal path becomes a transition from k to k + 1 on the perturbed path for both simulations. But now consider the next transition. In the AMC simulation, we look at min(X2, Y2) and make a transition up or down in exactly the same way as on the nominal path. Thus, if the nominal state sequence had been (k, k - 1,k,...) the perturbed state sequence would be (k, k + 1, k + 2,...). On the' other hand, if we look at the QN simulation we draw a different conclusion. The second transition is determined by min(X~ - Y1 + A1, Y2)- Then, since X t - ]'1 < 0, for A1 small enough, we have min(X1 - ]'1 + A1, Y2)ffi X1 - Yl + A~ so that, for ,~ in some interval, the perturbed state sequence will be (k, k + 1, k , . . . ) where the nominal state sequence was (k, k - 1, k,...). Furthermore, on the perturbed path the next transition will be determined by min(X 2, Y2 - (X1 - ]'1 + AI) - &2)(where A2 is the new perturbation generated in state k) and on the nominal path it is determined by rain(X2- (Y1 XO, ]'2). Thus, not only do the paths converge at k.after a brief separation, their transition mechanisms also stay close together: if A~ + A2 is sufficiently small then the next state will be the same on both paths. As more and more perturbations are accumulated, it becomes more likely that the nominal and perturbed paths will be quite different. Nevertheless, it should be clear from the example above that the QN simulation is much more tolerant of perturbations, tending to keep the nominal and perturbed paths together while the AMC simulation actually tends to keep them apart. The key observation is this: using residual lifetimes rather than generating new lifetimes at eacb transition restricts the possible ch~ges in event sequences for small parameter changes. 47
Volume 7, NumberI
OPERATIONSRESEARCHLETTERS
We have only analyzed a specific case but it is typical of the situation in general. Similar arguments can be made for the case where the transition rates are initially state-dependent. Indeed, similar insight can be obtained for larger networks where difference between the two methods is even more pronounced sin~ the number of possible ~ t i o u s out of any state is greater. If, in a network, a transition is determined by
mm(Xl",...,
(4)
and the next transition is determined by
(
- x,,,) )),
(5)
then, for A in some interval, the only effect of replacing Xa(1) with X~(1) + A is to reverse the order of these two transitions- that is, to change a single state in the state sequence. In light of the discussion above, let us reiterate a point made earlier. It might appear that the IPA-like algorithm in [4] fails because it ignores changes in the state sequence. Our results (as well as those in [1] and [9]) should indicate that the real issue is the effect of changes in the state sequence. Our algorithm also ignores changes in the state sequence but it works because it relies on a representation of the sample paths that limits the 'damage' caused by such changes. For some systems, order changes cannot be ignored regardless of the representation. New perturbation analysis techniques have been proposed to deal with such cases (Gong and Ho [3], Ho anti Li [6]).
February1988
parameter changes in the representation used is irrelevant. From this perspective, an IPA estimator is simply a sample path functional. Once its definition is fixed (using a particular representation) how the sample paths on which it acts are generated is immaterial. A similar point is made in [2] and [8], where it is shown that derivatives of real-valued random variables can often be generated using an expression derived in [8] - based on viewing the random variables as having been generated by the inversion method - even if the random variables are actually generated by a different method. In our birth and death example, the IPA algorithm obtained using the queueing network representation could be used on a simulation that follows the representation in [4]: given a sample path generated by any means we could apply our IPA estimator by viewing the path as though it had been generatedusing our representation. It is not hard to see from the algorithm in Section 3 that the sample derivatives calculated from a regenerative cycle depend only on the state sequence and the sojourn times. Provided these have the correct joint distribution, IPA estimates obtained from them are independent of how they were generated. t
Acknowledgments I am grateful to P~'~fessor Y.C. Ho, W.B. Gong, P. Vakili and Professor R. Suri for helpful discussions; and to an anonymous Associate Editor for a careful reading of the manuscript.
$. A note on implementation Thus far, we have emphasized that different representations lead to different IPA algorithms and different derivative estimators. As a final remark, we point out that although the choice of representation is important in the derivation of the algorithm, for the actual implementation it is usually the case that any simulation based on some valid representation will do - provided the IPA a l g o r i ~ is incorporated into the simulation as though the 'right' representation were being used, The crucial point is ~ a t in IPA no parameter change is ever really ~trodUced; hence, the actual sample path by sample path effect of 48
References
[1] X.R. Cao, "Convergenceof parametersensitivityestimates in a stochastic experiment", IEEE Transactions on Automatic Control AC-30,845-853(1985). [2] P.W. Glynn, "Process-differentiable representations of parametric families of random variables", Technical Report, Department of Industrial Engineering, Univet.~ityof Wisconsin-Madison,Wl, 1987. [3] W.B.Gongand Y.C. Ho, "Smoothedperturbation analysis of discrete-eventdynamicsystems", IEEE Transactions on Automatic Control AC-3Z856-868 (1987). [4] P. Heidelberger,X.R. Cao, M. Zazanis and R. Suri,"Convergence properties of infinitesimalperturbation analysis estimates", subn~tt=d for pubEca~n, 1987(revisedversion of IBM Research Report RCl1891 (53575)).
Volume 7, Number I
OPEI~ATIONS RESEARCH LETrERS
[5] Y.C. Ho and X.R. Cao, "Perturbation analysis and optimization of queueing networks", Journal of Optimization Theory and Applications 40 559-582 (1983). [6] Y.C. Ho and S. Li, "Extensions of perturbation analysis of discrete-event dynamic systems", IEEE Tra~sactions on Automatic Control, forthcoming May, 1988. [7] Y.C. Ho and P.Q. Yang, "Equivalent network, load-dependent servers and perturbation analysis - An experimental
February 1988
study", Proc. Int. Conf. Teletraffic Anal. and Comp. Perf. Ev., North-Holland, New York, 1986. [8] R. Suri, "Infinitesimal perturbation analysis for general discrete event systems", Journal of the Association for Computing Machinery 34 686-717 (1987). [9] M. Zazanis and R. Suri, "Estimating second derivatives of performance measures for G / G / 1 queues from a single sample path", submitted for publication, 1985.
49