Abstract
We consider a model for a single link in a circuit-switched network. The link has circuits, and the input consists of offered calls of two types, that we call primary and secondary traffic. Of the links, are reserved for primary traffic. We assume that both traffic types arrive as Poisson arrival streams. Assuming that is large and , the arrival rate of primary traffic is , while that of secondary traffic is smaller, of the order . The holding times of the primary calls are assumed to be exponentially distributed with unit mean. Those of the secondary calls are exponentially distributed with a large mean, that is, . Thus, the primary calls have fast arrivals and fast service, compared to the secondary calls. The loads for both traffic types are comparable , and we assume that the system is โcritically loadedโ; that is, the system's capacity is approximately equal to the total load. We analyze asymptotically the steady state probability that (resp., ) circuits are occupied by primary (resp., secondary) calls. In particular, we obtain two-term asymptotic approximations to the blocking probabilities for both traffic types.
1. Introduction
A classic model in teletraffic is the Erlang loss model. Here, we have servers (or circuits), and customers (telephone calls) arrive as a Poisson process with rate parameter . The arriving customer takes one of the circuits if one is available, and if they are all occupied then the call is blocked and lost. When occupying a circuit, the customer has an exponentially distributed holding time whose mean we take as the unit of time. It is well known that the steady state probability that circuits are occupied is the truncated Poisson distribution, that is, , , with . This model dates back to circa 1918 [1]. When , we obtain the steady state blocking probability. The transient probability distribution is much more complicated, but it can be computed in terms of special functions (see [2]).
Over the years, many generalizations of the basic model have been analyzed, including networks of such loss models (see [3, 4]). One important extension is that of trunk reservation, which is fundamental in the analysis of circuit-switched communication networks. Here, we consider a model with circuits that are used by the two types of customers (or offered calls). We refer to these as primary (or high-priority) calls and secondary (or low-priority) calls. They arrive as Poisson arrival streams with respective rates and . Of the circuits, are reserved for primary calls. Thus, if a high-priority call arrives, it is blocked if all circuits are busy, while a low-priority call is blocked if at least circuits are busy. All calls are assumed to have independent and exponentially distributed holding times, with respective means and . The total load on the system is . If this exceeds , then typically all the circuits are busy (an overloaded link), while if , typically some circuits are free (an underloaded link). An interesting situation is when is large and ; we refer to this as โcritical loading.โ
Previous work on this and related models includes Mitra and Gibbens [5] who considered the asymptotic regime with , , and (thus, secondary calls are less frequent than primary calls). We thus have ; so this is an example of critical loading. They analyzed a single link and used their results to obtain approximations for more complicated loss networks with a distributed, state-dependent, dynamic routing strategy. Related work appears in [6, 7], and optimization and control policies for such problems were analyzed by Hunt and Laws [8].
Of fundamental importance in this model is the probability (resp., ) that a primary (resp., secondary) call is blocked and lost in the steady state. Roberts [9, 10] obtained approximations to these blocking probabilities, which are based on a certain recursion which is exact for special cases of the model parameters, but not for all cases. Morrison [11] investigated this model for and , and obtained the blocking probabilities as asymptotic series in powers of . This led to a better understanding of the asymptotic validity of Roberts' approximation(s). However, the coefficients in the asymptotic series in [11] were not explicit, as their calculation still involves recursively solving an infinite system of differential equations. But, if it is further assumed that is small, the blocking probabilities were obtained more explicitly in terms of parabolic cylinder functions. Also, if rather than , explicit results are obtained without the small assumption.
In [12], we analyzed the case , with and with the arrival rates and both . Expressions for the blocking probabilities were obtained for the overloaded and underloaded cases. In the first case, both blocking probabilities remain as , while in the second case they are exponentially small. In this paper, we investigate the case of critical loading, where again but now with . We will also assume that (secondary calls are less frequent than primary ones) but now with ; that is, secondary calls have large holding times. Thus, primary calls have faster arrivals and faster service. Note that the loads due to the primary and secondary calls remain asymptotically comparable with this scaling. In this asymptotic regime, we are able to obtain explicit analytic expressions for the first two terms in the expansions in powers of for the blocking probabilities, which involve readily evaluated definite integrals.
We are currently investigating situations where the secondary calls are the ones with fast arrivals and service [13]. Here, the asymptotic structure of the problem turns out to be quite different.
We comment that the basic problem to be solved is a two-dimensional difference equation (cf. (2.1)), with discontinuities of the coefficient functions at various boundaries and an interface. Such a problem appears to be very difficult or impossible to solve exactly. From a numerical point of view, the problem corresponds to solving roughly linear equations. A good general method such as Gaussian elimination has computational complexity or . Some methods that use the sparseness of the system and some iteration procedures may improve this to or . The purpose of our asymptotic analysis is to obtain reasonable approximations whose numerical evaluation has computational complexity that is independent of , and also to obtain explicit formulas that show the dependence of the stationary distribution and blocking probabilities on the model parameters (i.e., , , and the arrival and service rates).
The paper is organized as follows. In Section 2, the problem is stated more precisely and the basic equations are obtained. Here, we summarize our main results, which are derived in detail in Sections 3 and 5. In Section 3, we obtain the leading terms for the blocking probabilities. In Section 4, we relate the present results to the ones in [11, 12] using asymptotic matching. The first-order correction terms to the blocking probabilities are derived in Section 5, while in Section 6 we present some numerical studies to assess the accuracy of the asymptotics.
2. Statement of the Problem and Summary of Results
We denote by the number of servers serving high-priority customers, and by the number of servers serving low-priority ones. The total number of servers (circuits) is of which are reserved for the high-priority customers. Thus, if , a newly arriving high-priority customer (call) is lost; if , then a newly arriving low-priority call is lost. The high- and low-priority customers arrive as independent Poisson processes, with respective rates and . The service times are exponentially distributed with respective means and . Thus, the unit of time is taken as the service rate of the high-priority customers.
We denote the steady state joint distribution of the numbers of servers used by the two priority classes by . We let be the indicator function on the event . Then, from the description of the model, we obtain the following balance equation:
This applies over the domain Thus, we may view the problem as solving a second-order difference equation in two variables, over the triangle and the oblique strip , with the two subdomains separated by the โinterfaceโ . There are also boundary conditions inherent in (2.1), along , , , and . The normalization condition is
Of particular interest are the blocking probabilities for the high-priority customers, defined by and for the low-priority ones, defined by Note that we clearly have .
We analyze the problem in the asymptotic limit where
We furthermore assume that the arrival rate of high-priority customers is large, of the same magnitude as , and then scale the other rate parameters as
Thus, the arrival rate of low-priority customers is large, but only of the order . The service times of these customers however are also large, and the total load due to low-priority customers is . Also, we have so that the total load due to all customers is roughly equal to the capacity of the system. Hence, this asymptotic limit may certainly be considered as โheavy trafficโ or โcritical loading.โ
Once we input the model parameters , , , , and , we can compute , , and from (2.7). For some of our numerical studies, it is desirable to fix , , , , and and then vary the original rate parameters, which are computed by inverting (2.7) using
and and are obtained from (2.7) once is known.
We next scale and as
with and and are taken as . We consider the scaled state space with . Then, since is large and is , the domain in the plane is the triangle ; . The scaling (2.9) corresponds to a small neighborhood of the point
However, this is where most of the probability mass accumulates in this asymptotic limit, and the analysis of this range is sufficient to obtain the blocking probabilities . We will obtain these as asymptotic series in powers of .
Using (2.9) and (2.10) in (2.1), we obtain
Note that in this asymptotic scaling, the indicator functions , , and may be replaced by one, since these correspond to boundaries that are far from the point in (2.11). However, the interface and the boundary are evident in (2.12) and will play a large part in the analysis. In Section 3, we will analyze (2.12), and then also consider a second scale where is large and negative. This second scale will lead to a diffusion equation in two variables.
The main results are as follows. For and , we obtain the approximation , with
where The correction term in (2.13) takes the form for and the form for , where For the latter function, we have , where is given by (5.5) in terms of and . Then, is given by (5.66)โ(5.68) (with and defined in (5.6) and (5.56)), and the constant in (5.66) can be obtained by using (5.78) in (5.83).
On the scale, where and , we find that
where and is as above.
The analysis of these two scales leads to two-term approximations to the blocking probabilities in (2.4) and (2.5). To leading order, these are
where the integral is evaluated using (2.14) and (2.15) (see also (3.49)). The correction terms follow from (5.80) and (5.81), using also (5.74).
We note that the numerical evaluation of the leading order asymptotic results involves only the integrals in (2.15) and (2.17). The correction terms involve numerically evaluating some double or triple integrals (cf. (5.78)), but the computational complexity of evaluating the asymptotic results is independent of (or ).
3. Asymptotic Analysis: Leading Terms
We consider (2.12) and assume that for the probabilities have the expansion
In this section, we focus on the leading term, but its calculation will necessitate that we also analyze the problem for . This correction term is calculated completely in Section 5. We will also need to couple the analysis of the scale to that where is large and negative, with .
Using (3.1) in (2.12) and equating coefficients of powers of , we obtain to leading order
and this applies for all . This is a simple difference equation with a boundary condition at . Its most general solution is
where is to be determined.
For , we obtain from (3.1), (3.3), and (2.12) the problem
whose solution is Here, is not yet determined. For , the terms in (2.12) yield
which simplifies to and hence
Here, we imposed continuity between (3.5) and (3.7) along , and is another function not yet determined. By setting in (2.12) and comparing terms of order , we obtain
Using (3.7) to compute and (3.5) for and , we obtain
so that (3.7) becomes We next consider the problem (2.12) for , with the scaling
Note that this still corresponds to a local approximation near the point in (2.11). In terms of , we let
and (2.12), upon multiplying by , becomes
We note that corresponds to , and that for the indicator functions in (2.12) can all be replaced by one.
We assume that has an expansion in the form
Using (3.14) in (3.13), we obtain to leading order the PDE This is a parabolic PDE whose solution is facilitated by the change of variables Using (3.16), (3.15) becomes The most general solution, that decays exponentially as , is
and hence We will determine shortly.
We observe that on the scale, with , expansion (3.14) becomes
Comparing this to (3.1), for we conclude from (3.3) that
and, from (3.10) that It follows that which is a boundary condition for the PDE (3.15) along . In terms of , this becomes , but this holds automatically (for any ) in view of (3.18). To determine , we must analyze the correction term in (3.14).
From (3.13) and (3.14), we find that the first correction term satisfies the PDE
Switching to the variables and then using (3.18), we get
Integrating (3.25) with respect to yields
Setting in (3.26), we obtain
We will show that (3.27) leads to a differential equation for . But, we must first consider . In view of (3.20), this term arises as a part of the term in the expansion on the scale. We therefore return to (3.1) and (2.12).
For , we obtain from (2.12) and (3.1), at order , the equation
Using (3.3) and (3.5), we find that (3.28) has a solution in the form where It follows that We note that and
For , the terms in (2.12), with (3.1), yield
Using (3.3) and (3.10), we find after some calculation that (3.34) simplifies to and hence Setting in (2.12), we obtain at order the interface relation
Using (3.36), with and , and (3.33), we obtain from (3.37) after some calculation
Now, comparing the terms in (3.20) and (3.1) with (3.36), we conclude that
But then (cf. (3.27)) involves only, and since and are related via (3.21), we obtain from (3.27) the following ODE for :
This equation may be written as a perfect derivative as Integrating once and requiring to vanish as yield
and hence where is a constant, which will be fixed by normalization.
We use (2.10) and (3.12) in the normalization sum (2.3), and then use the Euler-MacLaurin formula to approximate sums by integrals. To leading order, this yields
We use (3.19) and (3.43) and evaluate one of the two integrals in (3.44) using integration by parts, with the result To summarize the calculation, we have shown that on the scale in (3.11) we have
with given by (3.45). Note that is , but the probabilities are spread out over an range near the point given by (2.11). On the scale in (2.9), we have obtained
which applies for and is independent of . To calculate correction terms to (3.46) and (3.47), we would need to find in (3.5) and (3.10), and solve (3.24) for . This ultimately involves calculating the and terms in (3.1) and the term in (3.14); this is done in Section 5.
We next calculate the blocking probabilities in (2.4) and (2.5). Evaluating these sums requires the expansion on the scale. Again, approximating sums by integrals and using the scaling (2.9), we obtain
as . From (3.21), (3.43), and (3.45), we then obtain
The numerical accuracy of (3.48) is investigated in Section 6. This completes the analysis of the leading terms.
4. Consistency with Previous Results
In [11], Morrison studied the current model with the scaling and being either or . For the latter case, we define from (see [11, equations (7.16)โ(7.18)]) and then
where We show that (4.1) matches asymptotically to (3.48), in the limit where . When , we can simplify the integrals in both the numerator and denominator in (3.49), as the factor has the effect of freezing the remaining integrands at . Thus, by the Laplace method, we obtain
But, in view of (2.7), so that . Since , (4.1) agrees with (3.48)-(3.49) and (4.3).
In [12], Knessl and Morrison analyzed the model in the limit with and . The total load is thus and the cases (resp., ) correspond to an underloaded (resp., overloaded) link. The case of critical loading was not considered in [12].
We consider first the asymptotic matching of (3.48)-(3.49) to the underloaded case in [12]. We note that the parameters and are related to the current ones by
For the matching, we must thus let with , and . The results in [12] for and were as follows:
where In view of (4.4), we have, in the matching region, , , , and
Hence, and (4.5) becomes We show that (4.9) agrees with (3.48)-(3.49) when the latter is expanded for . In (3.49), the main contribution to the integral in the denominator comes from , and in the numerator from , where . Thus, we can approximate by everywhere, and obtain
With (4.10), (3.48) agrees with (4.9).
Next, we consider the overloaded case in [12], for which we obtained
where and is the solution of
We note that by using (4.12) in (4.13), is a particular root of a polynomial. For the matching, we will take in (4.11) and in (3.48)-(3.49). As , we will have and this will allow us to simplify (4.11). Let us set
with which , and thus
Furthermore, Using (4.14)โ(4.16) in (4.13), we find after some calculation that and then (4.11), for , simplify to
We expand (3.49) for . Scaling , we obtain
For , we use the asymptotic approximation
and conclude that both integrals in (4.19) have their major contribution from where But then by the Laplace method, we have and (3.48) agrees with (4.18).
This completes the matching verifications.
5. Correction Terms
We will compute the terms in expansions (3.1) and (3.14), and then obtain corrections to the blocking probabilities.
We first consider (3.13). Using the relations and , we rewrite (3.13) in terms of . Defining
we obtain from (3.13) so that satisfies
where are the operators
Before analyzing (5.3), we obtain a more complete description of . We recall that is known completely, in view of (3.18), (3.43), and (3.45). But, we computed only partially as the combination in (3.26). We integrate (3.26) to get
where We observe that satisfies
and decays as . We also have
The function will be determined shortly (actually, not very shortly, but only after a lengthy calculation).
We evaluate the right side of (5.3) more explicitly. Some terms are expressible as derivatives with respect to , while the ones involving derivatives in may be evaluated using (3.18) and (5.5). Then, (5.3) becomes
From (3.18), we have and by direct calculation With the above, we integrate (5.9) to get where We will show that the calculation of will require only that we evaluate (5.12) along . However, we must first reconsider the scale and analyze at least partly the term in (3.1) (i.e., the coefficient of in the series).
Returning to (2.12) with the expansion (3.1), we find that for , satisfies
We recall that for , is given by (3.3), by (3.5), and by (3.29). Let us write
We then rewrite (5.14) as which holds for . We sum (5.16) for and use the identities
After some rearrangement, we obtain where all derivatives are with respect to .
Setting in (2.12) and using again expansion (3.1), we obtain at the relation
Next, we consider (2.12) for , and recall that is given by (3.3) for all , is given by (3.7) or (3.10), and
where is in (3.27), and (3.36) yields
We subtract from both sides of (5.19) and substitute (5.18) into the resulting equation. We then use (3.3), (3.10), (5.20), (5.15), and (3.29), and after some calculation, obtain
From (2.12), for , we obtain the following problem for :
Here, we used in this range of , and all functions in (5.23) are evaluated at . After rearranging (5.23) and using (5.20), (3.3), and (3.7) to evaluate for , we are led to
We sum (5.24) from to (with ) to obtain
Note that (5.25) remains true if . Using (5.22) in (5.25) gives us
Using our previous results for , , and , we have
Adding (5.26) and (5.27), we obtain an explicit expression for (in terms of , , , , , and ) which is quadratic in . By summing from to , we get
This holds for all . We write (5.28) as
where , and , , and may be identified from (5.28).
We recall that for , and the coefficient of in the expansion (3.20) is
Comparing (5.29) to (5.30) yields In view of (3.19), we have while (5.28) shows that But, from (5.21) and (3.21), we get , and from (3.9), . Then, we can easily verify that (5.32) is consistent with (5.33). Also, from (5.28), we find that
We show that this is the same as . We recall that is given by (5.5) and is expressed in terms of in (3.38). From (3.38) and (3.21), it follows that
where is the operator in (5.13), with replaced by . The second equality in (5.35) follows from (3.40) and (5.6). Using (3.9), we obtain
which when combined with (5.35) gives We use (5.37) in (5.34), also noting that which follows from (5.5). Then, (5.34) becomes Here, we also used . Now, from (5.21) and (3.9), we obtain
Using (5.38) to compute and (5.40), we get
In view of (5.5), the above is the same as .
We next examine the relation . We will use this to ultimately obtain in (5.5), which will complete the determination of the correction terms in (3.1) and (3.14). We first note from (5.28) and (5.29) that
We solve (5.42) for the combination , that we rewrite as
where Next we note that and may be integrated explicitly and we write
where ? It follows that The right side of (5.48) was computed in (5.12).
Using , (5.38), and the identities , , and , we evaluate in terms of , and then use (5.12). After some simplification, this leads to
Since this must be equal to , we try to write the right side of (5.49) as a perfect derivative. To this end, we note that
Adding to (5.49), we rewrite that equation as
We next evaluate in terms of , and in terms of and , and then we integrate (5.51) (and thus explicitly obtain ). Since and is in (5.21), we have
We thus note that is canceled by the bracketed term that precedes it in (5.51).
Using (5.46), we explicitly calculate , recalling that is given by (5.15), and can be computed from (3.30) and (3.32). After some cancellation of terms, we obtain
Using (5.52) and (5.53), we integrate (5.51), subject to the condition that the solution decays exponentially as . Hence,
What remains is a linear first-order ordinary differential equation for , which is readily solved by multiplying by the integrating factor:
We introduce the notation and note that Then, we have From (3.43), we have and thus With (5.60), we have Furthermore, To obtain (5.64), we used and integrated by parts. Combining (5.58) and (5.61)โ(5.64), we integrate (5.54) to get
Here, is a constant that will be fixed by normalization.
We thus write as
where where is as in (5.13).
We next determine by normalization and then obtain correction terms to the blocking probabilities and . This requires that we evaluate the integrals and . From (5.5), we have
and from (3.43), we calculate and obtain
Consider the contribution to that comes from the term in (5.69) and the part of that is proportional to (cf. (5.66)). We obtain
Here, we wrote and integrated by parts. We also have Using (5.71) and (5.73), we integrate (5.69) and get Here, we also used (5.70) to eliminate from the expression.
Next, we consider
Integration by parts shows that From (5.68) and (5.70), we obtain
With (5.66), (5.68), and (5.75)โ(5.77), we obtain Using and (5.6), we can show that
which helps in the numerical evaluation of (5.78).
Finally, we calculate the blocking probabilities. Using (3.1), (3.3), and (3.5), we obtain where the last integral is given by (5.74) in terms of and . To obtain (5.80), we used the scaling (2.9) and (2.10) in (2.4), and approximated the sum by an integral. Since the integrand has exponentially small tails, the finite limits in (2.4) may be replaced by infinite ones, with an error, that is, for all .
Similarly, we obtain as
Finally, we determine from the normalization (2.3). Again, using (2.9), (2.10), (3.1), and (3.12), we obtain
Here, the second integral in (5.82) comes from the Euler-MacLaurin approximation as we go from a discrete sum to an integral over . Note that the expansion on the scale, for , leads to the third term in (5.82). For , the expansion on the scale contains that on the scale. The leading term in (5.82) regains (3.45) and determines . The terms lead to
In view of (5.78), (5.83) may be viewed as a linear equation for , and thus all the correction terms are now known fully.
To summarize the calculations in this section, we have determined in (5.66)โ(5.68), with computed from (5.78) and (5.83). In terms of , the second term in the expansion on the (or ) scale is given by (5.5). Then, and the second terms on the scale are given by (3.5) for , and by (3.10) for (with in (3.21) and (3.43)).
6. Numerical Studies
We test the numerical accuracy of our asymptotic expansions, focusing on the blocking probabilities and . The numerical results are obtained by solving the linear system (2.1) with the normalization (2.3). We simply omitted the equation with in (2.1) and replaced it by (2.3), thus obtaining an inhomogeneous problem with a unique solution.
We solved (2.1) by two different methods. First, we simply used the program MAPLE to solve the linear system numerically. We also tried an iteration method of the form where is the basic equation (2.1). Starting from some initial guess and iterating up to correspond to solving approximately for the transient solution for this model, from time to . We verified that choosing sufficiently large leads to the same results as the MAPLE solution of (2.1).
Since the asymptotic results are expansions in powers of with coefficients expressed in terms of , we input the five parameters , calculate from (2.8), then and from (2.7), and solve (2.1) numerically. In Table 1, we have , , , and , and we compare the exact (numerical) results for and with the one- and two-term asymptotic approximations. We give various values of and also tabulate the corresponding values of , as computed from (2.8). We see that the one-term approximations always overestimate the true values, while the two-term approximations underestimate them. The two-term approximations are more accurate especially for the second blocking probability and for larger values of . In Table 2, we have , and in Table 3, , with the other parameter values unchanged. With decreasing (which corresponds to increasing the total load (cf. (2.7)), we get similar results, but the overall asymptotics (both one- and two-term) are getting somewhat worse. Also note that the two-term approximations may sometimes lead to negative answers, and this is explained in what follows.
We next consider a different purely numerical approach to estimating the coefficients in the expansions of the . We choose some , and for , , and , we equate
Note that and , for fixed values of . Thus, (6.1) may be viewed as systems of linear equations for the and , respectively. This allows us to numerically estimate the first three coefficients in the asymptotic series. In Table 4, we consider in the range of 5 to 70, and give the and , fixing . We see that the sequence of and does appear to converge as ; the convergence of and is slower, and that of and is even slower. The asymptotic results in Sections 3 and 5 show that for these parameter values
This is in good agreement with Table 4. The data in Table 4 also give a rough estimate of the third terms in the expansions of the blocking probabilities. These can be computed analytically by continuing our expansions further, but the calculations are too foreboding.
In Tables 5 and 6, we again give the and for between 5 and 70, but now with and , respectively. For these values, our asymptotic analysis predicts that
Again this is in good agreement with the apparent limiting values of , , , and as . The data in Tables 4โ6 show that the expansions do indeed appear to be in powers of , and that we correctly computed the leading two terms. Note that in each case, the second coefficient ( and ) is negative, while the first and third ones are positive. This is consistent with the fact that in Tables 1โ3 the leading terms always overestimate the exact answer, while the two-term approximations underestimate it. As we decrease , the ratio increase, as do (though these are always larger). Hence, we expect that decreasing leads to further cancellation between the first and second terms in the asymptotic series, and this again is in agreement with the data in Tables 1โ3. It also explains why the two-term asymptotic approximations to sometimes lead to negative answers, for moderate values.
To summarize, we have shown that the asymptotic approximations are reasonably accurate, though certainly not excellent, and that there is merit to computing the correction terms unless is quite small. For small , however, the one-term approximations may be superior, as the two-term approximations may lead to negative answers. The accuracy of the asymptotic approximations presumably increases as increases further, and the two-term approximations are presumably better than the one-term approximations. However, limitations of the available computing facilities have so far prevented the evaluation of the exact numerical results for larger values of .
Acknowledgments
J. A. Morrison is a Consultant at Alcatel-Lucent Bell Laboratories. C. Knesslโs work was partially supported by NSF Grants DMS 02-02815 and DMS 05-03745.