Importance sampling for Markovian tandem queues using subsolutions: exploring the possibilities

We consider importance sampling simulation for estimating the probability of reaching large total number of customers in an M | M | 1 tandem queue, during a busy cycle of the system. Our main result is a procedure for obtaining a family of asymptotically efficient changes of measure based on subsolutions. We explicitly show these families for two-node tandem queues and we find that there exist more asymptotically efficient changes of measure based on subsolutions than currently available in literature.


Introduction
In this paper, we explore possibilities for importance sampling in a Markovian tandem queue. With importance sampling, the underlying probability distributions are changed to speed up a simulation. This change of the probability distributions is also called a change of measure. When exploring the possibilities for a change of measure, we focus on the probability of total buffer overflow during a busy cycle of an MjMj1 tandem queue.
This problem was first studied by Parekh and Walrand, 1 where a state-independent change of measure is suggested based on heuristics. In that paper, the authors note that their change of measure performs poorly in practice in some cases. Afterwards, in Glasserman and Kou, 2 necessary and sufficient conditions have been determined for d-node MjMj1 tandem queues in order for the same change of measure to be asymptotically efficient (which means that the relative error of the estimator grows less than exponentially fast). These conditions have been extended in the work of de Boer 3 for the case d = 2 and, in that paper, it also has been shown that the change of measure from Parekh and Walrand 1 is the only state-independent change of measure that may be asymptotically efficient when d = 2. The conclusion of these two papers is that the stateindependent change of measure from Parekh and Walrand 1 for the d-node MjMj1 tandem queue is not asymptotically efficient for all input parameters.
To resolve this issue, in the work of Dupuis et al., 4 a state-dependent change of measure has been developed for the two-node MjMj1 tandem queue such that this change of measure is asymptotically efficient for all input parameters. This change of measure is derived using socalled subsolutions, which have been introduced (in the context of importance sampling) by Dupuis and Wang. 5 Afterward, the work from Dupuis et al. 4 has been extended to Jackson networks in the work of Dupuis and Wang, 6 and to non-Markovian tandem queues in the work of Buijsrogge et al. 7 For more details on subsolutions, see also Budhiraja and Dupuis's, 8 Chapters 14 and 15. Another approach has been developed by Blanchet,9 where an algorithm is presented that gives bounded relative error. This approach uses the time reversed process of Jackson networks. While the time reversed process is known for Jackson networks, this is not the case in other contexts, e.g., for non-Markovian tandem queues. Thus, even though a more efficient simulation approach might exist for the model in the current paper, the approach that we use here is generalizable to processes where the time reversed process is unknown, see, for example, Buijsrogge et al. 7 Therefore, we believe that the insights obtained in this paper can be useful for those types of processes.
In addition, in the work of Sezer, 10 an approximation method has been developed for the two-node MjMj1 tandem queue using subsolutions. While the approximation of this probability could be very interesting in some cases, an unbiased simulation might be preferred in other cases or could be used to verify the results. Thus, in this paper, we focus on simulation.
Even though the problem that we consider has been studied previously in [1][2][3][4]6,9,11 and asymptotic efficiency of changes of measure based on subsolutions has been shown in some of these works, 4,6,11 still there are ''only'' three such changes of measure known that have been proven to be asymptotically efficient for the two-node MjMj1 tandem queue: two when queue 2 is the bottleneck queue (i.e. m 2 4 m 1 , where m i is the service rate of queue i), and one when queue 1 is the bottleneck queue (m 1 4 m 2 ), see Dupuis et al. 4 and Dupuis and Wang. 6 However, in the work of Buijsrogge et al., 7 we have shown that there are several possible changes of measure based on subsolutions yielding an asymptotically efficient estimator for GIjGIj1 tandem queues with bounded support. Since the exponential distribution clearly does not have bounded support, this raises some thoughts on whether there exist more asymptotically efficient changes of measure based on subsolutions for the two-node MjMj1 tandem queue.
In this paper, we give sufficient conditions for an asymptotically efficient change of measure based on subsolutions for the MjMj1 tandem queue. These conditions follow naturally when using the same method as in the works of Dupuis et al. 4 and Dupuis and Wang, 6 but have never been stated as such. We show how to find changes of measure satisfying these conditions and we believe the same could be done for other models. There is no numerical section as we do not aim for more efficient simulation schemes for the particular tandem model at hand (for which we could simply follow Blanchet 9 ), but we aim for providing a procedure to obtain possible schemes based on subsolutions that surely lead to an asymptotically efficient estimator. For more complex models, we do have numerical results, see Buijsrogge et al. 7 and Buijsrogge. 12 Moreover, in this paper, we consider both m 1 4 m 2 and m 2 4 m 1 . Even though this does not affect the probability of interest, since both queues are interchangeable, 13 it is particularly interesting to see the possibilities for a change of measure when m 1 4 m 2 . This is in line with Buijsrogge et al., 7 where we have considered both queue 1 and queue 2 being the bottleneck queue (for a two-node MjMj1 tandem queue this is equivalent to m 1 4 m 2 and m 2 4 m 1 , respectively). This is in slight contrast with most existing literature on importance sampling for two-node Markovian tandem queues, see for example, 2,3,4,11 where only the second queue being the bottleneck queue has been discussed (i.e., m 2 4 m 1 ).
The contributions of this paper are two-fold. After summarizing the subsolution method for importance sampling and stating the results from Dupuis et al. 4 and Dupuis and Wang 6 in Section 2, our first contribution is in Section 3, where we state conditions for a change of measure for the d-node MjMj1 tandem queue based on subsolutions to give an asymptotically efficient estimator, and we prove that these conditions are sufficient for d = 2. The other contribution, in Section 4, is that we provide a whole family of changes of measure for the two-node MjMj1 tandem queue that satisfy these conditions and hence result in an asymptotically efficient estimator. We conclude this paper in Section 5.

The model
In this paper, we consider a d-node MjMj1 tandem queue, with arrival rate l and service rates m 1 ,..., m d for queues 1, . . . , d, respectively, and we are interested in estimating the probability that the total number of customers in the system reaches some high level N during a busy cycle of the system. Since the case d = 2 has been studied mostly in literature, we mainly consider this case. However, it seems likely that all of the results can also be extended to d . 2 and so we will briefly comment on these extensions when necessary. Thus, we now let d = 2.
We consider the underlying embedded discrete time Markov chain and we assume without loss of generality l + m 1 + m 2 = 1. Furthermore, we assume that l \ minfm 1 , m 2 g, so that the system is stable. As in previous works, 4,6,11 we let the state description be the number of customers in each queue, denoted by Z i = (Z 1, i , Z 2, i ), where Z j, i is the number of customers in queue j after i transitions.
We let v k denote the possible transitions and Y(v k ) their corresponding probabilities, i.e., v 0 = (1, 0) corresponds to an arrival to the first queue, and has probability Y(v 0 ) = l. Similarly, we have v 1 = ( À 1, 1) and v 2 = (0, À 1), Y(v 1 ) = m 1 and Y(v 2 ) = m 2 . The transitions v 1 and v 2 can only occur when Z 1, i . 0 and Z 2, i . 0, respectively. If Z j, i = 0 for some queue j, thus, there are no customers in queue j, we add a self-loop transition with probability m j to make sure that the sum of all rates equals 1 (so that all rates are probabilities).
As in the work of Dupuis et al., 4 we let X i = 1 N Z i be the scaled state description, which has the advantage that its elements are in [0,1] and therefore do not increase as N increases. This allows us to make the following definitions: and this is sketched in Figure 1.
Using these definitions, we can define the first time that the process hits level N in a busy cycle as t N = inffi . 0 : X i 2 d e , X k 6 2 d 0 8k = 1, . . . , i À 1g, and we set t N = ' if the process hits d 0 before d e . Therefore, our probability of interest can be written as It is known that the asymptotic decay rate of this probability is given by where g j =À log (l=m j ), see Glasserman and Kou 2 and Buijsrogge et al. 14 The queue for which we have g j = g is called the bottleneck queue. For the MjMj1 tandem queue, the bottleneck queue is equivalent to the queue with the largest server utilization. In most papers on similar topics, see, for example, Glasserman and Kou, 2 de Boer, 3 Dupuis et al. 4 and de Boer and Scheinhardt, 11 it is assumed that m 2 4 m 1 , as for the probability of interest the queues are interchangeable. 13 In this paper, we both consider m 2 4 m 1 and m 1 4 m 2 .

Importance sampling simulation
To estimate our probability of interest using simulation, we use importance sampling. In importance sampling, we perform our simulation under some new measure Q. Under this new measure, we let Y(v k jx) denote the probability for transition v k given that we are in state x. While doing the simulation, we keep track of the likelihood ratio L(P) of a path P = (X i , i = 0, . . . , t N ): where to include the self-loop transition when one of the queues is empty. Let I(P) = 1ft N \ 'g indicate whether we have reached our event of interest during a busy cycle of the system or not.
Then, under the new measure Q, we have where E Q denotes the expectation under the new measure Q. Thus, L(P)I(P) is the estimator for our probability of interest. As in the work of Dupuis et al., 4 we construct a subsolution W (x) -formally defined in Section 2.3 -and we use this function to specify a change of measure in the following way: where DW (x) denotes the gradient of W (x) with respect to x, and where If we compare Equation (4) with the corresponding notation in the works of Dupuis et al. 4 and de Boer and Scheinhardt, 11 a factor 2 is missing. However, we will also scale the function W (x) accordingly, so that the change of measure remains the same.
As we are interested in finding a change of measure that gives an asymptotically efficient estimator for p N , we need lim sup  In previous works, 4,6,11 a subsolution is constructed (assuming m 2 4 m 1 in the works of Dupuis et al. 4 and de Boer and Scheinhardt 11 ) and asymptotic efficiency is proven for a specific choice of W (x). In this paper, we provide conditions on W (x) and prove that under these conditions, we get an asymptotically efficient estimator. Afterwards, we study the possibilities of the function W (x), and the corresponding change of measure in Equation (3).

Subsolution approach
In both Dupuis et al. 4 and Dupuis and Wang, 6 the change of measure for (two-node) tandem queues (and Jackson networks) has been studied. In those papers, the change of measure has been determined using subsolutions. We first briefly recap the ideas presented in those papers and we start with a formal definition of a subsolution, see also Dupuis et al. 4 Definition 2.1. A function W (x) is a classical subsolution if: In Dupuis et al., 4 there is an additional condition in this definition that needs to hold at the boundaries of the state space. Instead, we include the boundaries of the state space in H(x, DW (x)), i.e., H(x, DW (x)) differs along the boundaries (which is similar as in the work of Dupuis and Wang 6 ).
It is known from Glasserman and Kou 2 and de Boer 3 that for an asymptotically efficient change of measure it is not possible that DW (x) is constant throughout the whole state space as in the work of Parekh and Walrand, 1 and in the work of Dupuis et al. 4 this is ''confirmed,'' since the change of measure yielding an asymptotically efficient estimator differs from the change of measure in the work of Parekh and Walrand 1 along one of the boundaries of the state space and near the origin. Thus, to find an asymptotically efficient change of measure, we determine several (constant) changes of measure for various regions of the state space, in particular along the boundaries of the state space, and combine these so that we have a change of measure for the whole state space.
To determine such a change of measure that differs along various parts of the state space, in the work of Dupuis et al., 4 there are multiple -say r -affine functions W d k (x), k = 1, . . . , r, considered such that for each of these functions H(x, DW d k (x)) ø 0 for some part of the state space, so that all r functions cover the whole state space, and for at least one of these functions, we have W d k (x) 4 0 for x 2 d e . All these functions have the following form: where a k = (a k, 1 , a k, 2 ), which we assume to be finite, specifies the gradient of W d k (x) and both c k and d k . 0 are constants. Combining these functions to W d (x) by taking the minimum of these functions for all x, then results in a piecewise affine function, that unfortunately is not continuously differentiable. That is, we have, To satisfy the continuous differentiability, which is the first requirement for a function to be a classical subsolution, the functions W d k (x) are combined into a (continuous) function W e, d (x) by a similar mollification procedure as in the works of Dupuis et al. 4 and Dupuis and Wang, 5 where in this case we obtain such that W e, d (x) converges to W d (x) when e ! 0. Throughout this paper, we make the following assumptions on e and d (that depend on N), which also have been made in previous works. 4,6,11 For convenience, we do not write the explicit dependence of e and d on N.
Assumption 2.1. We choose e and d dependent on N, such that The gradient of Equation (6) is then used as change of measure in Equation (3). It can be expressed as The functions r k (x) are weight factors for the ''influence'' of each function W d k (x), and so of each different ''regional'' change of measure, in the final change of measure. They can also be used to define a change of measure slightly different than in Equation (3) as follows, see also Dupuis et al. 4 and de Boer and Scheinhardt 11 , which we will refer to in this paper later on. In fact, we will also show that asymptotic efficiency for the change of measure in Equation (3) implies asymptotic efficiency for the change of measure in Equation (8), similar as in the work of de Boer and Scheinhardt. 11 We note that, from an implementation perspective, the change of measure in Equation (8) and the function W d (x) is illustrated in Figure 2. We remark that also in the functions in Equation (9), we scaled the results from Dupuis et al. 4 (except for the constant in front of d) by a factor 1/2, but this does not influence the resulting changes of measure in Equations (3) and (8).

2.4.2.
Changes of measure from Dupuis and Wang. 6 . In the work of Dupuis and Wang, 6 the work of Dupuis et al. 4 is extended to Jackson networks and hence in the work of Dupuis and Wang, 6 all queues being the bottleneck queue are considered, as in this paper. Not only the probability that the total number of customers in the system reaches some high level N during a busy cycle of the system is considered in the work of Dupuis and Wang, 6 but the authors also consider buffer overflow in a single queue or in several queues at the same time. If we consider a two-node tandem queue with queue 2 being the bottleneck queue, we find that the following four functions are used in Dupuis and Wang 6 : and W d (x) is illustrated in Figure 3.
In the work of Dupuis and Wang, 6 the authors do not explicitly mention by which constant d is multiplied, (9). Queue 2 is the bottleneck queue (so g 2 4 g 1 ). though implicit requirements with a proof for asymptotic efficiency are given. For simplicity, we use kd throughout this paper. It turns out that these values are sufficient for asymptotic efficiency, as we show later in this paper, but we will also see that they are by no means unique.
Since there is no limitation to queue 2 being the bottleneck queue in Dupuis and Wang, 6 we also present the result from Dupuis and Wang 6 when queue 1 is the bottleneck queue. The result from Dupuis and Wang, 6 when queue 1 is the bottleneck queue, is different compared with when queue 2 is the bottleneck queue, even though for the probability of interest both queues are interchangeable. When queue 1 is the bottleneck queue, the following three functions are derived by Dupuis and Wang 6 : and the resulting function W d (x) is illustrated in Figure 4. These subsolutions are very similar to the ones in Equation (10), but with one function less. Also here, the multiplication factor of d is not explicitly mentioned, but we use the values above which are sufficient for asymptotic efficiency, as we will show later in this paper.

2.4.3.
Comparison of the existing changes of measure. In this section, we briefly comment on the similarities and differences of the changes of measure from Dupuis et al. 4 and Dupuis and Wang. 6 We will do so by comparing the functions W d (x) for all different cases (see . When queue 2 is the bottleneck queue, we see from Figures 2 and 3 that along the x 2 -axis the function W d (x) is roughly the same. The only difference is that in Figure 3 the function, along the x 2 -axis, is d lower. In particular, this part of the state space covers the most likely path. In all other parts of the state space, the function W d (x) in Figure 3 is slightly steeper than in Figure 2 since g 2 \ g 1 . These observations suggest that any change of measure based on some function W d (x) that somehow lies ''in between'' the functions in Figures 2 and 3 is also asymptotically efficient. In Section 4, we show that this is indeed the case.
When queue 1 is the bottleneck queue, there is not much to compare. However, since there are already two possibilities for the change of measure based on subsolutions to be asymptotically efficient, and even more to expect, when queue 2 is the bottleneck queue, also the case when queue 1 is the bottleneck queue is studied in Section 4.

Sufficient conditions for asymptotic efficiency
Similar to Dupuis et al. 4 and Dupuis and Wang, 6 the construction of the changes of measure in this paper is based on finding appropriate subsolutions W e, d (x). We start with a general proof for asymptotic efficiency for a change of measure based on the subsolution approach and the mollification procedure that is explained in Section 2.3. In our theorem, we provide sufficient conditions for the subsolution yielding an asymptotic efficient change of measure, which we use later for the derivation of the possibilities for the change of measure. Afterwards, we discuss some general observations with respect to some of the conditions in our theorem.

Main result
Theorem 3.1. Consider a two-node MjMj1 tandem queue. Let W d k (x) be affine functions for all k = 1, . . . , r, as in Equation (5), and let the classical subsolution W e, d (x) be constructed using Equation (6). Then, using the gradient of the function W e, d (x) as change of measure in Equations (3)   Proof. We start by showing that under the above conditions the change of measure in Equation (3) is asymptotically efficient, after which it follows that also the change of measure in Equation (8) is asymptotically efficient using a similar argument as in Theorem 2 from de Boer and Scheinhardt. 11 From Equations (2) and (3), it follows that the likelihood ratio of a path L(P) is given by We find, using Equation (7), that for all states x we have due to concavity of H(x, DW e, d (x)) in its second argument, see Proposition 3.2 in the work of Dupuis et al. 4 Combining the two expressions above, we arrive at where the last inequality follows from Condition 1. Similar to Lemma 2 in the work of de Boer and Scheinhardt, 11 also when using r regions, we can obtain the following bound. The idea of this Lemma in de Boer and Scheinhardt 11 is to replace the summation in Equation (13) by W e, d (X t N ) À W e, d (0) and to give an upper bound on the error that is introduced. Thus, since a k is finite by construction, and hence ja k j 4 c for some 0 Next, we follow similar steps as in Theorem 1 of the same paper. By combining Equations (13) and (14) we have where the second inequality follows from Conditions 2 and 3 when X t N 2 d e , and thus, we find, as in de Boer and Scheinhardt, 11 To conclude the proof, we need Lemma 3 from de Boer and Scheinhardt, 11 which states that for any sequence u N ø 0 such that lim Thus, taking limits in Equation (15) gives where the last equation follows using Equation (1), ! 0 when N ! ', we can apply Equation (16) to the fourth term of Equation (17). This concludes the proof for the change of measure in Equation (3).
For the change of measure in Equation (8), we note that similar to de Boer and Scheinhardt, 11 we have where the inequality follows by concavity of the logarithm and the last equality follows by definition of DW e, d (x), see Equation (7). Thus, we have the same bound as in Equation (12) and so we also find that the change of measure in Equation (8) is asymptotically efficient. h Remark 3.1. In Equation (17), we see that we end up with some term g(N ) À h(N) À g, which goes to Àg as N ! '. This term arises from bounding W e, d (X t N ) À W e, d (0) on a path that leads to the overflow level. However, we cannot have W e, d (X t N ) À W e, d (0) \ À g when N ! ', since this would contradict with Jensen's inequality: That is, it is impossible to obtain a tighter bound. As a result, we find that for an asymptotically efficient change of measure we need W e, d (X t N ) À W e, d (0) ! Àg when N ! ' on at least one path that leads to reaching the overflow level, e.g., the most likely path, see also Dupuis and Wang. 6 Thus, on such a path, we need lim In the sequel, we use a slightly stronger condition than Condition 1 of Theorem 3.1, namely that for each k separately, we have Remark 3.3. It seems likely that Theorem 3.1 can also be extended to a d-node MjMj1 tandem queue (and Jackson networks), with the same sufficient conditions as in the current statement for d = 2. To do so, observe that in the proof of Theorem 3.1, Lemma 3 from de Boer and Scheinhardt 11 is the only part that restricts to d = 2 (and hence to tandem queues). Thus, one could either extend this result to d . 2 (and Jackson networks), or use similar techniques as in the works of Dupuis et al. 4 and Dupuis and Wang 6 in order to show that the theorem holds in a more general setting.

General observations
Now that we have shown under which conditions we obtain an asymptotically efficient change of measure based on subsolutions, it remains to find a k , c k and d k for all k = 1, . . . , r such that Conditions 1, 2, and 3 of Theorem 3.1 are satisfied.
In this section, we make some general observations with respect to Conditions 1 and 3 of Theorem 3.1 when considering a two-node MjMj1 tandem queue, that are used later to construct the possibilities for the change of measure. These observations are independent of the bottleneck queue.
3.2.1. Observations with respect to Condition 1 of Theorem 3.1. We recall that the first condition is In this section, we state some observations with respect to H(x, a) for some general a, independent of k, after which we present some observations with respect to r k (x).
Proof. For the first statement, let a 2 =À g 2 and x j . 0 for j = 1, 2.
Then, Equation where we let x = e Àa 1 . Using elementary calculus we find that a 1 2 ½À max j g j , À min j g j . The other statements follow similarly. h In the following lemma, we consider Equation (19) at one of the boundaries given a certain choice for either a 1 or a 2 . The first equation considers x 1 . 0, x 2 = 0 and a 1 =À g 1 , and the second equation considers x 1 = 0, x 2 . 0 and a 2 =À g 2 . m 1 + le Àa 2 + m 2 4 1 iff a 2 ø 0, and le Àa 1 + m 1 + l 4 1 iff a 1 øÀ g 2 .
Proof. The statements follow directly by elementary calculus, combined with l + m 1 + m 2 = 1 and l . 0. h We conclude this section with a remark on r k (x). Using Equation (7) we find that for any ', where the inequality follows trivially.

Observations with respect to Condition 3 of Theorem
, such that lim N !' h k (N ) = 0 for some k, see Equation (6). By Equation (5), we have W d k (0) = c k À d k d. Therefore, we find c k = g for all k, as a sufficient condition to satisfy Condition 3 since d ! 0 when N ! ' by Assumption 2.1.
As a result of the observations above, to construct the possibilities for the change of measure, it remains to find a k and d k for all k = 1, . . . , r that satisfy Conditions 1, 2, and 3 of Theorem 3.1. This is the topic of the next section.

Construction of the subsolution W e, d (x)
In this section, we construct possible subsolutions, based on the approach mentioned in Section 2.3, that satisfy the conditions in Theorem 3.1 and thus yield an asymptotically efficient estimator. It may not be clear at first sight that the method below results in subsolutions that satisfy all these conditions, since the construction is partly based on intuition. However, we conclude all sections by showing that the conditions are indeed satisfied.
For the two-node MjMj1 tandem queue, we consider both possibilities for the bottleneck queue, i.e., we consider both queue 1 and queue 2 as bottleneck queue. In addition, we focus on a maximum of four different regions. More  regions may or may not be possible. However, this is undesirable from a practical point of view and does not contribute to an easier implementation of the change of measure.
We start using three regions and queue 2 being the bottleneck queue, since this case has been studied most in literature. Afterwards, we consider three regions and queue 1 the bottleneck queue. We conclude this section with four regions, for which we again consider both queue 2 and queue 1 as the bottleneck, respectively. For brevity, when considering four regions, we will only state the result of the construction (which is similar to the construction when having three regions) and show that indeed the conditions in Theorem 3.1 are satisfied.

Three regions and queue 2 bottleneck
In this section, our starting point is to consider three functions W d k (x), k = 1, 2, 3, and we let queue 2 be the bottleneck queue, thus g 2 4 g 1 . As the boundaries turned out to be crucial in designing an asymptotically efficient change of measure, we consider the following regions: (i) x 2 . 0, which also covers x 1 = 0; (ii) x 1 . 0, which also covers x 2 = 0; and (iii) x 1 ø 0 and x 2 ø 0, so that we have covered the whole state space. All regions overlap in the sense that they all cover the case in which both x 1 . 0 and x 2 . 0. However, by construction of the continuously differentiable subsolution, in that part of the state space the function W d k (x) of only one of the three regions will be used since we use the minimum of the functions W d k (x). Clearly, the third region covers the whole state space, but it is important to note that there is no non-trivial solution that satisfies Condition 1 from Theorem 3.1 for the whole state space. The most important part of region three is that it covers x 1 = x 2 = 0.
It turns out that the zero change of measure can be used when both x 1 = x 2 = 0. This may seem strange at first sight, since simulating under the zero change of measure is the same as simulating the original system. However, using this as a starting point for this part of the state space turns out to be a good choice (which is the same choice as in the works of Dupuis et al. 4 and Dupuis and Wang 6 ). In particular, the condition in Equation (19) for this part of the state space is equivalent to a 3, 1 ø 0. The zero change of measure gives us a 3 = (a 3, 1 , a 3, 2 ) = (0, 0) and hence satisfies the condition mentioned.
The ordering of the regions that we assign can be found in Table 1, the reasons for this ordering will become clear later in this section.
4.1.1. Finding a 1 . To find a 1 , we start with Condition 2, which is W e, d (x) 4 g(N ) for all x 2 d e , where lim N !' g(N ) = 0. In Remark 3.1, it is noted that for the most likely path equality should hold. Trivially, we have W e, d (x) 4 W d 1 (x), with equality on the most likely path when N ! ', since W d 1 (x) covers the most likely path. As a result, taking into account Section 3.2.2 and recalling that a k = (a k, 1 , a k, 2 ) for all k, we must have a 1, 2 =À g 2 .
Using Condition 1, we can now determine a 1, 1 . We only consider the term of the summation from this condition that involves k = 1, i.e., r 1 (x)H(x, a 1 ) should be nonnegative for large enough N. As W d 1 (x) is not designed for x 2 = 0, the intuition is that the weight factor r 1 (x) tends to 0 for large enough N for all states x such that x 2 = 0, see also Equation (28). Thus, using Equation (18) for a = a 1 , we see that for all x such that x 2 . 0, we can only satisfy Condition 1, if we have le where we used a 1, 2 =À g 2 . As a result of Lemma 3.2, the first bullet, we find that Equation (21) is satisfied when a 1, 1 2 ½Àg 1 , À g 2 . Using Lemma 3.3, we find that a 1, 1 øÀ g 2 is necessary to satisfy Equation (22), and hence, we must have a 1, 1 =À g 2 . Therefore, we need a 1 = ( À g 2 , À g 2 ), to get an asymptotically efficient change of measure based on Theorem 3.1. For future reference, we remark that as a result of this condition on a 1 we find, using Equation (18), since l + m 1 + m 2 = 1.

Finding a 2 .
Using the underlying idea of the construction of the subsolution -i.e., the idea to construct several functions, each for different parts of the state space, that are combined through mollification to obtain a classical subsolution -we determine conditions on a 2 and d k for k = 1, 2, 3. For some parts of the state space, it is known which function W d k (x) has to be the minimum of the three functions or cannot be the minimum of the three. For example, when x 1 = 0, it follows that W d 2 (x) cannot Table 1. Overview of proposed regions for the case r = 3.
Region k be the minimum function, since it is designed for x 1 . 0. As a result, we find that for some parts of the state space, some weight factors r k (x) must tend to 0 as N ! '.
To start with, we consider the origin of the state space, i.e., x 1 = x 2 = 0. Here, we want W d 3 (x) to be the minimum function, since this is the only function that is designed for this part of the state space. Thus, we need both Trivially, these inequalities result in Second, we consider the boundary x 2 = 0 (and so x 1 . 0). At this part of the state space, we want W d 2 (x) to be the minimum function. Thus, we want to have both for all x such that x 2 = 0. Clearly, the first inequality holds for all x 1 whenever As a result, we immediately have, for all x such that x 2 = 0, using Equation (20) for k = 1 and ' = 2, for all a 2, 1 4À g 2 and since d 1 \ d 2 , the right-hand side tends to 0 as N tends to infinity. This intuitively implies that the weight factor r 1 (x) for x such that x 2 = 0 tends to 0 when N ! ', as suggested in Section 4.1.1, and thus, the change of measure that is designed for x 2 . 0 hardly has any influence when x 2 = 0. The second inequality, Equation (26), is satisfied for all x 1 . (d 3 Àd 2 )d Àa 2, 1 , which is positive because a 2, 1 4À g 2 \ 0, see Equation (27), and d 2 \ d 3 . Thus, W d 2 (x) is the minimum function for all x 1 . (d 3 Àd 2 )d Àa 2, 1 , and note that the righthand side of this inequality tends to zero as N ! '. For all other (very small) x 1 , the function W d 3 (x) is the minimum function. It turns out in the sequel that this is not a problem for the resulting change of measure to be asymptotically efficient, since the function W d 3 (x) can be used throughout the whole state space. More importantly, W d 1 (x) is not the minimum function for all x such that x 1 . 0 and x 2 = 0 whenever Equation (27) holds.
Combining all conditions on d k , see Equations (24) and (27), we find d 1 \ d 2 \ d 3 . Here, we see that the choice d k = k, as in the work of Dupuis et al., 4 satisfies all requirements that we have imposed until now to get an asymptotically efficient change of measure based on Theorem 3.1, but it is by no means unique.
Next, we consider the boundary x 1 = 0 (and so x 2 . 0). Here, we want W d 1 (x) to be the minimum function, since all other functions are not designed for this part of the state space. Thus, for all x such that x 1 = 0, we want to have both The first inequality is equivalent to which unfortunately holds only for x 2 . (d 2 À d 1 )d= (a 2, 2 + g 2 ), provided that a 2, 2 . À g 2 . Since we do not want W d 2 (x) to be the minimum function for x such that This condition is equivalent to and if a 2, 2 ø 0 this inequality holds for all x 2 , so in particular for x 2 4 (d 2 À d 1 )d=(a 2, 2 + g 2 ). So we need As a result, using Equation (20) for k = 2 and ' = 3, we find for all x such that x 1 = 0 Since d 2 \ d 3 , the right-hand side tends to 0 as N ! ' and this implies that the weight factor r 2 (x) for all x such that x 1 = 0 tends to 0 as N ! '. Therefore, the change of measure that is designed for x 1 . 0, has hardly any influence when x 1 = 0.
The second inequality, Equation (30), is satisfied for all is the minimum function (which is not a limitation for the resulting change of measure to be asymptotically efficient since this function can be used for the whole state space). Recall that it is more important that W d 2 (x) is not the minimum function for all x 1 = 0 and x 2 . 0, and this requirement is satisfied.
To derive a lower bound on a 2, 1 and an upper bound on a 2, 2 , in contrast to the bounds in Equations (27) and (32), we use Condition 1 of Theorem 3.1. In this case, we only consider the term of the summation from this condition that considers k = 2, since this involves a 2 . That is, r 2 (x)H(x, a 2 ) should be non-negative for large enough N. Recall that we have derived an upper bound for r 2 (x) for all x such that x 1 = 0, see Equation (33), that tends to 0 as N ! '. Thus, for all x such that x 1 . 0 we see that, using Equation (18) for a = a 2 , we can only satisfy Condition 1, if we have le Àa 2, 1 + m 1 e a 2, 1 Àa 2, 2 + m 2 e a 2, 2 4 1, ð34Þ le Àa 2, 1 + m 1 e a 2, 1 Àa 2, 2 + m 2 41: It is clear that, as we have a 2, 2 ø 0 from Equation (32), the second inequality is implied by the first inequality. Using Lemma 3.2, Figure 5, and Equations (27) and (32), we find that the first inequality can only be satisfied when a 2, 1 2 ½Àg 1 , À g 2 and a 2, 2 2 ½0, g 1 À g 2 : Clearly, a 2, 1 and a 2, 2 have a dependence, e.g., a 2, 1 =À g 1 implies a 2, 2 = 0 and a 2, 1 =À g 2 implies a 2, 2 2 ½0, g 1 À g 2 , see Lemma 3.2 and Figure 5. The dependence can be found in Equation (34), however, this equation cannot be simplified. For future reference, we remark that as a result of these conditions on a 2 , we find from Equation (18), using a 2, 1 øÀ g 1 and a 2, 2 4 g 1 À g 2 , that since l + m 1 + m 2 = 1.

Summary and proof that all conditions are satisfied.
To summarize, we have found the following values for a 1 , a 2 , and a 3 that intuitively satisfy all conditions for an asymptotically efficient change of measure based on Theorem 3.1, see Table 2.
We show that these possibilities for a k , k = 1, 2, 3, indeed give an asymptotically efficient change of measure, by considering all conditions in Theorem 3.1. Recall that l + m 1 + m 2 = 1 implies H(x, a 3 ) = 0. To start with Condition 1, using Equations (23), (28), (33), and (36), we find the following lower bound ø À e maxfd 1 Àd 2 ;d 2 Àd 3 gd=e logð3m 1 Þ; where the last step follows since queue 2 is the bottleneck queue, and thus, m 2 4 m 1 . It follows that Condition 1 is satisfied. For Condition 2, we note that for all x 2 d e we have W d 1 (x) =À d 1 d, and thus, Condition 2 is satisfied To conclude with Condition 3, we find which goes to g 2 as N ! ' and hence all conditions are satisfied. Therefore, the change of measure for the given possible values of a 1 , a 2 , a 3 and d k , k = 1, 2, 3, is asymptotically efficient. Note that also the conditions in Remark 3.1 are satisfied, based on the construction.

Discussion.
It is clear that the choice of a 2 = ( À g 2 , 0) results in the change of measure from Dupuis et al., 4 using d k = k for k = 1, 2, 3. For this case, W d (x) is illustrated in Figure 2. We show some other examples of the piecewise affine function W d (x) below, where we let a 2 = ( À g 1 , 0) in Figure 7 and Table 2. Possibilities for a k when queue 2 is the bottleneck queue, provided that d 1 \ d 2 \ d 3 .
a k, 1 a k, 2 k Condition Figure 7. Display of W d (x) when queue 2 is the bottleneck queue, a 2 = ( À g 1 , 0) and d k = k, k = 1, 2, 3. a 2 = ( À g 2 , g 1 À g 2 ) in Figure 8. In both cases, we choose d k = k, k = 1, 2, 3. Choosing a 2, 2 = 0, we find a 2, 1 2 ½Àg 1 , À g 2 , and thus, W d (x) can be ''anything in between'' Figures 2 and 7. That is, the narrow area next to the x 1 -axis in Figure  2 can be as steep as in Figure 7, and anything in between, while the resulting change of measure still gives an asymptotically efficient estimator. Choosing a 2, 1 =À g 2 , we find a 2, 2 2 ½0, g 1 À g 2 , and thus, W d (x) can be ''anything in between' ' Figures 2 and 8. That is, the narrow area next to the x 1 -axis in Figure 2 can be slightly tilted, while the resulting change of measure still gives an asymptotically efficient estimator.
We remark that the change of measure that is used for x 2 . 0 is the state-independent change of measure from Parekh and Walrand. 1 The changes of measure that are found here, can be interpreted as ''protecting'' the x 1 -axis in the sense that we have to apply a different change of measure in that part of the state space. In Figure 7, we protect the x 1 -axis quite a lot, rather than in Figure 8, where we only protect it slightly. Recall that the most likely path goes along the x 2axis, where we have to apply the change of measure from Parekh and Walrand, 1 which also follows from the subsolution approach. It turns out that the change of measure along the most likely path is very important in the construction of an asymptotically efficient change of measure based on subsolutions and that along the most likely path there is no variation possible for the change of measure. However, it turns out that for all other parts of the state space it is possible to apply a (slightly) different change of measure than the one from Dupuis et al. 4 As we have seen in the work of Dupuis and Wang 6 and as we will see in Section 4.3, it is also possible to apply a different change of measure in the interior of the state space, rather than the same change of measure as along the most likely path or the change of measure that is applied along the x 1 -axis.

Three regions and queue 1 bottleneck
In this section, we again consider three functions W d k (x), k = 1, 2, 3, but in contrast to Section 4.1, we now let queue 1 be the bottleneck queue, so g 1 4 g 2 . The regions that we consider in this section are the same as in the previous section (see Table 1). By changing the bottleneck queue, the most likely path also changes and, therefore, we start the construction of an asymptotically efficient change of measure based on Theorem 3.1 by finding a 2 .
To find a 2 , we start with Condition 2, as in Section 4.1.1, and we use Remark 3.1 to determine a 2, 1 . When queue 1 is the bottleneck queue, the most likely path is covered by W d 2 (x). Therefore, we find a 2, 1 =À g 1 . Next, we use Condition 1 to determine a 2, 2 . Considering the term of the summation in this condition that involves k = 2, we need r 2 (x)H(x, a 2 ) to be nonnegative for large enough N. For those parts of the state space where x 1 = 0, we will find in the sequel that the weight factor r 2 (x) tends to 0 as N ! '. Since the function W d 2 (x) is designed for x 1 . 0, it is expected that in that case the weight factor r 2 (x) does not tend to 0 as N ! '. Thus, using Equation (18) for a = a 2 , for all x such that x 1 . 0 we see that we can only satisfy Condition 1 from Theorem 3.1 if we have m 1 + le Àa 2, 2 + m 2 e a 2, 2 4 1, where we used a 2, 1 =À g 1 . Using the second bullet of Lemma 3.2 and Lemma 3.3, we find that both of the above conditions hold when a 2, 2 = 0. Thus, to get an asymptotically efficient change of measure based on Theorem 3.1, we need a 2 = ( À g 1 , 0): For future reference, we remark that as a result of this condition we find, using Equation (18) and l + m 1 + m 2 = 1, 4.2.2. Finding a 1 . As in Section 4.1.2, we use the underlying idea of the construction of subsolutions to determine conditions on a 2 and d k for k = 1, 2, 3. Since this approach has been fully explained in Section 4.1.2, we will skip most of the details and highlight the results. Figure 8. Display of W d (x) when queue 2 is the bottleneck queue, a 2 = ( À g 2 , g 1 À g 2 ) and d k = k, k = 1, 2, 3.
By considering the origin of the state space, i.e., x 1 = x 2 = 0, we immediately find the same result as in Equation (24), since that result is independent of a 1 and a 2 .
Next, we consider the boundary x 2 = 0 (and so x 1 . 0) and recall that in that case, we want W d 2 (x) to be the minimum function. The conditions that follow from this observation can be found in Equations (25) and (26). Since we have determined that a 2 = ( À g 1 , 0), these conditions are equivalent to The first condition is satisfied when since then the right-hand side of the inequality is negative. Similarly to Equation (28) we find and so the weight factor r 1 (x) for x such that x 2 = 0 tends to zero when N ! '.
For the second condition, see Equation (39), we find that this is satisfied for all x 1 . (d 3 Àd 2 )d g 1 , and so for small values of x 1 the function W d 3 (x) is the minimum function. It turns out that this is not a problem for the resulting change of measure to be asymptotically efficient, since the function W d 3 (x) can be used throughout the whole state space and since W d for all x such that x 2 = 0. Therefore, it is ruled out that W d 1 (x) is the minimum function for these x, as desired.
At the boundary x 1 = 0 (and so x 2 . 0), we want W d 1 (x) to be the minimum function. As a result, we want Equations (29) and (30) to hold. These conditions are equivalent to Clearly, both conditions are satisfied when the second condition holds, since d 2 \ d 3 . In particular, we need a 1, 2 \ 0, since d 1 \ d 3 , so that both conditions are satisfied for all x 2 . (d 3 Àd 1 )d Àa 1, 2 . To prevent W d 2 (x) from being the minimum function for x such that x 1 = 0 and x 2 4 (d 3 Àd 1 )d Àa 1, 2 , we need W d 3 (x) \ W d 2 (x) for those states x. This condition is equivalent to Equation (31), and since a 2, 2 = 0 and d 2 \ d 3, this condition is always satisfied. As a result, we find for all x such that x 1 = 0 that Equation (33) holds.
To get a tighter condition for a 1, 2 we use Condition 2 from Theorem 3.1. That is, for all x on the exit boundary we need to have W e, d (x) 4 g(N ) such that g(N ) ! 0 when N ! '. As we have W e, d (x) ! W d (x) when N ! ', we in particular need W d 1 (x) to be non-positive for large enough N when x 1 = 0 and x 2 = 1, since W d 1 (x) is designed to be the minimum function at the boundary x 1 = 0 (and so x 2 . 0). Therefore we find Using the same condition of Theorem 3.1, we also derive an upper bound on a 1, 1 . Consider W d 2 (x) =À g 1 x 1 + g 1 À d 2 d, which is non-positive for all x 1 ø (g 1 À d 2 d)=g 1 . Thus, for all x 1 \ (g 1 À d 2 d)=g 1 , we need W d 1 (x) to be non-positive for all x 2 d e . On the exit boundary we have x 1 + x 2 = 1, so we find that for all x 2 d e , using Equation (42), which goes to zero as N ! ' if and only if a 1, 1 4À g 1 .
Combining with Equation (40) we find a 1, 1 =À g 1 : To conclude, we derive a lower bound on a 1, 2 using Condition 1. We only consider the term of the summation in Condition 1 that involves k = 1. That is, r 1 (x)H(x, a 1 ) should be non-negative for large enough N. We do not expect that r 1 (x) ! 0 when N ! ' for all x such that x 2 . 0, as W d 1 (x) is designed for x 2 . 0. Thus, using Equation (18) for a = a 1 , for all x such that x 2 . 0, we see that we can only satisfy Condition 1, if we have m 1 + le Àa 1, 2 + m 2 e a 1, 2 4 1, ð43Þ where we have used that a 1, 1 =À g 1 . Using Equation (42), it follows that the second inequality is implied by the first inequality. Using Lemma 3.2, the second bullet, we find that Equation (43) is satisfied when a 1, 2 2 ½Àg 2 , À g 1 , and so is Equation (44). For future reference, we remark that as a result of this condition we find, using Equation (18) and a 1, 2 øÀ g 2 , since l + m 1 + m 2 = 1.

4.2.3.
Summary and proof that all conditions are satisfied. Summarizing, we have found the following values for a 1 , a 2 , and a 3 that intuitively satisfy all conditions for an asymptotically efficient change of measure based on Theorem 3.1, see Table 3. Again, we show by considering all conditions in that theorem, indeed these possibilities for a k , k = 1, 2, 3, give an asymptotically efficient change of measure. We start with Condition 1. First, recall that H(x, a 3 ) = 0, since a 3 = 0 and l + m 1 + m 2 = 1. Then, from Equations (33),  (37), with g 2 replaced by g 1 , since queue 1 is the bottleneck queue. Therefore, the change of measure for the given possible values of a 1 , a 2 , a 3 , and d k , k = 1, 2, 3 is asymptotically efficient.

Discussion.
For d k = k for k = 1, 2, 3, it is clear that choosing a 1 = ( À g 1 , À g 2 ) results in the change of measure from Dupuis and Wang 6 and, for this case, W d (x) is illustrated in Figure 4. Another example for W d (x), where we choose a 1 = ( À g 1 , À g 1 ), is illustrated in Figure 9.
As can be seen in Table 3, a 1, 2 can range from Àg 2 to Àg 1 , and thus, W d (x) can be ''anything in between'' Figures 4 and 9. That is, the area next to the x 2 -axis and in the interior in Figure 9 can be as steep as in Figure 4 and anything in between. Note that Figure 9 is very similar to Figure 2 in the sense that g 2 is replaced by g 1 , and so in particular when g 1 = g 2 the change of measure based on this function W d (x) can be used.
When queue 1 is the bottleneck queue, we find that the state-independent change of measure from Parekh and Walrand 1 is only used around x 2 = 0. While when queue 2 is the bottleneck queue, the change of measure based on W e, d (x) can be interpreted as that we have to ''protect'' the x 1 -axis, see also Section 4.1, it seems natural that we have to ''protect'' the x 2 -axis when queue 1 is the bottleneck queue. Looking at Figures 4 and 9, we see that this is clearly not the case. One could argue that in this way we are ''overprotecting'' the x 2 -axis. However, it turns out that for an asymptotically efficient change of measure based on subsolutions we only need to have the state-independent change of measure from Parekh and Walrand 1 along the most likely path, which in this case is along the x 1 -axis.

Four regions and queue 2 bottleneck
In this section, the starting point is to use four functions W d k (x), k = 1, 2, 3, 4, and we let queue 2 be the bottleneck queue, thus g 2 4 g 1 . The regions that we consider can be found in Table 4, and the only difference compared with Section 4.1, where we considered three regions, is that we now have a dedicated region in case x j = 0 for some queue j.
The zero change of measure is used when both x 1 = x 2 = 0, for similar reasons as in Section 4.1. The only difference is that we now set a 4 = 0, instead of a 3 = 0 as in Section 4.1.
The derivation for using four regions is very similar as when using two regions and, therefore, it is omitted in this paper. A detailed derivation can be found in Buijsrogge. 12 Table 3. Possibilities for a k when queue 1 is the bottleneck queue, provided that d 1 \ d 2 \ d 3 .
4.3.1. Summary and proof that all conditions are satisfied. We have found the following values for a 1 , a 2 , a 3 , and a 4 that intuitively satisfy all requirements for an asymptotically efficient change of measure based on Theorem 3.1, see Table 5. We show that these possibilities for a k , k = 1, 2, 3, 4, indeed give an asymptotically efficient change of measure, by considering all conditions in Theorem 3.1. Recall that l + m 1 + m 2 = 1 implies H(x, a 4 ) = 0. To start with Condition 1 of Theorem 3.1, we find that where the first step follows by conveniently substituting the upper and lower bounds on a k from Table 5 (37). Therefore, the change of measure for the given possible values of a 1 , a 2 , a 3 , a 4 , and d k , k = 1, 2, 3, 4, is asymptotically efficient. Table 5, we again see that the change of measure from Parekh and Walrand 1 is used for x 1 = 0 and x 2 . 0, as expected, since it is the region in which the most likely path lies. Clearly, the result in the work of Dupuis and Wang 6 satisfies the conditions in Table  5. Recall that in that paper also different overflow probabilities are considered, which imposes additional constraints for an asymptotically efficient change of measure based on subsolutions. The change of measure from Dupuis and Wang, 6 where in particular a 1 = ( À g 1 , À g 2 ) and a 3 = ( À g 1 , 0), can be found in Figure 3.

Discussion. From
From the results in Table 5, it is also clear that we could adapt Figure 3 such that along the x 1 -axis we have a similar affine function as in Figure 7 or Figure 8. However, when choosing a 3 = ( À g 2 , g 1 À g 2 ), as we did for a 2 in Figure 8, we need a 1, 1 =À g 2 and so the functions W d 1 (x) and W d 2 (x) are almost the same. That is, only their constants d 1 and d 2 would differ.
Summarizing, at the x 1 -axis, we can make similar adaptations as when we considered three regions in Section 4.1. The difference is that now also in the interior we can ''push'' the function a bit more toward the origin when we compare with Figure 2

Four regions and queue 1 bottleneck
In this section, we again consider four functions W d k (x), k = 1, 2, 3, 4, but in contrast to Section 4.3, we let queue 1 be the bottleneck queue, i.e., g 1 4 g 2 . The regions that we consider in this section are the same as in the previous section (see Table 4). Again, the derivation for using four regions is very similar as when using two regions and therefore it is omitted in this paper. Recall that a detailed derivation can be found in the work of Buijsrogge. 12 4.4.1. Summary and proof that all conditions are satisfied. We have found the following values for a 1 , a 2 , a 3 , and a 4 that intuitively satisfy all conditions for an asymptotically efficient change of measure based on Theorem 3.1, see Table 6.
Again, we show by considering all conditions in that theorem, indeed these possibilities for a k , k = 1, 2, 3, 4, give an asymptotically efficient change of measure. We start with Condition 1, and recall that H(x, a 4 ) = 0. Then, we find that Table 4. Overview of proposed regions for the case r = 4.

Discussion. From
Comparing with Figure 9, we see that in the area close to the x 2 -axis, we can apply a slightly different change of measure. Of course, there are several other possibilities.

Conclusion
In this paper, we determined sufficient conditions for subsolution-based changes of measure to give asymptotically efficient estimators. As a result, for the two-node MjMj1 tandem queue, we explicitly gave a whole family (continuum) of changes of measure that all lead to asymptotically efficient estimators, and the previously known changes of measure are just three members of this family. For d-node tandem queues, it seems likely that we can use a similar analysis to find a family of changes of measure that are asymptotically efficient.
For the case d = 2, we like to highlight one particular change of measure based on the subsolution W e, d (x) in Equation (6) (via either Equation (3) or Equation (8)), that uses the following three functions: W d 3 (x) = g À 3d, Table 6. Possibilities for a k when queue 1 is the bottleneck queue, provided that d 1 \ d 2 \ d 3 \ d 4 . There is no simple expression for the maximum value that a 2, 1 could possibly attain and it is therefore denoted with *. It should satisfy the conditions stated (see also Figure 6).
where we recall that g = minfg 1 , g 2 g. We note that this is the same subsolution as in the work of Dupuis et al. 4 if queue 2 is the bottleneck queue, while it also works when queue 1 is the bottleneck queue (with g 2 replaced by g 1 ). This matches nicely with the known fact that interchanging the queues leaves our probability of interest unchanged. From an implementation point of view with respect to subsolutions, in general, it makes sense to use as few regions as possible, i.e., three regions (or d + 1 in the dnode case). However, when the event of interest is not total buffer overflow, but, e.g., individual buffer overflow or simultaneous buffer overflows, it may be more useful to implement the change of measure from Dupuis and Wang 6 that is based on four regions.
Finally, we mention that future work could aim at investigating whether or not the method generalizes to more general models, which we expect to be the case. In the work of Buijsrogge et al., 7 something similar has already been done for non-Markovian tandem queues, but one could also think about more general (non-Markovian) networks, for which we expect similar results to hold.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by the Netherlands Organisation for Scientific Research (NWO), project number 613.001.105.