Cyclic regenerative method of simulation

Volume I, Number 6

OPERATIONS RESEARCH LETTERS

December 1982

CYCLIC R E G E N E R A T I V E M E T H O D OF' S I M U L A T I O N

Robert G. SARGENT Deparm:eat of Industrial Enl~inecrinl~and Operations Research. Syracuse Unlucrsity. Syracuse. N Y 13210. U.S.A.

J.G. SHANTHIKUMAR gy.ffcnt~and IndustrialEnRmeerml¢ Umcer.tityo/Ari:ona. Tucson, A Z 85721, U.S.A. Received June 1982 Revised November 1982

In this paper we develop a cyclic regenerative method of simulation. As an example of its application, we consider a single server queuing system with N-control policy of Heyman. Numerical results are presented to illustrate the ability of this method to reduce the variance of the estimate over that of the traditional regenerative method of simulation. ReMnerative precesses, regenerative method of simulation, variance reduction techniques, controlled queue

I, Introduction The regenerative method of simulation is one of the major contributions made during the 1970's towards the theory of simulation output analysis (see, for example, Crane and lglehart [2], Fishman [4], 181ehart [7], Lavenber8 and Slutz [8], and others). This method has its own advantages and disadvantages (see Schruben [! I]) and several studies have been conducted to improve this hpproach, in general, attempts have been made to reduce the bias in estimates (see Meketon and Heideiberger [10], for example) and to reduce the variance of the estimates (see Lavenberg et al. [9], for example). Attempts to generalize the method of regenerative simulation are, however, limited. Efforts have been made to generalize this method to non.regenerative processes, thus introducing the almost-regenerative method of simulation (see [3l and [5]). The direct opposite to this case is the situation where the process under consideration has more than one regeneration point, it is our purpose in this paper to develop a method to utilize such a special property. in this paper we consider the cyclic regenerative process and develop a method of cyclic regenerative simulation. The cyclic regenerative process is described in Section 2, and the cyclic regenerative method of simulation is developed in Section 3. A numerical example illustrating the ability of this method to reduce the variance of the estimate over that of the traditional regenerative method of simulation is given in Section 4. 2, Cyclic repnerative process The cyclic regenerative process is a special case of regenerative processes. Some useful properties of cyclic regenerative processes are studied in Shanthikumar [12]. We first describe this process and then present some useful results that will be used later to develop the method of cyclic regenerative simulation. 2.1. Continuous time cyclic regenerative process

Let (V(t), t >~O) be a continuous time regenerative stochastic process with regeneration epochs ~'o-0 < % < .-. < % < . . . . For the cyclic regenerative process, we assume that there are r - ! (r > I) 222

0167-6377/82/0000-0000/$02.75 © 1982 North-Holland

Volume !, Number 6

OPERATIONS RESEARCH LETTERS

December 1982

epochs a~j, i = !, 2 , . . . , r - i , during the j t h regeneration cycle (,~_=, ~], V j = 1, 2 ..... such that ~_ ~ < a~j < . . . < a,_ ~,j < ~) and V(a~j),j - 1, 2,..., have identical probability distribution with distribution function (d.f.) A~(. ), i = !, 2 , . . . , r - I. It is also assumed that V ( t ) , t > a~j given V(a~j) is independent of V(t), t < a~j for all i - !, 2 , . . . , r - 1, a n d j = !, 2, .... Let V(t[z) be the conditional value of V(t + ao) given that 1/'(a~j) - z. Then V(ilz) a__ V(t + a~j)lV(a~j) = z, has identical probabilistic properties for all j = I, 2,..., for a fixed z, t > 0, and i. That is, these intermediate epochs a o serve as weak regenerative epochs and we will call (V(t), t >i 0) satisfying these properties as the cyclic regenerative process. The time intervals ( ~ - i , alj] will be called phase 1, (a~_l. j, a~j] will be called phase i (for i -- 2 , . . . , r - !) and the time interval (a,_ =,j~] will be called phase r, for all j = 1, 2, .... Let lim~_.=V(t)=* ° V and we wish to obtain E [ f ( V ) ] for some real function f ( . ) . Let Ya be the total of f ( V ( t ) ) during phase i (i = I, 2,..., r) of the first regeneration cycle. That is,

Ye-fa~".

i-1,2,.... r,

with a01 = 0 and at! = I"I. Also let X t = all -- ~), X~ = a~t - ai_ I,t, i - 2, .... r - I, and X, = ~t - a,_ I.l" Then, from the theory of regenerative processes (refer to Cohen [I ]), we have

(1) 2.2. Discrete time cyclic regenerative process A discrete time regenerative process (V,. n = I. 2. 3 .... ) with respect to the integer valued renewal process (N,}~ is a cyclic regenerative process if (N,)~ has an imbedded integer valued cyclic renewal process (N,~, i = !, 2,..., r)~ for some r > !, and the occurrence of a type i renewal process can be regarded as a weak regeneration point for (V~, n = !, 2 , . . . ) similar to the continuous time alternating regenerative process. Then, as in the previous case,

where E[ X~] = E[ N,,], i = I, 2,..., r.

3. Cyclic regenerative simulation In this section we develop the method of cyclic regenerative simulation for processes satisfying the above conditions, and for which the distribution function A~(.), for some values (say k, k >I I) of i (i = I, 2,..., r - !) are known. Then we will look at this process as a cyclic regenerative process with only k + I phases. Therefore, without loss of generality, we will assume that the distribution functions A~(.) for all i = 1, 2,..., r - I, are known. This enables us to simulate the phases independently from one another. That is, if we wish to simulate phase i + I independently of others, we initiate this phase by sampling a value from A~(.) and continuing the simulation until this phase is over. We will consider two cases. The first case is where we have analytical results for some phases of the alternating regenerative process. In the second case, we assume that we do not have analytical results for any phase.

Case 1. Hybrid simulation~analytic method Let D be the set of phases for which we have analytical results. Let F = (I, 2 , . . . , r ) - D be the set of phases for which we do not have analytical results. Let ((Yo, x , j ) , j = I, 2,...,n~; i ~ F ) be n~ samples of ((Y,, X~)}, i ~ F, obtained by independently simulating the ph~,~e~ for which we do not have analytical 223

OPERATIONSR E S E A R C H LETTERS

Volume I,Number 6

December 1982

results. Then, using the standard methods, it can be easily established from (!) or (2) that

Z E[Y,] + Z g(n,) i~ i =

lED

isF

F.. ~[x,l + I". ~ , ( . , ) iGD

•

'

(~)

iGF

where l

n,

yi( ni ) ~ii j~dlyi/ and ~i( n~) =

=

x O, l nj~ , I

--

i ~ F

(4)

n~

is a strongly consistent estimate of E[f(V)]. Now, noting that I : E[ Y,] + ~ E[~(.,)]

e[.f(v)]-"° I2

"~

(5)

~[x,l+ Y'. ~[ g(.,)]

itD

iGF

a,~d considering the expectation and variance of

y.. (e[r,l- e[ /(v)] ~[ x,]) + Z (~(.,)- e[ /(v)] ~,(.,)), i¢=D

iGF

we can show that an approximate (I - B) 100~ confidence interval for E[ f( V)] is

(~,- !,(., 8),,~,

+

!,(,,. 8)).

(6)

where

(7) iGD

~,.

iGF

"

~ { (~,,, + (~, .~ ,. 2~,4,)/.,}. , ~,.,-

(8)

~EF n,

s_.~ - n~ -'

! j~. ,

(y.-~,(.,))~,

(9)

n,

I: (~,j- ~,(.,))~, I

.,

(no)

n,

~,,, ".,=7 /-I I: (y,,-~,(.,))(~,j- ~,(.,)), and if Z is a unit normal random variable, P[-Zs/2

(nn)

< Z < Zs/2] = ! - B .

Case 2. Independentsimulation of phases in Case 2, we assume that we do not have analytical results for any phases. For this case, let ((y,j, x,j), j - I, 2,..., n~; i - I, 2,..., r) be n~ samples of flY, X~)), i = i, 2,...,r, obtained by simulating all r phases

independently. Also, let ((.~, x~), j - - I, 2,..., n ) , b e , s~mples of the total, (( Y = ]~.. t Y~), (X = ~ . nX,)) during a complete regeneration cycle of the r phases, obtained simulating the r phases together in one regeneration cycle for n times. Then, similar to the previous case, we have P

~ g(n,) +~(,,) /~./2 ~

i--I

•

i--I

224

,

(n2)

Volume !, Number 6

OPERATIONS RESEARCH LEITERS

December 1982

w h e r e .~(n,) and ~,(n+) are as defined in (4), a n d

.~(n) = n

~.,

~(n)-- n

j-t

(13)

xj. j=t

N o w noting that r

~:[/(v)]

E E[~,(.,)] +E:[ P(.)]

=

'~'

(14)

E E[ ~',(.,)] + ~:[ ~'(.)]

i-,,I

a n d considering the m e a n and variance of r

Z (~,(.,)-~.[/(v)] Y~,(.,))+ ~(.)-~.[f(v)] y~(.).

#,=1

we can show that an approximate (! - ~ ) ! 0 0 ~

confidence interval for £ [ f ( V ) ]

is

(,~,- ~(., ~), ,~, + t=(., ~)),

(t5)

where

(16) r

~-

E { ( s p , , + ( + ,) ~s==, ~ - 2,~,sp~,)/,,,)) + (+;, +

(~)'s~-

(17)

2m=sp+)/,,

i-I

Table I Variance reduction dt for Case I (hybrid simulation/analytic method) with n 2 = n and (variability in d I ) over 25 replications N

Inter+at+rival/Service time distribution

n

10

20

Ea/E=

M/M

H2/H2

P

P

0

0.5

0.7

0.5

0.7

0.5

0.7

$0

$7.9~ (10.6)

29.2~ (10.2)

63.4C$ (8.7)

19.61£ (9.2)

50.77 (16.4)

I I. I (i I.i)

100

57.3

26.0~

59.6~

19.6~

(8.8)

(9.6)

(8.3)

(9.7)

49.4~ ( ! 4.2)

(5.3)

200

$7.6~ (5.6)

26.7~ (7.4)

59.9¢£ (5.0)

17. i ~ (6.9)

49.0¢£ (9.9)

7.9~ (4.7)

50

$2.3~ (12.6)

30. I =£ (7.7)

$9.27 (8.5)

27.2~ ~9.2)

60.67 (9.0)

23.7=£ (I 1.7)

100

53.9 (7.3)

26.7=$ (6. I)

57.6~ (5.7)

26.4~ +'7.3)

58.0~ (7.5)

i 8.5~, (9.0)

200

$2.9~;

28.$~

57. i ~

25.2~

58.7~,

16.4~;

(5.7)

(4.1)

(5.3)

(5.2)

(5.5)

(5.8)

9.7~

E2 - Erlang-2 with squared coefficient of variation !/2. M - Exponential. H + - Hyper-exponential of two phases with squared coefficient of variation 2. 225

Volume I, Number 6

OPERATIONS RESEARCH LETTERS /./

tl

]

December 1982

E

S~l=n - I .

./--I 1

fl

and s21,, s=2~ , 2 and s22~ are as defined by equations (9), (10), and (I 1), respectively. Next we will consider three special cases of Case 2.

Case 2.1. No independent simulation of phases Suppose n t - n = - ... - n , - 0. Then this case will reduce to the traditional method of regenerative simulation. For this case, a point estimate rh~ and an approximate (I - 8 ) i 0 0 ~ confidence interval of (n'l~- I~(n, ~), fli.~+ i~(n, 8)) follows directly from (12) and (15), respectively. Table 2 Variance reduction d4 for Case 2.2 (only independent simulation of phases with unequal number of repetitions n l and n2) with n - min(n,, ha) and (variability in d4) over 25 replications

N

!0

nI

n:

Inter-arrival/Service time distribution

E2/E2

M/M

H=/H:

P

p

p

0.5

0.7

0.$

0.7

0.S

0.7

$0

100

i 1,3 (2,5,6)

34.5 (32.6)

23..5 (22.2)

13.4 (140.0)

2.5.2 (27.71

46.3 ( 19.1)

.SO

200

29,5 (30,3)

64.6 (12.1)

24.9 (24.1)

.58.9 (25.8)

34.2 (23.4)

39. I (68.1)

|00

200

15, ! (14,3)

47.2 (I !.3)

30.6 (14.6)

4.5. I (19.9)

29.9 (21.8)

38. I (27.0)

100

50

24,5 (15.5)

18..5 (8.10)

31.4 (16.3)

! 3. I (10.1)

2.5.4 (16.$)

200

SO

42,2

28.2

(10.8)

(8.8)

43.0 (12.,5)

16.0 (19.3)

43.4 (12.9)

13.6 (I 1.4)

14.5 (27.8) 31,86 (23,1)

18.2 (7.,5) 40.7 (I 3.4) 65.4 (11.7)

33.7 (10.8) 19.4 (24.8) 36.6 (21.2)

14.9 (7.8) 41.9 (I 3,8) 64.8 (13.1)

27.8 (10.7) 15.7 (24.7) 33. I (24.7)

7.$6 (4.9) 30.2 (50.8) $4.0 (32.7)

200

!00

25,9

(10.2) 20

SO

I00

50

200

i00

200

20.9 ( 14. I )

43.5 (12.0)

2S.7 (15.4)

38.6 (25.7)

24.7 ( i 5.7)

41.4 (36.6)

100

50

23.2 (20.7)

17.4 (11.6)

29.3 (IS.I)

16.6 (S.S)

24.6 (14.4)

16.2 (9.2)

200

SO

36,8 (13.7)

25.9 (9.2)

40.2 (13.7)

25.9 (10.2)

44..5 (12.2)

18.01 (10.2)

I O0

26.9 (I !, ! )

17.7 (6.8)

18. i (5.4)

30.6 ( I0. I )

14.9 (?. I )

200 :

31.9 (7.80)

Ez - EdanS-2 with s q ~ coefficient of variation I/2. M - Exponential. H2 - Hyper-exponential of two phases with squared coefficient of v~riation 2.

226

6.32 (6.63)

Volume i, Number 6

O P E R A T I O N S RESEARCH LETTERS

December 1982

Case 2.2. Only independent simulation of phases In this case, we set n - 0. That is, we have a pure cyclic regenerative simulation, where all phases are simulated independently. Then, a point estimate rh4 and an approximate (! - 8 ) 1 0 0 ~ confidence interval of ( m 4 - 14(n, 8), rh4 + 14(n, 8)) follows directly from (12) and (15), respectively.

Case 2.3. Only independent simulation of phases with equal numbers of repetitions In order to make a meaningful comparison of the above case to the traditional method of regenerative simulation (Case 2.1), we should let n I - n 2 - . . . . n, - n in Case 2.2. Thus, this case is the same as Case 2.2, except we set n I - n 2 = . . . = n , - n. Then a point estimate ms and an approximate (I - 8 ) 1 0 0 ~ confidence interval of ( ~ s - Is(n, 8), ms + Is(n, 8)) follows directly from (12) and (15), respectively. For the purp(,se of comparison of these different cases, an example of a single server queuing system with N-control policy of Heyman [6] is considered.

4. Example The example considered here is a single server queue with N-control policy of Heyman [6]. Under an N.control policy, the server is deactivated when there is no customer in the system and activated when there are N customers in the system. The arrival of the first N customers after server deactivation forms the control phase (call it phase 1) and the Nth arrival will activate the server. The server is next deactivated as soon as the system becomes empty. The time between the server activation and the next deactivation is the busy phase (call it phase 2), during which the server is active. Thus, the system alternates between control and busy phases.

Table 3 Variance reduction ds for Case 2.3 (only independent simulation of phases with equal number of repetitions n I = n 2 -. n) and (variability in d s) over 25 replications N

10

20

Inter.arrival/Service time distribution

n

E,/E,

M/M

H2/H2

P

P

P 0.5

0.7

8.5c~ (8.6)

7.2~o (22.9)

3.79/, (I 1.0)

10.6~ (14.9)

10.4~ (6.8)

6.0~

5.4~

(1.5.4)

(4.5)

12.1~o (5.6)

9.19/, (9.6)

9.1~, (4.6)

8.8~ (9.0)

4.5~ (3.7)

I !.5% (9.2)

4.9~ (20.8)

I 1.9~ (10.9)

6.0~ (! 9.0)

12.59/, (9.8)

3.9~

I 1.0fo

0.5

0.7

0.5

50

- 7.4~ (25.7)

i0.4~ (10.9)

10.2~ (20.7)

I00

- 5.8~

I !. I ~

(15.3)

(6.7)

200

- 3.5~ (14.5)

50

- !.6~ (27.4)

100

- 1.3fo

200

9.3c£

0.7

3. i ~

9.9~

(15.2)

(8.9)

(15.3)

(5.4)

(15.3)

(5.9)

- 1.0~ (12.4)

12.4~ (3.5)

5. i ~ (10.7)

I 1.5~ (3.4)

5.5~ (12.0)

9.5~ (3.9)

E 2 - Erlang-2 with squared coefficient of variation I / 2 . M - Exponential. H 2 - Hyper-exponential of two phases with squared coefficient of variation 2.

227

Volume !, Number 6

OPERATIONS RESEARCH LETTERS

December 1982

Assume that the system is empty at time zero and let IV, be the time spent in the system by the n th customer. Then (see Shanthikumar and Balci [13]), N

W,,=

E jmn÷|

n

Aj+

EBb, n-l,2,...,N-l.j--|

N

and

W.- )".Bj, j----|

where Bj is the service time of the nth customer and Aj is the interarrival time between the ( j - l)th and j t h customer, it is easily seen that the time spent by the first and Nth customer to enter the system after every server.deactivation has the same probability distribution as that of Wt and W,, respectively. Hence it is easily verified that (W,, n ,, !, 2,... ) is a discrete time cyclic regenerative process with two phases. The first phase is formed by the first N customers to arrive after server deactivation and the second phase is formed by the customers who arrive during the busy phase. This system is modelled and simulated using the different cases of the cyclic regenerative method of simulation to estimate the mean time spent in the system. Define the percentage variance reduction dr as ( I - d~/d~)× 100% over the traditional method of regenerative simulation. The average and standard deviation of the variance reduction dr, obtained over 25 replications, for different cases of interarrival/service time distributions and server utilization p and control value N are given in Tables ! to 3, The standard deviation of dr is given within ( ) .

S. Conclusion

The numerical results presented in Tables I to 3 clearly show that the cyclic regenerative method of simulation performs better than the tradi',ional method of regenerative simulation, in almost all cases considered. The hybrid simulation/analytic method is clearly a very good approach (see Table 1) and this approach has been discussed in Shanthikumar and Balci [I 3] and Shanthikumar and Sargent [14]. The other two cases, Case 2.2 and Case 2,3, have their potential (see Tables 2 and 3), and promise to be a fruitful area for future reseat'oh. Particularly, the application of other variance reduction techniques to this cyclic regenerative method of simulation should be explored, it should also be noted that these ideas can be extended to regenerative processes with underlying Markov-renewal processes as opposed to the cyclic renewal process considered in this paper.

Aeknowledlpnent Appreciation is given to Mr. R.V. Arun for his computatioaal assistance in the numerical example.

References [I] J,W. Cohen, On ileaenemtieeProcessesin Queuing Tkeory. Lecture Notes in Economics and Mathematics 121. Springer, Berlin (1976), [2] M,A, Crime and D,L, 18lehart, "Simulating stable stochastic systems !11: Regenerative processes and discrete event simulation", Opeeutions Rer, a3, 34-4,5 (1975).

[3] M,A, Crane and D.L, ll0ehurt, "Simulating stable stochastic systems, IV: Approximate techniques", Management $ci. 21, 1215--12,'M(1975), |4l G.S. Fishn~__~n, Principlesof Discrete Et~ent Simulation, Wiley,New York (1978). [5] F.L Guntherand R.W.Wolff,"The almostregenerativemethod for ~tochasticsystemsimulations", Operations Res. 28, 375-386 (1980). [6] D.P. Heymun, "OW~mal operating policies for M / G / I queueins systems", OperationsRes. 16, 362-382 (1968). [7J D.L I81ehart, "The regenerative method for simulation analysis", in: K.M. Chandy and R.T. Yen, Eds., Current Trends in Programming Metkodolo~w. Vol. I!!: Software Mode#in& Prem~ce-Hall, Enalewood Cliffs, NJ (1978) 52-71. [8] S.S. Lavenburg and D.R. Slutz, "Introduction to'regenerative simulation", IBM J. Res. Development 19, 458-46'2 (1975). 228

Volume !, Number 6

OPERATIONS RESEARCH LETTERS

December 1982

[9] S.S. Lavenburg, T.L. Moeller and C.H. Sauer, "Concomitant control variables applied to the regenerative simulation of queueing systems", Operations Res. 27, 134-160 (1979). [10] M.S. Meketon r,nd P. Heidelberger, "A renewal theoretic approach to bias reduction in regenerative ~mulation", Management $ci. 28, 173-181 (1982). [I I] L.W. Schruben, "A coverage function for interval estimation of ~imulation response", Management $ci. 26, 18-27 (1980). [121 J.G. Shanthikumar, "Some properties of the alternating regenerative processes", Working Paper O79-018, Department of Industrial Engineering, University of Toronto, Toronto (! 979). [13] J.G. Shanthikumar and O. Baici, "A variance reduction technique for the simulation of a single server G I / G / ! controlled queue", Proc. Summer Computer Simulation Conference (1980) 105- 110. [14] 3.(3. Shanthikumar and R.G. Sargent, "A unifying view of hybrid sJ.~.ulation/analytic models and modelling", Working Paper #82-003, Department of Industrial Engineering and Operations Research, Syracuse University, Syracuse (1982).

229