ISRN Applied MathematicsVolume 2012 (2012), Article ID 340306, 22 pagesdoi:10.5402/2012/340306

Research Article

Application of an Queuing Model to WiMAX-Based Wireless Mesh Networks: Grid Topology Case Study

Vahid Sattari Naeini and Naser Movahhedinia

Department of Computer Engineering, University of Isfahan, Hezar Jerib Avenue, Isfahan 8174673441, Iran

Received 25 November 2011; Accepted 4 January 2012

Academic Editors: Z. Duan and A. El-Sayed

Copyright © 2012 Vahid Sattari Naeini and Naser Movahhedinia. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Wireless mesh networking is an effective approach to reach high performance in the last mile of broadband Internet access. The mesh structure is the basic step toward providing cost-effective, dynamic, and high-bandwidth wireless connection. In this paper, WiMAX-like wireless mesh network is considered, emphasizing the grid arrangement which is the general topology described in the literature. To evaluate the performance of the conventional and proposed scheduling algorithms, each link is modeled using an M/D/1 queue and a virtual node concept is introduced to describe comparable performance metrics for the system. Performance measures of the system in addition to the simulation results are assessed in terms of the network length and the arrival rates.

1. Introduction

Wireless mesh network (WMN) is a cost-effective multihop network, usually comprising some subscriber stations (SS) and one base station (BS). The multihop feature of WMN can increase its range of accessibility, in addition to reliability due to dynamic self-organizing and self-configuring features. Each SS can communicate with other SSs in its transmission range and the BS that connects the mesh network to the backbone. While some SSs only relay traffic flows of other SSs, some others carry their underlying subnetwork traffic flows as well.

IEEE 802.16d has a potential to replace its tenacious rivals [1]. Its optional mesh mode transmits and routes various traffic streams among SSs directly in a distributed manner. Moreover, in its centralized mode, traffic stream routes to the network backbone through the BS. Since most of the traffic streams are between SSs and wired backbone network, the centralized mode is the most important [2–4] and we evaluate the system in this mode.

In this paper, we consider an IEEE 802.16-like WMN and analyze the performance of a grid topology, which is generally used in WMNs. We introduce a queuing model to evaluate the average delay and also investigate some other issues influencing the delay constraint such as scheduling algorithm and position of the BS.

The paper is organized as follows. In Section 2, a brief survey of the previous research in this area is presented. In Section 3, the topologies used for introducing the model are described. Section 4 elaborates the IEEE 802.16d mesh mode with respect to its central scheduling procedure. The analytical model is presented in the Section 5. In Section 6, the solution is proposed and verified with two simple examples. Impact of the scheduling algorithm is also considered in this section. Performance evaluation of the proposed model and simulation results are discussed in Sections 7 and 8, respectively. Finally, the paper is concluded in Section 9.

2. Related Works

A comprehensive survey of WMNs has been presented in [5]. Performance analysis of WMN has been considered in a few recent papers. Reference [4] analyzes a WMN by combining M/M/C/∞, M/M/N(s)/K(s), and M/M/1/∞ queuing models. The network throughput, average packet delay, and packet loss are evaluated using the model. Performance modeling and analysis of integrated WLANs and mesh networks are considered in [2]. The analytical model utilizes the Poisson arrival and service rates at each mobile station. The average queuing delay in multichannel multi-interface WMNs is minimized in [6] based on an M/M/C/∞ queue.

None of these models are suitable for IEEE 802.16 mesh network in which the frame duration is fixed. So, the main contribution of the paper is finding a suitable queuing model that can be considered in WiMAX-based WMN.

Time division multiple access (TDMA) scheduling in central mesh mode of IEEE 802.16 has been considered in [7]. The authors of [8] have represented a survey on scheduling in WiMAX mesh mode. Since the basic scheduling algorithm of IEEE 802.16 has two drawbacks (stated below), some works have been carried out, challenging the basic algorithm. In [9], the Bellman-Ford scheduling algorithm shows worthy results on the return path delay, which is also bandwidth optimal. An interference-aware research framework introduced in [10] uses a blocking metric to improve spectral utilization. Scheduling problem in WiMAX mesh network is formulated and solved in [11] as a binary linear programming problem. The authors of [12] introduce the uniform slot allocation and uplink and downlink slot allocation algorithms to achieve higher spatial reuse and network throughput both in uplink and downlink paths.

3. Assumed Topologies

We consider a general form of mesh networks (Figure 1), gird topology, which consists of some SSs and one BS. The length of the network () is assumed to be the square root of the number of nodes ().

Figure 1: General grid topology used for WMNs.

Each one of the nodes has a potential to be the BS; however, two more interesting geographical points will be elaborated in this paper. First, assuming the length of the grid to be an odd number, the BS is considered to be placed at the center of the grid (Figure 2). The main advantage of this topology is its minimized hop count to other nodes. In the second topology (Figure 3), the BS is placed at one of the four corners of the grid. This topology is suitable for special situations such as earthquake-stricken or flood-prone areas, where the BS is placed at the margin of the area and SSs are distributed inside the area. However, with this topology the average hop count to the BS is maximized. We call the former topology BS centered (BSCe) and the later one BS cornered (BSCo).

Figure 2: BS-centered scenario.

Figure 3: BS-cornered scenario.

4. IEEE 802.16d Mesh Mode

Each SS is interconnected with other SSs in its transmission area. This structure forms a graph including some nodes () and some links (). Each time an SS joins the network, it should experience an entry process. This process includes coarse synchronization establishment with the network. From all the potential sponsor nodes that the new node may connect to, it selects a sponsor node having the best signal quality metrics. After joining the network, the new node requests its data rate (service rate) from the BS through other SSs in the path to the BS. The BS collects all bandwidth requests from all SSs and assigns link bandwidths in a fixed frame length. The frame length is fixed at ; since each time its length is changed, the network should be set up again and this imposes some extra overheads to the network. As each node transmits in both uplink and downlink directions, node scheduling is replaced with link scheduling. Since this paper considers only uplink traffic analysis, node and link scheduling are equivalent and interchangeable throughout the paper.

If total duration of bandwidth assignment procedure exceeds the frame length, link durations are fairly scaled down to be fit in the frame length. Scaling down the link durations means that some of the traffic needs to be queued for transmission in the next frames.

A routing algorithm multicasts the end-to-end bandwidth granted to each mesh node by forming a spanning tree in the network. As stated in [13], this spanning tree can be constructed after a new node joins the network. We form one spanning tree for each of the two scenarios shown in Figures 4 and 5. Each tree has some leaf nodes and some middle nodes. Leaf nodes only generate traffic flows; however, middle nodes not only relay packets generated by leaf nodes but also may transmit their own traffic streams.

Figure 4: BS-centered spanning tree.

Figure 5: BS-cornered spanning tree.

In the basic IEEE 802.16d scheduling algorithm, the transmission order of each link (in this paper, node) in the frame is computed during the breadth-first traversal of the routing algorithm. The link transmission times are sequential based on their ranking obtained from the previous stage. Once a transmission ends, the next transmission is started so there is no overlapping or gap between the two consecutive transmissions. This means that the IEEE 802.16 scheduling algorithm does not benefit the spatial reuse mechanism that allows simultaneous transmissions of nonconflicting links.

Link durations are expressed in terms of minislots so that each minislot consists of some OFDM symbols. The maximum number of data minislots in a frame is limited to 256 [7]. Each OFDM symbol carries a number of bits (), which is dependent to the modulation. Therefore, the number of bits transmitted in an OFDM symbol belongs to the , where 96 is the number of bits carried in a symbol with the BPSK-1/2 modulation and 864 is the number of bits carried in a symbol with 64QAM-3/4 [7]. The total bandwidth occupied by the carriers is considered to be fixed at 20 MHz, so the duration of OFDM symbol will be . Assuming that each minislot consists of OFDM symbols and each transmission has an overhead of OFDM symbols, then, for the requested bandwidth of by node , in terms of bits per second, the number of requested minislots is given by

Since may be scaled down due to inadequate space in the frame, the number of granted minislots will be

where is the scaling factor. As spatial reuse may be used in the scheduling algorithm, is greater or equal to 256; however, in this paper spatial reuse is neglected for simplicity; so, . Moreover, we exclude the overhead from analysis (). We assume that ; so .

Granted link durations are fixed unless the bandwidth request of one or more nodes changes. Under such circumstances, the requested bandwidth propagates through the network and the BS recalculates the bandwidth granted to all the nodes again. This condition imposes some extra overheads to the system, which should be avoided [9]. We assume that the link durations or node service rates have been fixed for a long time.

5. Analysis

5.1. Mean Hop Count

Let be a random variable denoting node hop counts to the BS. Probability distribution of node hop counts to the BS can be obtained by means of its probability generation function (PGF):

where is the th hop count and denotes the probability of number of nodes in the th hop count. The mean hop count to the BS is defined as The expected number of hops to the BS in the BSCe scenario () is computed as follows:

With some simple algebra, can be simplified as By using (5.2) to calculate the expected number of hops to the BS in the BSCo scenario (), we have

Again, can be simplified as

Observation 1. The ratio of the expected number of hops to the BS in the BSCo scenario to the expected number of hops to the BS in the BSCe scenario, for large number of nodes in the system, is 2.

Proof. From (5.4) and (5.6), it is clearly observed that

5.2. Mean Service Rate

We assume that the service rate of each node follows a constant distribution with parameter . This assumption is equivalent to what is stated in the IEEE 802.16d central mesh mode described in previous sections.

The service process of the outgoing link of a node can be found by adding up its own service rate and the incoming service rates of that node:

where indicates the service rate of the th link located at the th hop count. indicates the set of all incoming links to node .

The mean service rate of links is obtained by the following relation:

This parameter, for the BSCe and the BSCo topologies, is obtained by the following relations, respectively:

Observation 2. It can be seen from (5.4), (5.6), and (5.10) that:

5.3. Mean Arrival Rate

The arrival rate of each node is assumed to follow a Poisson distribution with parameter . With the spanning trees of Figures 4 and 5, a queuing network is constructed. The arrival rate of the last hop link (farthest link related to the BS) is . Since each link carries not only its own traffic flows but also its subtree traffic flows, the links closer to the BS get more arrival and service rates. Based on the Burke theorem, this arrival rate passes through the next link (because of stability conditions), and, next link, with respect to its own node, has a Poisson arrival rate of . The following relation is obtained for a general topology:

where indicates the arrival rate of the th link located at the th hop count. The mean arrival rate of the links is obtained by the following relation:

The mean arrival rate for the BSCe and the BSCo topologies is obtained as follows, respectively:

Observation 3. It can be seen from (5.4), (5.6), and (5.14) that

5.4. Queuing Model

From the previous sections, it is concluded that the queuing model for the system follows an M/D/1 model for each link. The arrival process follows the Poisson model with parameter and the service model is assumed to be constant and equal to for the link , respectively. The queuing metrics for this model are summarized as follows [14]: where is the queuing delay of the link , is the sojourn time of the link , and represents the queue length of the link with .

6. Virtual Node Concept

In this section, we take advantages of a virtual node () to describe the overall properties of the mesh network, including

(i) the average hop counts to the BS,

(ii) their arrival and service rates.

This could be a candidate node replacing all the nodes other than the BS and is located at hop count from the BS (Figure 6). The arrival and service rates of its connected link are and , respectively, both defined by Observation 2 and Observation 3. The stability condition for this model is

Figure 6: The solution.

The aggregate service rate of all links should be smaller than the maximum service rate of a frame of length 256 minislots. The service rate of a frame can be obtained by

Since there are links in a frame, it can be concluded that

With some algebra, the following relations can be derived to obtain for the BSCe and the BSCo topologies, respectively:

Assuming the total bandwidth usage of the frame, the less than or equal sign is turned to equal sign. So

Since, for the stability conditions, should be less than , we have Each node should be granted at least one minislot per frame; hence: From (6.6), (6.7), and (6.9) it can be concluded that

The average end-to-end waiting time, delivering traffic streams of all the nodes to the BS, is obtained by

which is turned to Calculating , based on the data units, turns (6.12) to In which is obtained by (6.5), or (6.6) and is calculated by (5.4) or (5.6), respectively. It is worth noting that in this case .

6.1. Examples

In this section, two simple examples are used to validate and illustrate the above mentioned model.

6.1.1. Simple Linear Topology

Figure 7 shows the linear topology used for our experience. Each node is fed with the Poisson arrival rate with parameter . The arrival and service rates of each link are illustrated on that link. Each link is assumed to have an infinite buffer size.

Figure 7: Simple linear topology.

Mean end-to-end delay in this topology is calculated as follows:

The subscript in indicates the th . Since ,

The actual mean end-to-end delay in this topology is obtained by

From the solution described in this section, we have

Then, in the model is calculated as

which is equal to .

6.1.2. Simple Tree Topology

In simple tree topology of Figure 8, the waiting times in the system are obtained by the following relations:

Figure 8: Simple tree topology.

which results in

From the model, we haveThen, in the model is calculated as

which is equal to .

6.2. Impact of Scheduling Algorithm (Scheduling Delay)

Observation 4. The maximum number of required frames, in which the BS can process all the nodes, is the height of the spanning tree.

Proof. In a spanning tree the node with maximum number of hops to the root (one of the leaf nodes as an example) has the worst condition to reach the root. This situation gets worse if the transmission order of links is conversely from root to that leaf. In the worst condition, the scheduling algorithm does not exploit the spatial reuse mechanism. So each link of the path to the root, on which the leaf node is scheduled, should be scheduled in a separate frame.

The analytical model described in the previous section is based on the fact that the transmission order of the links is recursively from leaf nodes to the BS. For example in Figure 8, links 2 and 3 start the transmission; then, transmission of link 1 will be commenced.

Each link should be scheduled in a frame such that all the links take some amount of bandwidth in each frame. The transmission order of each link in the frame has a direct influence on delay consideration. For example, suppose that the transmission order of links in the linear topology of Figure 7 follows . Since each link should be scheduled in each frame, all the links deliver their traffic in one frame (Scheduling Algorithm 1 (SA1)). This situation is the best in terms of the total delay for delivering the traffic streams to the BS and has been analyzed in the previous sections. On the other hand, suppose the transmission order follows (Scheduling Algorithm 2 (SA2)) meaning that delivers its data to the BS in the first frame, delivers its data in the second frame, and delivers its data in the third frame. Therefore, the average delay in terms of the number of frames is 2. This situation, on the other hand, is the worst in delivering traffic streams to the BS. All the other situations lie between these two scenarios.

We estimate scheduling delay in terms of the frame length in the solution. As stated in the previous sections, can be obtained by (5.4) or (5.6). This number of hops (in this case, links) can be scheduled in maximum frames, based on Observation 4. In the best case (by using SA1) all of them can be scheduled in one frame; on the other hand, in the worst case (by using SA2) all of them can be scheduled in frame. It is worthy to note that we defined this metric comparing various scenarios; however, the accurate analysis of delay in this case is out of the scope of this paper and in left for another work in this area.

PGF of , a random variable denoting the number of required frames to deliver all the packets of nodes to the BS, is defined as

where denotes the number of required frames and is the probability of frames required to deliver all the packets to the BS. So, the mean number of required frames is defined as

Equation (6.13), with respect to (6.5) and (6.6) and with some simple algebra, can be rewritten as follows:

Separating queuing and service delay results in

The first term in the right-hand side of (6.26) is queuing delay and the second term is serving delay of each unit of data in the solution. This equation implies that end-to-end delay in both topologies of Figures 4 and 5 is equal; since their arrival and service rates are different based on (6.5) and (6.6); however, this is only true in the case of SA1 deployment. Providing such a relation in the case of SA2 deployment is out of the scope of this paper; however, it is possible to estimate this metric as well.

The serving delay in (6.26) can be expressed as scheduling delay. With respect to the solution and Observation 4, it can be concluded that, in the case of SA2, the scheduling delay is turned to . So

Finding is out of the scope of this paper; however, we can put a maximum bound on for the stability condition as well.

Since must be less than 1 and the serving delay is equal to , we have

7. Performance Evaluation

The system described in the previous sections has two essential delay constraints: queuing delay () and scheduling delay (). The mean end-to-end delay of the system is given by

To evaluate the performance of the system, the two scheduling algorithms, SA1 and SA2, which are the best and the worst cases, respectively, have been deployed to the analysis. is computed in these two cases by the following relations:

7.1. Delay Analysis

Deploying SA1 and SA2 scheduling algorithms results in the following relations obtaining the end-to-end delay:

Observation 5. Assuming the same work load in the cases SA1 and SA2, it can be concluded that

We set the parameters , , and to 3, , and 10 ms, respectively. At first, we compare queuing and scheduling delay constraints. Figure 9 shows this comparison for the case of SA1. It can be seen that for less than the scheduling delay is dominant. When the amount of work load passes , the queuing delay will be the main factor and for work loads close to 1 the scheduling delay can be neglected.

Figure 9: Queuing and scheduling delay comparison using SA1 scheduling algorithm.

Figure 10 compares the scheduling and queuing delay constraints for the case SA2 while the BSCe topology is exploited. In this case the queuing delay is dominant for greater than , too.

Figure 10: Queuing and scheduling delay comparison using SA2 scheduling algorithm for the BSCe topology.

Figures 11 and 12 illustrate the end-to-end delay for the case SA2 with possible .

Figure 11: End-to-end delay comparison using SA2 scheduling algorithm for the BSCe topology.

Figure 12: End-to-end delay comparison using SA2 scheduling algorithm for the BSCo topology.

Next, we compare end-to-end delay of the two topologies, the BSCe and the BSCo, in the cases of SA1 and SA2. Figure 13 compares end-to-end delay for the four scenarios.

Figure 13: End-to-end delay comparison for the four scenarios.

The main results derived from Figure 13 are summarized as follows.(i)In the case of SA2, delay of the BSCe topology for is greater than that of the BSCo topology for .(ii)In the case of SA2, delay of the BSCe topology for is greater than that of the BSCo topology for .(iii)In the case of SA2, delay of the BSCo topology and that of the BSCe topology for , and , respectively, are such that they are not suitable for real-time applications (maximum delay bound for real-time applications is about 150 ms).(iv)The BSCe topology with SA2 deployment and is suitable for real-time applications; however, the BSCo topology with SA2 deployment and is suitable for real-time applications for the network lengths of less than 4.(v)Delay bound for the two topologies and is greater than that for the BSCe topology with SA2 scheduling algorithm and .

7.2. Throughput Evaluation

Since each individual node is granted minislots, the throughput of the system is calculated according to the following relation, in the case of SA1, by means of the concept: In the case of SA2 deployment, the system throughput is computed as follows:

As previously stated, is equal to ; so:

Observation 6. Assuming the same work load in the two cases SA1 and SA2, it can be concluded that

Figure 14 illustrates the throughput comparison of the four scenarios regarding three values of for each scenario. The main results derived from the figure are summarized as follows. (i)As it can be indicated from Figure 13, the delay bounds for the two topologies in the case of SA1 are equivalent; however, the throughput of the BSCo topology is less than that of the BSCe topology in the case of SA2. This degradation is because of lower data rates of the BSCo topology.(ii)The throughput of the SA2 scheduling algorithm is greater than the throughput of the SA1 scheduling algorithm for small network lengths; however, this measure is lower as the length of the network grows. Since the network length for the BSCe topology could not be set to 2, this situation occurs.

Figure 14: Throughput of the system comparison for the four scenarios.

As the throughput is dependent on the modulation used, Figure 15 shows this dependency for for the BSCe topology in the case of SA1. It can be indicated from the figure that the maximum throughput reaches 39.8131 when and ; this measure is 1.8959 when and .

Figure 15: Throughput of the system comparison using SA1 scheduling algorithm for the BSCe topology.

8. Simulation Results

We have developed a simulation scenario in the NS-2 simulation environment. It is assumed that all links exploit the BPSK-1/2 modulation ( and each minislot contains 3 OFDM symbols. So, minislot size has been fixed at bits. The buffer size of each node has been set so large such that no packet loss is occurred. Transmission and interference ranges of each node are equal with respect to the modulation used.

For each network length we ran the simulation about 50 times and set the average as the simulation result. We took two scenarios for the simulation. The work load of the system in these two scenarios has been set to 0.9. In the first scenario, the BSCe topology along with SA1 scheduling algorithm with has been considered. In the second scenario, the BSCo topology along with SA2 scheduling algorithm with has been considered. In both cases, the analytical model has been illustrated along with the simulation result to verify the analytical model.

Figure 16 illustrates end-to-end delay calculated through simulations for various lengths of the network. As it can be indicated from the figure, the simulation and analytical results follow the same curve in the two cases; however, for the case of simulation because of some extra overheads, the end-to-end delay is more than that of analytical results. As network length grows up, the overhead gets bigger and the simulation results get a little far from the analytical approach. As a result, it can be seen that the analytical model can predict the system delay as well.

Figure 16: End-to-end delay comparison of the analytical and simulation approaches (, , , and ).

Figure 17 shows the throughput calculated through simulations for various lengths of the network. As it can be seen in this figure, the throughput degradation is completely evident in the case of simulation results. A simple indication is that, as each minislot carries some overhead, the aggregated minislots flown in the network are smaller than those of the analytical model. Again, the analytical results follow the same curves as simulation results and it can be concluded that the analytical approach predicts the throughput of the system as well.

Figure 17: Throughput comparison of analytical and simulation approaches (, , , and ).

9. Conclusions

The main contribution of this paper was devoted to introducing a model for estimating the performance measures of WiMAX-like WMNs. Exploiting the benefits of the virtual node solution and introducing two types of scheduling algorithms, delay and throughput metrics of the two topologies were obtained. The analytical study shows that the throughput of the system is highly dependent on the position of the BS and the length of the network. Although the analytical approach was restricted by some simplifying assumptions, its verification with simulation results confirms that this approach estimates the system parameters well.

References

IEEE Std 802.16-2004, “IEEE standard for local and metropolitan area networks part 16: air interface for fixed broadband wireless access systems,” 2004

G. Min, Y. Wu, K. Li, and A. Y. Al-Dubai, “Performance modelling and optimization of integrated wireless LANs and multi-hop mesh networks,” International Journal of Communication Systems, vol. 23, no. 9-10, pp. 1111–1126, 2010.

Z. Liu, M. Yang, and J. Dai, “Performance improvement based on path delay analysis in WiMax mesh networks,” in Proceedings of the 2nd International Conference on Communications and Networking in China (ChinaCom '07), pp. 958–962, Shanghai, China, August 2007.

Y. Feng, X. Shen, Z. Gao, and G. Dai, “Queuing based traffic model for wireless mesh networks,” in Proceedings of the 15th International Conference on Parallel and Distributed Systems (ICPADS '09), pp. 648–654, Shenzhen, China, December 2009.

I. F. Akyildiz, X. Wang, and W. Wang, “Wireless mesh networks: a survey,” Computer Networks, vol. 47, no. 4, pp. 445–487, 2005.

Y. Liu, T. Wu, and X. Xian, “Queue control based delay analysis and optimization for wireless mesh networks,” in Proceedings of the 2nd International Conference on Education Technology and Computer (ICETC '10), pp. 104–108, Shanghai, China, June 2010.

P. Djukic and S. Valaee, “Scheduling algorithms for 802.16 mesh networks,” in WiMax/MobileFi: Advanced Research and Technology, Y. Xiao, Ed., pp. 267–288, Auerbach Publications, 2007.

M. Kas, B. Yargicoglu, I. Korpeoglu, and E. Karasan, “A survey on scheduling in IEEE 802.16 mesh mode,” IEEE Communications Surveys and Tutorials, vol. 12, no. 2, pp. 205–221, 2010.

P. Djukic and S. Valaee, “Link scheduling for minimum delay in spatial re-use TDMA,” in Proceedings of the 26th IEEE International Conference on Computer Communications (IEEE INFOCOM '07), pp. 28–36, Anchorage, Alaska, USA, May 2007.

H.-Y. Wei, S. Ganguly, R. Izmailov, and Z. Haas, “Interference-aware IEEE 802.16 WiMax mesh networks,” in Proceedings of the IEEE 61st Vehicular Technology Conference (VTC '05), pp. 3102–3106, Stockholm, Sweden, June 2005.

J. Zou and D. Zhao, “Real-time CBR traffic scheduling in IEEE 802.16-based wireless mesh networks,” Wireless Networks, vol. 15, no. 1, pp. 65–72, 2009.

S. Liu, S. Feng, W. Ye, and H. Zhuang, “Slot allocation algorithms in centralized scheduling scheme for IEEE 802.16 based wireless mesh networks,” Computer Communications, vol. 32, no. 5, pp. 943–953, 2009.

Y. Zhang, H. Hu, and H.-H. Chen, “QoS differentiation for IEEE 802.16 WiMAX mesh networking,” Mobile Networks and Applications, vol. 13, no. 1-2, pp. 19–37, 2008.

L. Kleinrock, Queueing Systems: Volume 1: Theory, Wiley, New York, NY, USA, 1975.