Efficient implementations of the randomization method with control of the relative error
Abstract
Randomization is a well-known numerical method for the transient analysis of continuous-time Markov chains. The main advantages of the method are numerical stability, well-controlled computation error and ability to specify the computation error in advance. Typical implementations of the method control the truncation error in absolute value, which is not completely satisfactory in some cases. Based on a theoretical result regarding the dependence on the parameter of the Poisson distribution of the relative error introduced when a weighted sum of Poisson probabilities is truncated by the right, in this paper we develop efficient and numerically stable implementations of the randomization method for the computation of two measures on rewarded continuous-time Markovchains with control of the relative error. The numerical stability of those implementations is analyzed using a small example. We also discuss the computational efficiency of the implementations with respect to simpler alternatives.
Full text
Computers & Operations Research 32 (2005) 1089–1114 www.elsevier.com/locate/dsw E$cient implementations of the randomization method with control of the relative error V( )ctor Su˜n(e∗, Juan A. Carrasco Departament d’Enginyeria Electr onica, Universitat Polit ecnica de Catalunya, Diagonal 647, plta. 9, Barcelona 08028, Spain Abstract Randomization is a well-known numerical method for the transient analysis of continuous-time Markov chains. The main advantages of the method are numerical stability, well-controlled computation error and ability to specify the computation error in advance. Typical implementations of the method control the truncation error in absolute value, which is not completely satisfactory in some cases. Based on a theoretical result regarding the dependence on the parameter of the Poisson distribution of the relative error introduced when a weighted sum of Poisson probabilities is truncated by the right, in this paper we develop e$cient and numerically stable implementations of the randomization method for the computation of two measures on rewarded continuous-time Markov chains with control of the relative error. The numerical stability of those implementations is analyzed using a small example. We also discuss the computational e$ciency of the implementations with respect to simpler alternatives. ?2003 Elsevier Ltd. All rights reserved. Keywords: Rewarded continuous-time Markov chains; Transient analysis; Randomization; Relative error 1. Introduction Rewarded continuous-time Markov chains (CTMCs) are a powerful modeling formalism. A rewarded CTMC is a CTMC X={X(t); t¿0}with a reward structure imposed over it. The reward structure may include reward rates associated with states and impulse rewards associated with transitions. In this paper, we will consider rewarded CTMCs with a reward structure including only reward rates associated with states. Let be the state space of Xand let ri,i∈be the reward rate associated with state i. We will assume ri¿0, i∈. The quantity rihas the meaning of “rate at which reward is earned while the CTMC is in state i.” The behavior with time of the random ∗Corresponding author. E-mail addresses: [email protected] (V. Su˜n(e), [email protected] (J.A. Carrasco). 0305-0548/$ - see front matter ?2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2003.09.014
1090 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 variable rX(t)can be quantiFed using several measures. In this paper, we will consider two such measures: the “expected transient reward rate,” ETRR(t)=E[rX(t)], and the “expected averaged reward rate,” EARR(t)=E[t 0rX()d=t]. Those measures with particular instances of the CTMC X and the reward rate structure ri,i∈have important applications, particularly in the dependability and performability analysis of fault-tolerant computer systems. Consider, for instance, a fault-tolerant computer system which can be “up” or “down” whose evolution is modeled by a CTMC Xwith state space =U∪D, where Uincludes the states in which the system is up and Dincludes the states in which the system is down. Then, with a reward rate structure ri=0, i∈Uand ri=1, i∈D, the ETRR(t) measure would be the unavailability of the fault-tolerant computer system at time t, i.e. the probability that the system is down at time t. As another example, consider a degradable fault-tolerant multiprocessor system whose behavior is modeled by a CTMC X. Then, with a reward rate structure assigning to each state of Xthe speedup of the multiprocessor in that particular state, the ETRR(t) measure would be the expected speedup of the multiprocessor at time tand the EARR(t) measure would be the expected averaged speedup of the multiprocessor during the time interval [0;t]. Computation of the measures ETRR(t) and EARR(t) requires the transient analysis of the CTMC X. That transient analysis can be performed using ordinary diIerential equation solvers and randomization [1,2] (also called uniformization). Good properties of the randomization method are excellent numerical stability, well-controlled computation error, and ability to specify the computation error in advance. The randomization method is based on the randomization result. To review that result, let i;j,i; j ∈,i=j,i=j∈−{i}i;j,i∈, and A=(ai;j)i;j∈,ai;j =i;j,i=j,ai;i =−ibe, respectively, the transition rates, output rates and transition rate matrix of X. Let ¿maxi∈iand consider the discrete-time Markov chain (DTMC) ˆ X={ˆ Xk;k=0;1;2;:::}with same state space and initial probability distribution as Xand transition probability matrix P=I+(1=)A, where I denotes an identity matrix of appropriate dimension. The DTMC ˆ Xis said to be the randomized DTMC of Xwith randomization rate . Let Q={Q(t); t¿0}be a Poisson process with arrival rate independent of ˆ X. We have P[Q(t)=k]=Pk(t), with Pk()=e −k=k!, where Pk()isthe probability mass function of a Poisson random variable with parameter . Then, the randomization result states [3, Theorem 4.19] that Xis probabilistically identical to {ˆ XQ(t);t¿0}, which means that the probability of any event deFned over the values of the random variables X(t), t¿0isthe same as the probability of the corresponding event over the random variables ˆ XQ(t),t¿0. We will review next typical implementations of the randomization method for the computation of the ETRR(t) and EARR(t) measures with control of the absolute truncation error. Since the computational cost of the method increases with the randomization rate ,is taken equal to maxi∈i. We will start by the ETRR(t) measure. Using the randomization result, the transient regime of Xcan be expressed in terms of the transient regime of ˆ Xas P[X(t)=i]= ∞ k=0 P[ˆ Xk=i|Q(t)=k]P[Q(t)=k]= ∞ k=0 P[ˆ Xk=i]Pk(t); and using ETRR(t)=i∈riP[X(t)=i]: ETRR(t)= i∈ ri ∞ k=0 P[ˆ Xk=i]Pk(t)= ∞ k=0 d(k)Pk(t);(1)
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1091 with d(k)= i∈ riP[ˆ Xk=i]: In the implementation, the inFnite series (1) is truncated by the right to obtain an approximate value for ETRR(t): ETRRa R(t)= R k=0 d(k)Pk(t):(2) Letting rmax = maxi∈ri, and using d(k)6rmax, it follows that rmax ∞ k=R+1 Pk(t) upper bounds the truncation error and, then, Ris chosen using: R= min m¿0: rmax ∞ k=m+1 Pk(t)6; where ¿0 is the requested absolute error. The required probability row vectors of ˆ X,q(k)= (P[ˆ Xk=i])i∈,k=0;1;:::;R are obtained using q(0) = , where =(i)i∈,i=P[X(0) = i] and q(k+1)=q(k)P:(3) The truncation error upper bound increases with t, because it is rmax times the probability that the number of arrivals in the time interval [0;t] in a Poisson process with arrival rate is ¿R+1. Then, if ETRR(t) has to be computed for several values of twith absolute error 6, it is enough to control the truncation error for the largest of them. We will review next a typical implementation of the randomization method for the EARR(t) measure with control of the absolute truncation error. Using EARR(t)=(1=t)t 0E[rX()]d=(1=t) t 0ETRR()d,(1) and t 0Pk()d=(1=)∞ l=k+1 Pl(t) (see, for instance, [4, Formula 14.512]), we can obtain EARR(t)=1 tt 0 ∞ k=0 d(k)Pk()d=1 t ∞ k=0 d(k)t 0 Pk()d=1 t ∞ k=0 d(k)∞ l=k+1 Pl(t) =1 t ∞ k=1 k−1 l=0 d(l)Pk(t)= 1 t ∞ k=1 k−1 l=0 d(l)t kPk−1(t) =∞ k=0 1 k+1 k l=0 d(l)Pk(t):(4) In the implementation, the summatory is truncated by the right to obtain an approximate value for EARR(t): EARRa R(t)= R k=0 1 k+1 k l=0 d(l)Pk(t):
1092 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 Since (4) is formally identical to (1) with d(k) replaced by d(k)=k l=0 d(l)=(k+ 1) and, therefore, d(k)6rmax, the truncation point Ris chosen as in the method for the ETRR(t) measure. The required probability row vectors q(k) are obtained also as for the ETRR(t) measure. Similarly, if the measure has to be computed at several time points, it is enough to control the truncation error for the largest of them. Alternative, similar implementations of the randomization method for the computation of EARR(t) can be found in [5]. Computation of the Poisson probabilities Pk(t)=e −t(t)k=k! is a delicate issue due to the possibility of intermediate overQows and underQows. An approach to compute “exact” Poisson probabilities is to use the method proposed in [6, pp. 1028–1029] (see also [7]), which avoids intermediate overQows and underQows and has good numerical stability. The approach is, however, relatively expensive. The computational cost of the (standard) randomization method increases with R. Since (see, for instance, [8, Theorem 3.3.5]) the random variable Q(t) has for t →∞ an asymptotic normal distribution with mean and variance t, for large t and 1, the truncation point Rwill be approximately equal to t, implying that, for large t and 1, the method will be expensive. The high computational cost of the method for large t and 1 has motivated the development in the last years of several variants which can outperform the standard randomization method: selective randomization [9,10], multistepping [11, Section 3.1.2], adaptive uniformization [12], adaptive/standard uniformization [13], uniformization with steady-state detection [14,15], and regenerative randomization [16,17]. Randomization-based methods computing bounds which can be much more e$cient than “exact” methods have also been proposed [18]. Despite all these variants, the standard randomization method is still competitive for many rewarded CTMC models. The implementations of the standard randomization just reviewed are, in some cases, not completely satisfactory. The reason is that, in some cases, the user is interested in computing the measure with some reasonably small relative error, while the implementations control the absolute truncation error. In those cases, unless the user knows in advance a good estimate for the order of magnitude of the measure, he/she will be confronted with the dilemma between choosing a very small error control parameter to ensure that the desired accuracy will be satisFed, with the consequent increase in the computational cost (the truncation point Rincreases for decreasing ), or choosing a greater value for and running the risk of obtaining not enough accurate computations and having to run the method again. For very large t, the price paid by using an smaller than strictly required is relatively small. However, for not too large t, the price can be relatively important. This is illustrated in Table 1which gives the required Ras a function of and t, assuming rmax =1. In this paper, we will develop e$cient and numerically stable implementations of the standard randomization method for the ETRR(t) and EARR(t) measures with control of the relative error. The rest of the paper is organized as follows. Section 2will obtain a theoretical result concerning the dependence on the parameter of the Poisson distribution of the relative error introduced when a weighted sum of Poisson probabilities is truncated by the right. Based on that result, Section 3will derive the implementations. Section 4will analyze the numerical stability of the implementations using a small rewarded CTMC model of a repairable fault-tolerant computer system and will discuss the computational e$ciency of the implementations with respect to simpler alternatives. Finally, Section 5will present some conclusions. The Appendix includes some proofs.
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1093 Table 1 Truncation point Ras a function of and t for rmax = 1 in the reviewed implementations of the standard randomization method with control of the absolute truncation error t =10 −4=10 −6=10 −8=10 −10 =10 −12 0.01 1234 4 0.05 2345 6 0.13456 7 0.5 5 7 8 10 11 1 6 91112 14 515192225 27 10 24 28 32 36 39 50 78 87 94 101 107 100 139 151 161 170 178 500 585 610 630 649 665 1000 1120 1154 1182 1208 1230 2. A theoretical result In this section we will develop a theoretical result concerning the dependence of the relative error with respect to the parameter introduced when a weighted sum of Poisson probabilities ∞ k=0 !(k)Pk() is truncated by the right, i.e. letting Wa R()= R k=0 !(k)Pk() (5) and We R()= ∞ k=R+1 !(k)Pk() (6) we will show that, for R¿0 and assuming !(k)¿0 and !(k)¿0 for some k,06k6R, eR()=We R()=W a R() (7) increases with ,¿0. We will start by the following lemma. Lemma 1. Let 1¿2¿0, x=1=2,j¿i¿0, and !(k)¿0, k¿0. Then 1 xje1−1 x1 j k=i !(k)Pk(1)6 j k=i !(k)Pk(2)61 xie1−1 x1 j k=i !(k)Pk(1): Moreover,if !(k)is uniformly upper bounded,the right inequality also holds for j=∞.
1094 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 Proof. If !(k)=0 for k,i6k6j, the result is trivial. Therefore, let us assume !(k)¿0 for some k,i6k6j. Substituting 1=x2into j k=i!(k)Pk(1) yields j k=i !(k)Pk(1)= j k=i !(k)e−1k 1 k!= j k=i !(k)xke−x2k 2 k! =e −(x−1)2 j k=i xk!(k)e−2k 2 k!=e −(x−1)2 j k=i xk!(k)Pk(2): Then, since x¿1 and !(k)¿0, xie−(x−1)2 j k=i !(k)Pk(2)6 j k=i !(k)Pk(1)6xje−(x−1)2 j k=i !(k)Pk(2): Inverting the members of the previous inequalities (this can be done because !(k)¿0 for some k,i6k6jimplies j k=i!(k)Pk(1)¿0 and j k=i!(k)Pk(2)¿0 and, therefore, that all three members are ¿0): e(x−1)2 xjj k=i!(k)Pk(2)61 j k=i!(k)Pk(1)6e(x−1)2 xij k=i!(k)Pk(2): Multiplying by j k=i!(k)Pk(2)j k=i!(k)Pk(1)¿0, 1 xje(x−1)2 j k=i !(k)Pk(1)6 j k=i !(k)Pk(2)61 xie(x−1)2 j k=i !(k)Pk(1); and, noting that (x−1)2=(1−1=x)1, 1 xje1−1 x1 j k=i !(k)Pk(1)6 j k=i !(k)Pk(2)61 xie1−1 x1 j k=i !(k)Pk(1): Finally, if !(k) is uniformly upper bounded, from j k=iPk(), j→∞,¿0 being upper bounded and !(k)¿0, it is easy to prove that j k=i!(k)Pk(), j→∞,¿0 is upper bounded and, being increasing, convergent, implying that both j k=i!(k)Pk(2), j→∞ and (1=xi)e(1−1=x)1j k=i!(k) Pk(1), j→∞will converge and that the previous right inequality will also hold for j=∞. The following theorem asserts the result which is the purpose of this section. Theorem 1. Let ¿¿0, R¿0and !(k)¿0, k¿0. Assume !(k)¿0for some k,06k6R and assume that !(k)is uniformly upper bounded. Then,eR()¿eR().
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1095 Proof. Let x== ¿1. Using (5) and the left inequality of Lemma 1with 1=,2=,i=0, and j=R, Wa R()= R k=0 !(k)Pk()¿1 xRe1−1 x R k=0 !(k)Pk()= 1 xRe1−1 x Wa R():(8) Similarly, using (6) and the right inequality of Lemma 1with 1=,2=,i=R+ 1, and j=∞, We R()= ∞ k=R+1 !(k)Pk()61 xR+1 e1−1 x∞ k=R+1 !(k)Pk()= 1 xR+1 e1−1 x We R():(9) Finally, using (7), (8), (9), x¿1, noting that !(k)¿0 for some k,06k6Rand ; ¿0, implies Wa R()¿0 and Wa R()¿0: eR()=We R() Wa R()6 1 xR+1 e1−1 x We R() 1 xRe1−1 x Wa R() =1 x We R() Wa R()6We R() Wa R()=eR(): 3. The implementations In this section, based on Theorem 1, we will derive implementations of the standard randomization method for the computation of the measures ETRR(t) and EARR(t) at several time points t1;t 2;:::;t n with control of the relative error. We will make the following assumptions: A1) ri¿0 for some state ireachable from some state jwith non-null initial probability. A2) ti¿0, i=1;2;:::;n. Assumption A1 implies d(k)¿0 for some k¿0. Note that if assumption A1 is not satisFed, d(k)=0 for all kand, then, ETRR(t)=EARR(t)=0 for all t. Further, for t=0, we have ETRR(t)=EARR(t)= d(0). Therefore, when any of the assumptions A1 or A2 is not satisFed, computation of the ETRR(t) and EARR(t) measures is straightforward. We will only discuss in detail the implementation for the measure ETRR(t). We will start by, using Theorem 1, developing in Section 3.1 an e$cient method for obtaining a global truncation point by the right Rguaranteeing a relative truncation error 61at every time point ti, where 1¿0isan error control parameter. Next, in Section 3.2, we will discuss the introduction of a truncation point by the left Liand a truncation point by the right R∗ i, beyond the truncation point by the right R, for each particular time point tiwith the purpose of reducing the number of Poisson probabilities which have to be computed. Those truncations will introduce a further relative truncation error 62, where 2¿0 is another error control parameter. In Section 3.3, we will show how weights proportional to the Poisson probabilities can be used to determine the truncation points Liand R∗ iand how, at the price of introducing an additional relative error, normalized weights can be used instead of “exact” Poisson probabilities to perform the actual computation of the truncated summatories in a
1096 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 numerically stable way. Also in that section, we will show how the error control parameters 1and 2can be adjusted so as to guarantee an absolute relative error 6for every time point ti, where ¿0 is the requested relative error. Finally, in Section 3.4, we will give algorithmic descriptions of the implementations for both the ETRR(t) and the EARR(t) measures. 3.1. Determination of a global truncation point by the right R According to the review of the implementation of the standard randomization method for the ETRR(t) measure performed in Section 1, a truncation point Rby the right guaranteeing a relative truncation error, ∞ k=R+1 d(k)Pk(t)=ETRRa R(t), 61at a given time point t¿0 can be chosen using R= min m¿0: ETRRa m(t)¿0∧rmax ∞ k=m+1 Pk(t)ETRRa m(t)61:(10) That a Fnite Rcan be determined in this way follows from: 1) because d(k)¿0 for some kand t¿0, ETRRa m(t)is(2)¿0 for some m,2)ETRRa m(t) is increasing with m, and 3) limm→∞ ∞ k=m+1 Pk(t)=0. Using Theorem 1with !(k)=d(k) and =t, the relative truncation error increases with tand, therefore, a global truncation point Rguaranteeing a relative truncation error 61at each time point tican be obtained using (10) with t=tmax = max{t1;t 2;:::;t n}¿0. This is formally established by the following theorem, where ∞ k=R+1 d(k)Pk(ti)=ETRRa R(ti)¿0 follows from d(k)¿0. Theorem 2. Let R= min m¿0: ETRRa m(tmax)¿0∧rmax ∞ k=m+1 Pk(tmax)ETRRa m(tmax)61: Then,for i=1;2;:::;n, 06∞ k=R+1 d(k)Pk(ti) ETRRa R(ti)61: In the rest of the paper we will refer to the RdeFned in Theorem 2simply as R. Let Sm()= ∞ k=m+1 Pk(): Determination of Rcould then be performed by computing Sm(tmax) for increasing values of m. The quantity Sm(tmax) could be computed using Sm(tmax)=1−m k=0 Pk(tmax), but this is numerically unstable when Sm(tmax)is1. Another, better, approach is to use tight upper bounds for Sm(tmax). Let M¿0 satisfying M+1¿t max, which implies 0 ¡t max=(M+1)¡1, and M¿m+2.
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1097 We have: Sm(tmax)= M−1 k=m+1 Pk(tmax)+ ∞ k=M Pk(tmax) ¡ M−1 k=m+1 Pk(tmax)+ ∞ k=M PM(tmax)tmax M+1k−M = M−1 k=m+1 Pk(tmax)+PM(tmax)∞ k=0 tmax M+1k =˜ Sm;M (tmax)+Se M(tmax); with ˜ Sm;M (tmax)= M−1 k=m+1 Pk(tmax); Se M(tmax)=PM(tmax)∞ k=0 tmax M+1k =PM(tmax)1 1−tmax=(M+1): The quantity ˜ Sm;M (tmax)+Se M(tmax) upper bounds Sm(tmax) with relative error 6Se M(tmax)=˜ Sm;M (tmax). Then, being &a reasonable small number, say &=10−6, we can obtain a tight upper bound for Sm(tmax) by choosing Mas the smallest nonnegative integer satisfying M+1¿t max,M¿m+2 and Se M(tmax)=˜ Sm;M (tmax)6&. Computing Sm(tmax) in that way for increasing m, starting at m=0, until ETRRa m(t)¿0 and rmaxSm(tmax)=ETRRa m(t)61would be expensive when the required Ris large, which will happen if tmax is large and 11. The computational cost of the determination of the truncation point Rwill be reduced using two improvements. The Frst improvement consists in starting the computation of Sm(tmax) at a value m0of mpotentially larger than 0. First, observe that from d(k)6rmax, we can easily obtain ETRRa m(tmax)6rmax m k=0 Pk(tmax)=rmax1−∞ k=m+1 Pk(tmax)=rmax(1 −Sm(tmax)), which implies rmaxSm(tmax)= ETRRa m(tmax)¿Sm(tmax)=(1 −Sm(tmax)) and, to determine the truncation point R, we can restrict our attention to values of mfor which Sm(tmax)=(1 −Sm(tmax)) 61, i.e. Sm(tmax)61=(1+1) and do not examine values of mfor which Sm(tmax)¿ 1=(1 + 1). Second, we have the following result. Lemma 2. Let ¿1and ',0¡'¡1. Let =−and '∗=e −0:5+1=(8)=(2(). Then,for m¿0, Sm()=∞ k=m+1 Pk()¿' if any of the following conditions hold: (a) P()¿'e2= and m6−2+2log(P()='), (b) P()¿'or e−(1=2+(1=3)e−=!) ¿',and m6−1, (c) '61−'∗and m6−3=2−√, (d) '¿1−'∗and m6−3=2−1=4−log(2((1 −')2=). Proof. See Appendix A.
1104 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 Fig. 1. Implementation of the standard randomization method for the computation of the ETRR(t) measure with control of the absolute relative error.
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1105 Fig. 2. Description of the procedure Compute ETRR invoked in Figs. 1and 3. Liand R∗ i, the computation on the variables wkof the weights wi k, the computation on the variables ˜wkof the normalized weights ˜wi k,Li6k6R∗ i, and the computation on ] ETRR(ti) of the approximate value of ETRR(ti) using R∗ i k=Lid(k)˜wi k. The computation of d(k), 0 6k6Rand the computation of approximate values for Sm(tmax) for increasing m, starting at m0, as discussed in Section 3.1 are embedded in the pseudocode determining the truncation point R. In that pseudocode, computation of ETRRa m(tmax) from m= 0 up to m=Ris done on the variable trunc sum. In the call to the function Compute ETRR(), the variable trunc sum is used Frst to compute min{R;ti} k=md(k)wi kfrom m= min{R; ti} down to m=Liand, then, to compute m k=Lid(k)wi kfrom m= min{R; ti} up to m=R∗ i. The variable may take the value 0 and that case is properly dealt with in the pseudocode. Using (1) and (4), it follows that the EARR(t) measure has a formalization in terms of Pk(t) identical to that of the ETRR(t) measure with d(k) replaced by d(k)=k l=0 d(l)=(k+ 1). This makes all the developments in the previous sections for the measure ETRR(t) to carry over the EARR(t) measure with the only diIerence that d(k) has to be replaced by d(k). This allows to adapt easily to the EARR(t) measure the implementation of the standard randomization method for the ETRR(t) measure with control of the absolute relative error. The implementation of the method for the EARR(t) measure is described in Figs. 3and 2. 4. Analysis and discussion In this section we will analyze the numerical stability of the implementations of the standard randomization method for the measures ETRR(t) and EARR(t) derived in the previous section using
1106 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 Fig. 3. Implementation of the standard randomization method for the computation of the EARR(t) measure with control of the absolute relative error. a small rewarded CTMC model of a repairable fault-tolerant system. We will also discuss the computational e$ciency of the implementations with respect to simpler alternatives.
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1107 Fig. 4. State transition diagram of the CTMC describing the behavior of the example fault-tolerant system. The rewarded CTMC model which will be used to analyze the numerical stability of the implementations corresponds to a repairable fault-tolerant computer system using the triple modular redundancy (TMR) technique including three identical processing modules and a voter. The system is up if the voter and at least two processing modules are unfailed. It is assumed that components do not fail when the system is down. Failure and repair times are assumed to have exponential distributions. The failure rate of a processing module is Mand the failure rate of the voter is V. The repair rate of a processing module is 1Mand the repair rate of the voter is 1V. Repairs are performed by a single repairman who gives preemptive priority to the voter. It is assumed that initially all processing modules and the voter are unfailed. Fig. 4depicts the state transition diagram of the resulting CTMC. The states are labeled by (PM;V), where PMis the number of unfailed processing modules and Vhas value 1 if the voter is unfailed and value 0 otherwise. The initial state of the CTMC is the state (3;1). We will consider the reward rate structure: ri=1 for i∈{(3;0);(2;0);(1;1)}; 0 for i∈{(3;1);(2;1)}: Under that reward rate structure, the ETRR(t) measure is the unavailability at time t, i.e. the probability that the system is down at time t, and the EARR(t) measure is the expected interval unavailability at time t, i.e. the expected value of the fraction of time that the system is down in the time interval [0;t]. The numerical experiments will be performed using M=10 −3h−1,V=10 −4h−1, 1M=0:5h −1, and 1V=1h−1, which yields =1h−1. The steady-state unavailability of the system is 1:2384 ×10−4. The implementations used double precision arithmetic and were run in a Sun-Blade-1000 processor. The actual errors of the numerical results given by the implementations were obtained by comparing those results with an “exact” solution of the model computed using a mathematical software package with 100 digits of accuracy, which was enough to determine with good accuracy the actual errors. Fig. 5plots the actual relative error in the numerical solution given by the implementation for the ETRR(t) measure when that implementation is run with a single target time tand an absolute relative error requirement , for several values of and t. Fig. 6plots the results obtained by the implementation for the EARR(t) measure. We can note that the actual relative error, m, is always smaller than , even for such stringent values of as 10−10. For the ETRR(t) measure and not small t =t,mis extremely small and almost independent
1108 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 Fig. 5. Actual relative errors in the implementation for the ETRR(t) measure as a function of and t(h). Fig. 6. Actual relative errors in the implementation for the EARR(t) measure as a function of and t(h). of . The explanation of that behavior is the following. For not small t, both ETRR(t) and the d(k) corresponding to non-negligible Poisson probabilities Pk(t) are almost identical to the steady-state unavailability and, then, since the implementation computes ETRR(t) by averaging d(k) with weights proportional to the Poisson probabilities Pk(t) , the actual error is basically due to round-oI errors and, therefore, extremely small and almost independent of the requested absolute relative error, . Overall, the implementations for both methods seem to have a very high numerical stability. Simpler alternatives to the implementations of the standard randomization method developed in this paper would only use the truncation point RdeFned in Theorem 2with 1=,being the requested relative error, and would estimate the ETRR(t) measure at each time point tiusing R k=0 d(k)Pk(ti) and the EARR(t) measure at each time point tiusing R k=0 d(k)Pk(ti). This would require to
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1109 Table 2 Truncation point Rused by the implementations developed in the paper and the truncation point Rin simpler alternatives for the example when the implementations for both measures are run with a single target time tand =10 −4 ETRR(t)EARR(t) t(h) RR RR 0.01 4444 0.16666 111111111 10 33 32 33 32 100 163 161 163 161 1031,189 1,181 1,189 1,181 10410,587 10,562 10,587 10,562 105101,843 101,768 101,843 101,768 1061,005,817 1,005,580 1,005,817 1,005,580 10710,018,383 10,017,634 10,018,383 10,017,634 compute “exact” Poisson probabilities. Call Rthe truncation point Rin those alternatives. Then, since Rdecreases with 1,Rwould be smaller than the truncation point Rof the implementations developed in the paper. This would allow a reduction on the required number of vector-matrix multiplications (3), which for large rewarded CTMC models is an important component of the computational cost of the method. However, being Ra smooth function of 1, the reduction will not be, in general, signiFcant. This is illustrated in Table 2which compares Rand Rfor the example when the implementations for both measures are run with a single target time tand =10 −4.At the price of increasing slightly the truncation point Rand the required number of vector-matrix multiplications, the implementations developed in this paper have computational advantages. One of them is a better numerical stability, since available approaches to compute “exact” Poisson probabilities, e.g. [6, pp. 1028–1029], are less stable numerically than the use of normalized weights. The implementations developed in this paper have also computational advantages when tmax is large and the number of time points nat which the measure has to be computed is large. In that case, the improvements developed in Section 3.1 to determine the truncation point Rwhich are incorporated in the implementations reduce signiFcantly the number of “exact” Poisson probabilities which have to be computed and the truncation points Liand R∗ iintroduced in Section 3.2 make the number of weights and normalized weights which have to be computed for each time point tito be signiFcantly smaller than the number of “exact” Poisson probabilities which would have to be computed in the simpler alternative implementations. 5. Conclusions Based on a theoretical result regarding the dependence on the parameter of the Poisson distribution of the relative error introduced when a weighted sum of Poisson probabilities is truncated by the right, we have developed implementations of the standard randomization method for the computation of the “expected transient reward rate” and the “expected averaged reward rate” measures over rewarded
1110 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 CTMCs with control of the relative error. The methods seem to exhibit an excellent numerical stability and have computational advantages over simpler alternatives. The implementations developed in the paper are sometimes preferable from a practical point of view over the standard randomization method as it is usually implemented, i.e. with control of the absolute truncation error. Appendix A. Proof of Lemma 2.Let m∗=−2+2log(P()='). Then, for 0 6m6m∗, which implies m+16m∗+ 1, we have Sm()= ∞ k=m+1 Pk()¿P m∗+1():(A.1) From [19, Proposition 5] it follows that, for any strictly positive integer i, P+i()¿P() exp−(i+1) 2 2:(A.2) Noting that P()¿'e2= implies m∗¿and, therefore, m∗¿, we have m∗+1−¿1. Then, using (A.1), using (A.2) with i=m∗+1−, noting that e−x2is decreasing, and substituting m∗, we obtain, for 0 6m6m∗: Sm()¿P m∗+1()¿P() exp−(m∗+2−)2 2 ¿P() exp−(m∗+2−)2 2 =P() exp−2log(P()=') 2=': (b) Assume Frst P()¿'. Then, for 0 6m6−1, which implies ¿m+1, Sm()= ∞ k=m+1 Pk()¿P ()¿': Assume next e−(1=2+(1=3)e−=!) ¿'. Let Zbe a Poisson random variable with mean , and let 4such that 1=2=P[Z¡]+4P[Z=]. Using a result given in [20, p. 780], which depends on ¿1, 4¿1=3 and, therefore, 1=2−P[Z¡] P[Z=]=4¿1 3:(A.3) From (A.3) we can obtain 1 −P[Z¡]¿1=2+(1=3)P[Z=] and P[Z¿]¿1 2+1 3P[Z=]:(A.4)
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1111 Moreover, since 0 6¡1 and ¿1, for k¿0, Pk()=e −k k!=e −e−(+)k k!¿e−e−k k!=e −P[Z=k]: Then, for 0 6m6−1, which implies ¿m+ 1, using (A.4), and using e−(1=2+(1=3)e− =!) ¿', Sm()= ∞ k=m+1 Pk()¿∞ k= Pk()¿e−∞ k= P[Z=k]=e −P[Z¿] ¿e−1 2+1 3P[Z=]=e −1 2+1 3e− !¿': (c) Let U 5(x)=∞ x6(u)du, with 6(u)=(1=√2()e−u2=2. It is proved in [21, p. 13] that m k=0 Pk()6 e1=(8)U 5−m−3=2 √;¿1;06m6−2: From [7, Formula 26.2.12], U 5(x)=6(x) x−6(x) x3−15 ∞ x 6(u) u6du¡6(x) x;x¿0: Then, since m6−2 implies −m¿2 and (−m−3=2)=√¿0, m k=0 Pk()¡ 2(e1=(8)exp −1 2−m−3=2 √2 ×√ −m−3=2=Tm();¿1;06m6−2:(A.5) Also, extending the deFnition of Tm()in(A.5)tomreal, we have '∗=T−3=2−√():(A.6) Let f(x)=e −x2=2=x. Since f(x) is decreasing for x¿0, Tm()= 2(e1=(8)f−m−3=2 √ 6 2(e1=(8)f−m−3=2 √ =Tm();¿1;06m6m6−2:(A.7) Since ¿1 implies −3=2−√6−2, for 0 6m6−3=2−√we can use (A.5) and use (A.7) with m=−3=2−√, yielding m k=0 Pk()¡T m()6T−3=2−√():
1112 V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 Then, using (A.6) and recalling that, by assumption, '∗61−', m k=0 Pk()¡' ∗61−': Finally, Sm()= ∞ k=m+1 Pk()=1− m k=0 Pk()¿1−(1 −')=': (d) Let g(x)=e −x2=2. Since g(x) is decreasing for x¿0, retaking the deFnition of Tm() given in (A.5) extended to mreal, gm()=Tm()−m−3=2 √= 2(e1=(8)g−m−3=2 √ 6 2(e1=(8)g−m−3=2 √=gm();¿1;06m6m6−2:(A.8) A solution m(')ofgm()=1−'is m(')=−3 2−1 4−log2((1 −')2 ; which is real for '¿1−'∗. Trivially, since m(') is decreasing on 'for 0 ¡'¡1, and, by assumption, '¿1−'∗,m(')6m(1 −'∗). It is also easy to check that m(1 −'∗)=−3=2−√. Then, using ¿1, for m6−3=2−1=4−log(2((1 −')2=)=m('), m6m(')6−3 2−√¡−2;(A.9) from which we get √ −m−3=261: Using, then, the deFnition of gm()in(A.8) we obtain Tm()= √ −m−3=2gm()6gm():(A.10) Moreover, since (A.9)m(')¡−2 and ¿1, for 0 6m6m('), we have (A.8) gm()6gm(')()=1−': (A.11) Finally, since (A.9)m¡−2, we can invoke (A.5), which together with (A.10) and (A.11) yields m k=0 Pk()¡T m()61−';
V. Su˜ n* e, J.A. Carrasco / Computers & Operations Research 32 (2005) 1089–1114 1113 and, therefore, Sm()= ∞ k=m+1 Pk()=1− m k=0 Pk()¿1−(1 −')=': Proof of Proposition 2.Assume Li¿0. Taking into account Lemma 4, part a), that, because of Lemma 4, part a), min{R;ti} k=LiPk(ti)¿0, that, as noted in Section 3.2,Pk(ti) is increasing on kfor 0 6k6Li−1, that 0 6d(k)6dmax,Li6k6R∗ i, and the way in which Liis selected (13): Li−1 k=0 Pk(ti) min{R;ti} k=LiPk(ti) 6LiPLi−1(ti) min{R;ti} k=LiPk(ti) =LidmaxPLi−1(ti) dmax min{R;ti} k=LiPk(ti) 6LidmaxPLi−1(ti) min{R;ti} k=Lid(k)Pk(ti) 62 2; implying because 0 ¡min{R;ti} k=LiPk(ti)¡1 Li−1 k=0 Pk(ti)62 2:(A.12) The previous inequality can be trivially extended to the case Li=0. Assume R∗ i¡R. Taking into account Lemma 5, that, because of Lemma 5,R∗ i k=LiPk(ti)¿0, that, as noted in Section 3.2,Pk(ti) is decreasing on kfor k¿R∗ i+ 1, that 0 6d(k)6dmax, Li6k6R∗ i, and the way in which R∗ iis selected (14): R k=R∗ i+1 Pk(ti) R∗ i k=LiPk(ti) 6(R−R∗ i)PR∗ i+1(ti) R∗ i k=LiPk(ti)=(R−R∗ i)dmaxPR∗ i+1(ti) dmax R∗ i k=LiPk(ti) 6(R−R∗ i)dmaxPR∗ i+1(ti) R∗ i k=Lid(k)Pk(ti) 62 2; implying again, because 0 ¡R∗ i k=LiPk(ti)¡1 R k=R∗ i+1 Pk(ti)62 2:(A.13) The previous inequality can be trivially extended to the case R∗ i=R. Finally, since ∞ k=R+1 Pk(t) increases with t, because it is the probability that the number of arrivals in a Poisson process with arrival rate in the time interval [0;t]is¿R+ 1, taking into account the way in which Ris selected (10) and using 0 ¡ ETRRa R(tmax)6rmax, ∞ k=R+1 Pk(ti)6∞ k=R+1 Pk(tmax)61 rmax ETRRa R(tmax)61:(A.14) Finally, combining (A.12), (A.13) and(A.14) and using ∞ k=0 Pk(t)=1: R∗ i k=Li Pk(ti)=1− Li−1 k=0 Pk(ti)− R k=R∗ i+1 Pk(ti)−∞ k=R+1 Pk(ti)¿1−(1+2):