# Product-Form Solutions for Integrated Services Packet Networks and Cloud Computing Systems

Authors: Wanyang Dai
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We iteratively derive the product-form solutions of stationary distributions for a type of preemptive priority multiclass queueing networks with multiserver stations. This type of queueing systems can typically be used to model the stochastic dynamics of some large scale backbone networks with multiprocessor shared-memory switches or local (edge) cloud computing centers with parallel-server pools. The queueing networks are Markovian with exponential interarrival and service time distributions. The obtained iterative solutions can be used to conduct performance analysis or as comparison criteria for approximation and simulation studies. Numerical comparisons with existing Brownian approximating model (BAM) related to general interarrival and service times are provided to show the effectiveness of our current designed algorithm and our previous derived BAM. Furthermore, based on the iterative solutions, we can also give some analysis concerning network stability for some cases of these queueing systems, which provides some insight for more general study.


---

## Body

## 1. Introduction

At present, integrated services packet networks (ISPN) are widely used to transport a wide range of information such as voice, video, and data. It is foreseeable that this integrated services pattern will be one of the major techniques in the future cloud computing based communication systems and Internet (see, e.g., Mullender and Wood [1]). The introduction of concept and architecture of cloud computing can be found in, for example, Mell and Grance [2] and Rhoton and Haukioja [3]. A possible cloud computing based telecommunication network architecture (i.e., a large-scale network infrastructure as a service) is designed in Figure 1, where an end-user may require service (or services) from single local cloud computing center or multiple local and remote cloud computing centers. Among these centers, they communicate each other by using core switching network system (note that the switching system itself can also be independently viewed and handled as a cloud computing system with multiple service pools).Figure 1
An integrated services cloud computing based network.The speed and efficiency of core switching network systems are the bottleneck in realizing high speed owing to the drastic improvement in transmission speed and reliability of optical fiber. In an ISPN network, information is partitioned into packets depending on the employed protocol such as Internet protocol (IP). For the purpose of transmission, each packet consists of user’s data payload to be transmitted, a header containing control information (e.g., source and destination addresses, packet type, priority, etc.), and a checksum used for error control. The high-speed ISPN networks require fast packet switches to move packets along their respective virtual paths. The switches are computers with processing and storage capabilities. The main objective of a switch is to route an incoming packet arriving along a particular input link to a specified output link. More precisely, once the incoming packet is entirely received, the switch examines its header and selects the next outgoing link. In other words, packets are transmitted in astore-and-forward manner.Various pieces of information can be classified into a fixed number of different types with separate requirements of quality of service (e.g., different end-to-end delays and packet loss ratios). Real time traffic packets with stringent delay constraint (e.g., interactive audio/video) are endowed service priority. In the meanwhile, it is imperative to size buffers for nonreal time traffic packets (e.g., data), which can tolerate higher delay but demand much smaller packet loss ratios. Hence, efficient switching and buffer management schemes are needed for switches.Currently, three basic techniques (see, e.g., Tobagi [4] and Rao and Tripathi [5]) are designed to realize the switching function:space-division,shared-medium, andshared-memory. In our switches, the shared-memory technique is employed, which is comprised of a single dual-ported memory shared by all input and output lines. Packets arriving along input lines are multiplexed into a single stream that is fed to the common memory for storage. Inside the memory, packets are organized into different output buffers, one for each output line. In each output buffer, memory can be further divided into priority output queues, one for each type of traffic. In the meantime, an output stream of packets is formed by retrieving packets from the output queues sequentially, one per queue. Among different traffic types, packets are retrieved according to their priorities. Inside each type, packets are retrieved under the first-in first-out (FIFO) service discipline. Then, the output stream is demultiplexed and packets are transmitted on the output lines. There are some other ways to deal with the output queues such as processor sharing among output lines (see, e.g., Tobagi [4] and Rao and Tripathi [5]). We will address these issues elsewhere (see, e.g., Dai [6]).In addition to switching, queueing is another main functionality of packet switches. The introduction of queueing function is owing to multiple packets arriving simultaneously from separate input lines and owing to the randomness of packet arrivals from both outside routing and internal routing of the network. There are three possibilities for queueing in a switch:input queueing,output queueing, andshared-memory queueing (see, e.g., Schwarz [7] and Rao and Tripathi [5]). Shared-memory-queueing mechanism is employed in our switches since it has numerous advantages with respect to switching and queueing over other schemes. For example, both switching and queueing functions can be implemented together by controlling the memory read and write properly. Furthermore, modifying the memory read and write control circuit makes the shared-memory switch sufficiently flexible to perform functions such as priority control and other operations.Nevertheless, no matter what technology is employed in implementing the switch, it places certain limitations on the size of the switch and line speeds. Currently, two possible ways can be used to build large switches to match the transmission speed of optical fiber. The first one is to adopt parallel processors (see, e.g., Figure2). The second one is to interconnect many switches (known as switch modules) in amultistage configuration to build large switches (see, e.g., Figure 3). The remaining issue is how to reasonably allocate resources of these switches and efficiently evaluate the system performance.Figure 2
A two-type and four-class network with two multiprocessor switches. Type 1 includes class 1 and class 2. Type 2 consists of class 3 and class 4. Type 1 owns the service priority.Figure 3
A two-stage switch with four dual-ported shared-memory switching modules.The statistical characteristics of packet interarrival times and packet sizes have a major impact on switch hardware and software designs owing to the consideration of network performance. How to more effectively identify packet traffic patterns is a very active and involved research field (see, e.g., Nikolaidis and Akyildiz [8] and Dai [6, 9, 10]). Independent and identical distribution (i.i.d.) is the popular assumption for these times and packet sizes. Doubly stochastic renewal process introduced in Dai [6] is the latest definition and generalization related to input traffic and service processes for a wireless network under random environment. The effectiveness of these characteristics is supported by recent discoveries in Cao et al. [11] and Dai [9, 10].In all circumstances, it is imperative to find product-form solutions for those queueing networks under suitable conditions to conduct performance analysis or provide comparison criteria to show the effectiveness of approximation and/or simulation studies (see, e.g., Dai et al. [6, 12–15]). Furthermore, note that the stochastic dynamics of the backbone networks with multiprocessor shared-memory switches and the local (edge) cloud computing centers with parallel-server pools in Figure 1 can both be described by a multiclass queueing networks with parallel servers at each station. Hence, in this paper, we aim to derive the product-form solutions iteratively for one particular type of related queueing networks, that is, a type of particular preemptive priority multiclass queueing networks when the interarrival and service times are exponentially distributed. In addition, we also aim to provide some numerical comparisons with existing Brownian approximating model (BAM) related to general interarrival and service times to show the effectiveness of our current designed algorithm and our former derived BAM.Next, we provide some review of the existing literature associated with the current study. Comparisons between the existing studies and our current discussion are also presented.Under a general Whittle network framework, the product-form solutions are presented in Serfozo [16] for some multiclass networks, which include those with sector-dependent (e.g., Example 3.3) and class-station-dependent service rates (e.g., BCMP networks in Section 3.3, which are introduced by Baskett et al. [17]). Without considering the interaction among different stations, the distinguishing feature of these networks is that a job’s service rate at a station may consist of two intensities: one referred to asstation intensity is a function of the total queue length at the station and the other one referred to asclass intensity is a function of the queue length in the same class as the job being served. Some specific single-class queueing systems of BCMP networks in Section 3.3 of Serfozo [16] are revisited in Harrison [18, 19] by developing some Reversed Compound Agent Theorem (RCAT) based method. The corresponding product-form and non-product-form solutions are derived. Nevertheless, as claimed in Harrison [18, 19], heavier notation will be involved as long as a multiclass BCMP network is concerned.Although our networks are of a form of those networks with multiple types of units as introduced in Serfozo [16], our networks are beyond those with sector-dependent and class-station-dependent service rates as introduced previously. More precisely, for a station in our networks, the station intensity is not only a function of the total queue length but also a function of combinations of queue lengths of various classes. The class intensity depends on not only the queue length of itself and/or the total queue length but also the numbers of jobs in other classes at the station. Therefore, how to find suitable function Φ and how to prove our service rates to be Φ-balanced as defined in Chapter 3 of Serfozo [16] are not obvious. Moreover, how to apply the RCAT based method developed in Harrison [18, 19] to our multiclass network cases is also not trivial. Thus, in this paper, we use the method of solving Kolmogorov (balance) equations to get the product-form solutions iteratively, which are more engineering and computationally friendly. Furthermore, by this method, we can also give some analysis concerning network stability for some cases of these systems, which provides some insight for more general study.The remainder of this paper is organized as follows. The open priority multiclass queueing network associated with high-speed ISPN is described in Section2. Our main results including product-form solutions and performance comparisons are presented in Section 3. Numerical comparisons are given in Section 4. The proofs of our main theorems are provided in Section 5. The conclusion of the paper is presented in Section 6.


## 2. The Queueing Network Model

Note that the stochastic dynamics of the backbone core networks with multiprocessor shared-memory switches and the local (edge) cloud computing centers with parallel-server pools in Figure1 can both be described by multiclass queueing networks with parallel servers at each station. Therefore, we consider a queueing network that has J multiserver stations. Each station indexed by j∈{1,…,J} owns cj servers and has an infinite capacity waiting buffer. In the network, there are I job types. Each type consists of J job classes that are distributed at different stations. Therefore, the network is populated by K(=IJ) job classes which are labeled by k∈{1,…,K}. Upon the arrival of a job of a type from outside the network, it may only receive service for part of J classes and may visit a particular class more than once (but at most finite many times). Then, it leaves the network (i.e., the network is open). At any given time during its lifetime in the network, the job belongs to one of the classes and changes classes at each time a service is completed. All jobs within a class are served at a unique station and more than one class might be served at a station (so-calledmulticlass queueing network). The ordered sequence of classes that a job visits in the network is nameda route. Interrouting among different job types is not allowed throughout the entire network.We useC(j) to denote the set of classes belonging to station j. Let s(k) denote the station to which class k belongs. We implicitly set j=s(k) when j and k appear together. Associated with each class k, there are two i.i.d sequences of random variables (r.v.), uk={uk(i),i≥1} and vk={vk(i),i≥1}, an i.i.d sequence of K-dimensional random vectors, ϕk={ϕk(i),i≥1}, and two real numbers, αk≥0 and mk=1/μk>0. We suppose that the 3K sequences
(1)u1,…,uK,v1,…,vK,ϕ1,…,ϕK
are mutually independent. The initial r.v.s uk(1) and vk(1) have means E[uk(1)]=1/αk and E[vk(1)]=mk, respectively. For each i∈{1,2,…}, uk(i) denotes the interarrival time between the (i-1)th and the ithexternally arrival job at class k. Furthermore, vk(i) denotes the service time for the ith class k job. In addition, ϕk(i) denotes the routing vector of the ith class k job. We allow αk=0 for some classes k∈E≡{k:αk≠0}. Then, it follows that αk and μk are the external arrival rate and service rate for class k, respectively. We assume that the routing vector ϕk(i) takes values in {e0,e1,…,eK}, where e0 is the K-dimensional vector of all 0’s. For l=1,…,K, el is the K-dimensional vector with lth component 1 and other components 0. When ϕk(i)=el, the ith job departing class k becomes a class l job. Let pkl=P{ϕk(i)=el} be the probability that a job departing class k becomes a class l job (of the same type). Thus, the corresponding K×K matrix P=(pkl) isrouting matrix of the network. Furthermore, the matrix
(2)Q=I+P′+(P′)2+⋯
is finite; that is, (I-P′) is invertible with Q=(I-P′)-1 since the network is open. The symbol ′ on a vector or a matrix denotes the transpose and I denotes the identity matrix.We useλk for k∈{1,…,K} to denote the overall arrival rate of class k, including both external arrivals and internal transitions. Then, we have the followingtraffic equation:
(3)λk=αk+∑l=1K‍λlplk,
or in its vector form: λ=α+P′λ (all vectors in this paper are to be interpreted as column vectors unless explicitly stated otherwise). Note that the unique solution λ of (3) is given by λ=Qα. For each λk, if there is a long-run average rate of flow into the class which equals the long-run average rate out of that class, this rate will equal λk. Furthermore, we define thetraffic intensityρj for station j as follows:
(4)ρj=∑k∈C(j)‍λkcjμk,
where ρj with ρj≤1 is also referred to as the nominal fraction of time that station j is nonidle.The order of jobs being served at each station is dictated by a service discipline. In the current research, we restrict our attention to static buffer priority (SBP) service disciplines under which the classes at each station are assigned a fixed rank (with no ties). In our queueing network, each type of jobs is assigned the same priority rank at every station where it possibly receives service. When a server within a station switches from one job to another, the new job will be taken from the leading (or longest waiting) job at the highest rank nonempty class at the server’s station. Within each class, it is assumed that jobs are served on the first-in first-out (FIFO) basis. We suppose that the discipline employed isnonidling; that is, a server is never idle when there are jobs waiting to be served at its station. We also assume that the discipline ispreemptive resume; that is, when a job arrives at a station with all servers busy and if the job is with a higher rank than at least one of the jobs currently being served, one of the lower rank job services is interrupted; when there is a server available to the interrupted service, it continues from where it left off. For convenience and without loss of generality, we use consecutive numbers to index the classes that have the same priority rank at stations 1 to J. In other words, the highest priority classes for station 1 to J are indexed by 1 to J, the second highest priority classes are indexed by J+1 to 2J,…, and the lowest priority classes are indexed by K-J+1 to K. An example of such a two-station network is given in Figure 2. In the network, type 1 traffic possibly requires class 1 and class 2 services, type 2 traffic possibly requires class 3 and class 4 services, and classes in type 1 have the higher priority at their corresponding stations.Finally, we define the cumulative arrival, cumulative service, and cumulative routing processes by the following sums:(5)Uk(n)=∑i=1n‍uk(i),Vk(n)=∑i=1n‍vk(i),Φk(n)=∑i=1n‍ϕk(i),
where the lth component Φlk(n) of Φk(n) is the cumulative number of jobs to class l for the first nth jobs leaving class k with n=1,2,… and k,l∈{1,…,K}. Then, we define
(6)Ek(t)≡max⁡{n≥0,Uk(n)≤t},Sk(t)≡max⁡{n≥0,Vk(n)≤t},Ak(t)≡Ek(t)+∑l=1K‍Φkl(Sl(t)).
Note that Ek(t) denotes the total number of external arrivals to class k in time interval [0,t]. Sk(t) represents the total number of class k jobs which have finished service in [0,t]. Ak(t) is the total arrivals to class k in [0,t] including both external arrivals and internal transitions.


## 3. Steady-State Queue Length Distributions

We useQk(t) to denote the number of class k jobs in station j=s(k) with j∈{1,2,…,J} and k∈{1,2,…,K} at time t. It is called thequeue length process for class k jobs. For convenience, let Q(i,j)(t) and n(i,j) be the (j-i+1)-dimensional queue length process and (j-i+1)-dimensional state vectors, respectively. They are given by
(7)Q(i,j)(t)=(Qi(t),…,Qj(t)),n(i,j)=(ni,…,nj)
for j>i and i,j∈{1,2,…,K} and nonnegative integers ni,…,nj. Then, we use Pn(1,K)(t) to denote the probability of system at state n(1,K)=(n1,…,nK) and at time t; that is,
(8)Pn(1,K)(t)=P{Q(1,K)(t)=n(1,K)}.
Furthermore, let
(9)Pn(1,K)=P{Q(1,K)(∞)=n(1,K)}
denote the corresponding steady-state probability of system at state (n1,…,nK) if the network exists in a stationary distribution. In addition, let PniJ+ks(iJ+k) denote the probability at state niJ+k for class iJ+k with i∈{0,1,…,I-1} and k∈{1,2,…,J}.Under the usual convention, letx∧y denote the smaller one of any two real numbers x and y. Let x∨y denote the larger one of x and y; that is, x∧y≡min⁡⁡{x,y} and x∨y≡max⁡{x,y}. Then, for each k∈{1,…,K}, we have the following notation:
(10)a(x,y)≡(x∧(cs(k)-y))∨0.
Finally, for i>0, define
(11)κ(niJ+k)≡{0ifniJ+k=0,∑∑u=0i-1nuJ+k<cs(k)‍a(niJ+k,∑u=0i-1‍nuJ+k)×∏u=0i-1‍PnuJ+ks(uJ+k)ifniJ+k>0.
Furthermore, for i=0, define κ(nk)≡nk∧cs(k).Theorem 1 (steady-state distribution).
Assume that all of the service times and external interarrival times are exponentially distributed with rates as before and the traffic intensityρj<1 for all j∈{1,…,J}. Furthermore, suppose that
(12)λiJ+kκ(niJ+k)μiJ+k<1
for niJ+k>0, i∈{0,…,I-1}, and k∈{1,…,J}. Then, for each h∈{1,2,…,I}, the steady-state distribution is given by the following product form:
(13)Pn(1,hJ)=∏i=0h-1∏k=1JPniJ+ks(iJ+k).
More precisely, for niJ+k≥1,
(14)PniJ+ks(iJ+k)=P0s(iJ+k)∏riJ+k=1niJ+kλiJ+kκ(riJ+k)μiJ+k‍,
and, for niJ+k=0, the initial distribution P0s(iJ+k) is determined by
(15)P0s(iJ+k)=11+∑niJ+k=1∞∏riJ+k=1niJ+k(λiJ+k/κ(riJ+k)μiJ+k).The following theorem relates condition (12) to primitive interarrival time and service rates for some cases of these systems, which provides some insight for more general study.Proposition 2 (network stability condition).
Under the exponential assumptions as stated in Theorem1, if ρj<1, for each j∈{1,…,J}, the stability condition (12) holds for the following networks.
Net I. Multiclass networks with single-server stations, that is, the number cj of servers is one for all stations j∈{1,…,J} while the number I of job types can be arbitrary.
Net II. The number cj of servers can be arbitrary for all stations j∈{1,…,J} while the number of job types equals two (I=2).We conjecture thatρj<1 for j∈{1,…,J} implies the condition (12) for our general network. Nevertheless, owing to complex computation involved, the corresponding analytical illustration is not a trivial task.Example 3.
Consider a network with three job types (I=3) and at least one station j=s(k) having two servers (cs(k)=2 for k∈{1,…,J}) while other stations have at most two servers (cs(l)≤2 for l∈{1,2,…,J}∖{k}). For a station s(k) with three servers, the condition (12) can be explicitly expressed as follows:
(16)((1+λkμk+λJ+k2μJ+k-(λk2μk)2-14(λkμk)3+λkλJ+k4μkμJ+k)λ2J+kμ2J+kk×(1+λkμk+λJ+k2μJ+k-(λk2μk)2-14(λkμk)3+λkλJ+k4μkμJ+k))×(2(1-λk2μk-λJ+k2μJ+k)kkkklk×(1+λkμk+λJ+k2μJ+k-(λk2μk)2-14(λkμk)3kkkkkkkk+2(1-λk2μk)λkλJ+k4μkμJ+k))-1<1
for n2J+k≥2. Under ρj<1, it is easy to see that the above inequality is true if λk/μk≤1. Numerical tests in Table 1 have also been conducted and show that the inequality is true even for 1<λk/μk<2, but the corresponding analytic demonstration could be nontrivial. The detailed illustration of the example will be provided at the end of this paper.Table 1
Numerical tests for network stability condition (12).
1
<
λ
k
/
μ
k
<
2, ρs(k)<1
λ
2
J
+
k
/
κ
(
n
2
J
+
k
)
μ
2
J
+
k
λ
k
/
μ
k
λ
J
+
k
/
μ
J
+
k
λ
2
J
+
k
/
μ
2
J
+
k
0.16748
1.10
0.30
0.10
0.98814
1.10
0.30
0.59
0.99111
1.50
0.30
0.19
0.85890
1.90
0.05
0.04Remark 4.
The justifications of Theorem1 and Proposition 2 are postponed to Section 5. Instead, we will first use these results as comparison criteria to illustrate the effectiveness of the diffusion approximation models developed in Dai [13] and answer the question on when these approximation models can be employed.


## 4. Numerical Comparisons

First of all, we note that partial results presented in this section were briefly reported in the short conference version (Dai [20]) of this paper. More precisely, we consider a network with single-server station and under preemptive priority service discipline. By employing the exact solutions developed in previous sections, we conduct performance comparisons between these product-form solutions and the approximating ones of Brownian network models.Brownian network models, which are also known as semimartingale reflecting Brownian motions (SRBM), have been widely employed as approximating models for multiclass queueing networks with general interarrival and service time distributions when the traffic intensity defined in (4) is suitably large or close to one (see, e.g., Dai [13]). The effectiveness of the Brownian network models has been justified in numerical cases by comparing their solutions with either exact solutions or simulation results (see, e.g., Dai [13], Chen and Yao [21], and Dai et al. [22]). Thus, it is meaningful to illustrate the correctness of our newly derived formula by comparing it with the Brownian network models.For the network with exponential interarrival and service time distributions, by using Theorem1, we get the steady-state mean queue length for each class k∈{1,…,K} as
(17)EQk(∞)=∑nk=0∞‍nkPnks(k)=λkmk1-∑i=1kλimi=αkmk/(1-pkk)1-∑i=1kαimi/(1-pii)
and the expected total time (sojourn time) a job had to spend in the system as
(18)ETk(∞)=mk1-∑i=1kλimi=mk1-∑i=1kαimi/(1-pii).For the network with general interarrival and service time distributions, owing to the nature of our network routing structure, the higher priority classes are independent of lower ones. Then, it follows from the studies in Dai [13], Harrison [23], and Chen and Yao [21] that the steady-state mean sojourn time and mean queue length for each class k can iteratively be calculated with respect to priority order as follows:
(19)ETk(∞)=EWk(∞)+mk1-∑i=1k-1λimi,EQk(∞)=λkETk(∞),
where EWk(∞) is the steady-state mean total workload for all classes i∈{1,…,k} and k=1,…,K. More precisely, it is given by
(20)EWk(∞)=σk22|θk|,
where
(21)θk=(1-pkk)(∑i=1k‍λimi-1),σk2=(1-pkk)2∑i=1k‍(λibi+mi2(αi3ai+λipii(1-pii))(1-pii)2),
and ai, bi are the variances of interarrival and service time sequences for each class i.In our numerical comparisons, we consider an exponential network withK=3. For this case, our corresponding data are listed in Table 2. In the table, ERROR=SRBM-EXACT and RATIO=(|ERROR|/EXACT)*100%, EQ^(∞)=∑i=13‍EQi(∞) and ET^k(∞)=∑i=1k‍Ti(∞) for k=1,2,3. From the table, we can see that SRBMs are more reasonable approximations of their physical queueing counterparts when traffic intensities for the lowest priority jobs are relatively large.Table 2
Performance comparisons for a priority multiclass network with three job types.
K
=
3, m1=1, m2=2, m3=3, p11=0.1, p22=0.2, p33=0.3
λ
1
λ
2
λ
3
ρ
E
Q
^
(
∞
)
E
T
^
1
(
∞
)
E
T
^
2
(
∞
)
E
T
^
3
(
∞
)
EXACT
0.20
0.20
0.06
0.78
2.0682
1.2500
6.2500
14.8864
SRBM
0.20
0.20
0.06
0.78
2.2410
1.2500
6.4722
17.0253
ERROR
−0.1728
0.0000
0.2222
2.1389
RATIO
8.35%
0.00%
3.56%
14.37%
EXACT
0.20
0.20
0.085
0.855
3.0086
1.2500
6.2500
21.9397
SRBM
0.20
0.20
0.085
0.855
3.2349
1.2500
6.4722
24.0785
ERROR
−0.2263
0.0000
0.2222
2.1389
RATIO
7.52%
0.00%
3.56%
9.75%
EXACT
0.20
0.20
0.115
0.945
7.5227
1.2500
6.2500
55.7955
SRBM
0.20
0.20
0.115
0.945
7.8131
1.2500
6.4722
57.9344
ERROR
−0.2904
0.0000
0.2222
2.1389
RATIO
3.86%
0.00%
3.56%
3.83%
EXACT
0.10
0.15
0.115
0.745
1.9641
1.1111
4.4444
12.8758
SRBM
0.10
0.15
0.115
0.745
2.0944
1.1111
4.5432
13.8804
ERROR
−0.1303
0.0000
0.0988
1.0046
RATIO
6.64%
0.00%
2.22%
7.80%
EXACT
0.10
0.15
0.15
0.85
3.6111
1.2500
4.4444
21.1111
SRBM
0.10
0.15
0.15
0.85
3.7766
1.2500
4.5432
22.1157
ERROR
−0.1655
0.0000
0.0988
1.0046
RATIO
4.58%
0.00%
2.22%
4.76%
EXACT
0.10
0.15
0.185
0.955
12.9445
1.2500
4.4444
67.7778
SRBM
0.10
0.15
0.185
0.955
13.1451
1.2500
4.5432
68.7824
ERROR
−0.2007
0.0000
0.0988
1.0046
RATIO
1.55%
0.00%
2.22%
1.48%
EXACT
0.45
0.22
0.02
0.95
6.0182
1.8182
20.0000
61.8182
SRBM
0.45
0.22
0.02
0.95
6.3818
1.8182
20.7273
72.0000
ERROR
−0.3636
0.0000
0.7273
10.1818
RATIO
6.04%
0.00%
3.64%
16.47%
EXACT
0.45
0.22
0.03
0.98
9.3182
1.8182
20.0000
151.8181
SRBM
0.45
0.22
0.03
0.98
9.7836
1.8182
20.7273
161.9999
ERROR
−0.4655
0.0000
0.7273
10.1818
RATIO
5.00%
0.00%
3.64%
6.71%


## 5. Proofs of Theorem1 and Proposition 2

### 5.1. Proof of Theorem1

For convenience, we introduce some additional notations. LetΔEk(t), ΔSk(t), and ΔAk(t) be defined by
(22)ΔEk(t)≡Ek(t+Δt)-Ek(t),ΔSk(t)≡Sk(t+Δt)-Sk(t),ΔAk(t)≡Ak(t+Δt)-Ak(t),
which denote the cumulative external arrivals to class k, the cumulative number of jobs finished services at class k, and the total arrivals to class k in [t,t+Δt]. Then, we can justify Theorem 1 by induction as in the following several steps.Step 1.
We consider the steady-state distribution for the highest rank classes with each indexk∈{1,2,…,J}. In this case, the type index i=0 is used in Theorem 1. Owing to the preemptive-resume service discipline and the class routing structure, we know that these J classes form a Jackson network. Then, by the theorem in Jackson [24], we have the following product form:
(23)Pn(1,J)=Pn1s(1)Pn2s(2)⋯PnJs(J),
where Pnks(k) denotes the steady-state probability at state nk for class k(∈{1,2,…,J}) at station s(k). More precisely, it is given by
(24)Pnks(k)={P0s(k)(λk/μk)nknk!,if0≤nk<ck,P0s(k)(λk/μk)nknk!(nk)nk-ck,ifnk≥ck
with P0s(k) being determined by equation ∑nk=0∞‍Pnks(k)=1; that is,
(25)P0s(k)=(∑nk=0cs(k)-1‍1nk!(λkμk)nk+(λk/μk)cs(k)cs(k)!(1-λk/cs(k)μk))-1.Step 2.
We derive the steady-state distribution for the highest rank and the second highest rank classes withk∈{1,2,…,2J}. Owing to the preemptive-resume service discipline and the job routing structure, we have
(26)Pn(1,2J)(t)=Pn(1,J)/n(J+1,2J)(t)Pn(J+1,2J)(t)=Pn(1,J)(t)Pn(J+1,2J)(t),
where the second equality follows from the preemptive assumption and Pn(1,J)/n(J+1,2J)(t) is the conditional probability at time t for classes with k∈{1,…,J} at state n(1,J) in terms of classes with k∈{J+1,…,2J} at state n(J+1,2J); that is,
(27)Pn(1,J)/n(J+1,2J)(t)=P{Q(1,J)(t)=n(1,J)∣Q(J+1,2J)(t)=n(J+1,2J)}.
In order to get the steady state distribution for Pn(J+1,2J)(t), we consider each state n(J+1,2J) at time t for the second highest rank class jobs. There are several ways in which the system can reach it. They can be summarized in the following formula:
(28)Pn(J+1,2J)(t+Δt)=(1-∑k=1J‍αJ+kΔtkkkkkl-∑k=1J∑n1,…,nJ=0∞‍a(nJ+k,nk)μJ+k(1-pJ+k,J+k)kkkkkkkkkkl×Pn(1,J)(t)Δt∑k=1J‍)×Pn(J+1,2J)(t)+∑k=1J‍(∑n1,…,nJ=0∞‍a(nJ+k+1,nk)μJ+kkkkkkkkkkkkkkkkkll×(1-∑l=1J‍pJ+k,J+l)Pn(1,J)(t)Δt)kkkkkkk×PnJ+1,…,nJ+k+1,…,n2J(t)+∑k=1J‍(nJ+k∧1)αJ+kΔtPnJ+1,…,nJ+k-1,…,n2J(t)+∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)kkkkkkkkkkkkkkkkkk×μJ+rpJ+r,J+sPn(1,J)(t)∑n1,…,nJ=0∞‍)kkkkkkkkk×(nJ+s∧1)ΔtPnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J(t)+o(Δt),
where a(nJ+k,nk)=(nJ+k∧(cs(k)-nk))∨0 as defined in (10). o(Δt) is infinitesimal in terms of Δt; that is, o(Δt)/Δt→0 as Δt→∞. Equation (28) can be illustrated in the following disjoint events.Event
A. The J-dimensional queue length process Q(J+1,2J)(t) keeps at state n(J+1,2J) unchanging at times t and t+Δt; that is,
(29)Q(J+1,2J)(t)=n(J+1,2J),Q(J+1,2J)(t+Δt)=n(J+1,2J).
This event involves the following two parts.Part One. Suppose ΔSJ+k(t)=ΔAJ+k(t)=0 for all k∈{1,2,…,J}; that is, no external jobs arrive to classes with indices belonging to {J+1,…,2J} during [t,t+Δt] while no jobs finish their services, either because the jobs being served at time t require longer service times than Δt or because the services are blocked or interrupted by higher rank class jobs. Then, the probability for Part One can be expressed in the following product form:
(30)P{∑k=1J‍ΔSJ+k(t)=∑k=1J‍ΔAJ+k(t)=0}=P{∑k=1J‍ΔSJ+k(t)=∑k=1J‍ΔEJ+k(t)=0}=∏k=1J‍P{=∑s=kJ‍ΔEJ+s+1(t)=0ΔSJ+k(t)=ΔEJ+k(t)=0∣kkkkkkkkkk∑s=kJΔSJ+s+1(t)=∑s=kJ‍ΔEJ+s+1(t)=00∑s=kJΔSJ+s+1(t)}=∏k=1J‍P{EJ+k(t)}=∏k=1J‍(1-αJ+kΔt-aJ+k(t)μJ+kΔt+o(Δt))=1-∑k=1J‍αJ+kΔt-∑k=1J‍aJ+k(t)μJ+kΔt+o(Δt),
where aJ+k(t) is defined to be
(31)aJ+k(t)≡(QJ+k(t)∧(cs(k)-Qk(t)))∨0.
The event EJ+k(t) in the third equality of (30) is defined as follows:
(32)EJ+k(t)={∑s=kJ‍ΔSJ+k(t)=ΔEJ+k(t)=0∣kkkk∑s=kJΔSJ+s+1(t)=∑s=kJ‍ΔEJ+s+1(t)=00∑s=kJΔSJ+s+1(t)}.
To explain the fourth equality of (30), we introduce more notations. Let bk(t)Δt denote the probability that ΔAk(t)=1 for k∈{1,…,J}, and let cJ+k(t)Δt be the probability that ∑s=1k-1‍ΔSJ+s(t)=1 for k∈{2,3,…,J}; that is,
(33)bk(t)Δt≡P{ΔAk(t)=1},cJ+k(t)Δt≡P{∑s=1k-1‍ΔSJ+s(t)=1}.
These probabilities can be explicitly expressed in terms of the external arrival rates, service rates, and network states; for example,
(34)bk(t)Δt=(αk+∑l∈{1,…,J}‍plkμl(cs(l)∧Ql(t)))Δt.
Thus, by the independent assumptions on external arrival and service processes among different stations and classes, and for each k∈{2,3,…,J}, we have
(35)P{EJ+k(t)}=∑n=0∞∑m=0∞‍P{EJ+k(t)∣ΔAk(t)=n,∑s=1k-1‍ΔSJ+s(t)=m}kkkkkkkkkkkk×P{ΔAk(t)=n}P{∑s=1k-1‍ΔSJ+s(t)=n}=P{EJ+k(t)∣ΔAk(t)=0,∑s=1k-1‍ΔSJ+s(t)=0}×P{ΔAk(t)=0}P{∑s=1k-1‍ΔSJ+s(t)=0}+P{EJ+k(t)∣ΔAk(t)=1,∑s=1k-1‍ΔSJ+s(t)=0}×P{ΔAk(t)=1}P{∑s=1k-1‍ΔSJ+s(t)=0}+P{EJ+k(t)∣ΔAk(t)=0,∑s=1k-1‍ΔSJ+s(t)=1}×P{ΔAk(t)=0}P{∑s=0k-1‍ΔSJ+s(t)=1}+o(Δt)=(1-aJ+k(t)μJ+kΔt)(1-αJ+kΔt)×(1-bk(t)Δt)(1-cJ+k(t)Δt)+(1-((QJ+k(t)∧(cs(k)-(Qk(t)+1)))∨0)μJ+kΔt)×(1-αJ+kΔt)×(bk(t)Δt-pkkμk(cs(k)∧Qk(t))Δt)×(1-cJ+k(t)Δt)+(1-((QJ+k(t)∧(cs(k)-Qk(t)))∨0)μJ+kΔt)×(1-αJ+kΔt)×pkkμk(cs(k)∧Qk(t))Δt(1-cJ+k(t)Δt)+(1-(((QJ+k(t)+1)∧(cs(k)-Qk(t)))∨0)μJ+kΔt)×(1-αJ+kΔt)(1-bk(t)Δt)cJ+k(t)Δt+o(Δt)=1-αJ+kΔt-aJ+k(t)μJ+kΔt+o(Δt).
For the case k=1, one can use the similar way to check that the result in the above equation is also true.Part Two. Suppose ΔSJ+k(t)=ΔAJ+k(t)≥1 for at least one k∈{1,2,…,j}; that is, the number (at least one) of jobs that finished their services for class J+k in [t,t+Δt] equals that of jobs that arrived at the class. It is easy to see that the probability corresponding part one is given by
(36)P{ΔSJ+k(t)=ΔAJ+k(t)≥1klforatleastonek∈{1,2,…,j}ΔSJ+k(t)}=∑k=1J‍aJ+k(t)μJ+kpJ+k,J+kΔt+o(Δt),
where aJ+k(t) is defined in (31).Then, it follows from (30) and (36) that the probability for Event A is given by
(37)P{Q(J+1,2J)(t)=n(J+1,2J);Q(J+1,2J)(t+Δt)=n(J+1,2J)}=∑n1,…,nJ=0∞‍P{Q(J+1,2J)(t+Δt)=n(J+1,2J)∣kkkkkkkkkkklQ(1,2J)(t)=n(1,2J)Q(1,2J)(t)}kkkkkkkkk×P{Q(1,J)(t)=n(1,J)}kkkkkkkkk×P{Q(J+1,2J)(t)=n(J+1,2J)}+o(Δt)=(1-∑k=1J‍αJ+kΔtkkkk-∑k=1J‍∑n1,…,nJ=0∞‍a(nJ+k,nk)μJ+k(1-pJ+k,J+k)Pn(1,J)(t)Δt)×Pn(J+1,2J)(t)+o(Δt).Event
B. There is a k∈{1,…,J} such that QJ+k(t)=nJ+1+1, QJ+k(t+Δt)=nJ+1 and for all l∈{1,…,J}∖{k}, QJ+l(t)=QJ+l(t+Δt)=nJ+l. Similar to the discussion in Event A, we can obtain the probability for Event B as follows:
(38)∑k=1J‍P{QJ+k(t)=nJ+k+1,QJ+k(t+Δt)=nJ+k,kkkklkQJ+l(t)=QJ+l(t+Δt)=nJ+l,l∈{1,…,J}∖{k}}=∑k=1J‍P{∑l=1,l≠kJ‍ΔSJ+k(t)=1,ΔAJ+k(t)=0,kkkkklkk∑l=1,l≠kJ‍ΔSJ+l(t)=∑l=1,l≠kJ‍ΔAJ+l(t)=0,kkkkkklkkkkQJ+k(t)=nJ+k+1,QJ+l(t)=nJ+lkkkkkklkkkk∀l∈{1,…,J}∖{k}∑l=1,l≠kJ‍}+o(Δt)=∑k=1J‍P{∑l=1,l≠kJ‍ΔSJ+l(t)=ΔAJ+l(t)=0,ΔAJ+k(t)=0,kkkkkkkkkkkQJ+k(t)=nJ+k+1,QJ+l(t)=nJ+l,kkkkkkkkkkkl∈{1,…,J}∖{k}∣ΔSJ+k(t)=1∑l=1,l≠kJ‍}×P{ΔSJ+k(t)=1}+o(Δt)=∑k=1J‍(∑n1,…,nJ=0∞‍a(nJ+k+1,nk)μJ+kkkkkkkkkkkkkkkk×(1-∑l=1J‍pJ+k,J+l)Pn(1,J)(t)Δt)kkkkkkk×PnJ+1,…,nJ+k+1,…,n2J(t)+o(Δt).Event
C. There is a k∈{1,…,J} such that QJ+k(t)=nJ+k-1, QJ+k(t+Δt)=nJ+k and for all l∈{1,…,J}∖{k}, QJ+l(t)=QJ+l(t+Δt)=nJ+l. Then, we can get the probability for Event C as follows:
(39)∑k=1J‍P{QJ+k(t)=nJ+k-1,QJ+k(t+Δt)=nJ+k,kkkkkQJ+l(t)=QJ+l(t+Δt)=nJ+l,l∈{1,…,J}∖{k}}=∑k=1J‍P{∑l=1,l≠kJ‍ΔSJ+k(t)=0,ΔEJ+k(t)=1,kkkkkkk∑l=1,l≠kJ‍ΔSJ+l(t)=∑l=1,l≠kJ‍ΔAJ+l(t)=0,kkkkkkkkkQJ+k(t)=nJ+k-1,kkkkkkkkkQJ+l(t)=nJ+l∀l∈{1,…,J}∖{k}∑l=1,l≠kJ‍}+o(Δt)=∑k=1J‍(nJ+k∧1)αJ+kΔtPnJ+1,…,nJ+k-1,…,n2J(t)+o(Δt).Event
D. There exist r,s∈{1,…,J} such that QJ+r(t)=nJ+r+1, QJ+r(t+Δt)=nJ+s, QJ+s(t)=nJ+s-1, and QJ+s(t+Δt)=nJ+s and for all l∈{1,…,J}∖{r,s}, QJ+l(t)=QJ+l(t+Δt)=nJ+l. Then, we can get the probability for Event D as follows:
(40)∑r,s=1,s≠rJ‍P{QJ+r(t)=nJ+r+1,QJ+r(t+Δt)=nJ+r;kkkkkkkkkQJ+s(t)=nJ+s-1,QJ+s(t+Δt)=nJ+s;kkkkkkkkkQJ+l(t)=QJ+l(t+Δt)=nJ+lkkkkkkkkk∀l∈{1,…,J}∖{r,s}QJ+r(t)=nJ+r}=∑r,s=1,s≠rJ‍P{∑l=1,l≠r,sJ‍ΔSJ+r(t)=1,ΔAJ+r(t)=0;ΔSJ+s(t)=0,kkkkkkkkkkkkkΔAJ+s(t)=1;kkkkkkkkkkkkk∑l=1,l≠r,sJ‍ΔSJ+l(t)=∑l=1,l≠r,sJ‍ΔAJ+l(t)=0;kkkkkkkkkkkkkQJ+r(t)=nJ+r+1,QJ+s(t)=nJ+s-1,kkkkkkkkkkkkkQJ+l(t)=nJ+l∀l∈{1,…,J}∖{r,s}∑l=1,l≠r,sJ‍}+o(Δt)=∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)μJ+rpJ+r,J+sPn(1,J)(t))kkkkkkkkkk×(nJ+s∧1)PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J(t)+o(Δt).Now we go back to discuss (28). After transferring the Pn(J+1,2J)(t) from right to left, dividing by Δt, and taking the limit as Δt tends to zero, we get the following differential equations:
(41)P˙n(J+1,2J)(t)=(-∑k=1J‍αJ+kkkkk-∑k=1J‍∑n1,…,nJ=0∞‍a(nJ+k,nk)μJ+k(1-pJ+k,J+k)Pn(1,J)(t))×Pn(J+1,2J)(t)+∑k=1J‍(∑n1,…,nJ=0∞‍a(nJ+k+1,nk)μJ+kkkkkkkkkkkkkkkk×(1-∑l=1J‍pJ+k,J+l)Pn(1,J)(t))kkkkk×PnJ+1,…,nJ+k+1,…,n2J(t)+∑k=1J‍(nJ+k∧1)αJ+kPnJ+1,…,nJ+k-1,…,n2J(t)+∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)μJ+rpJ+r,J+sPn(1,J)(t))kkkkkkkkkkk×(nJ+s∧1)PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J(t).
Next we show that the given distribution corresponding h=2 in Theorem 1 is a steady-state solution of those equations described by (41). It is enough to demonstrate that the derivatives in the above equations are all made zero by setting Pn(1,2J)(t)=Pn(1,2J), that is, to prove that
(42)(∑k=1J‍αJ+k+∑k=1J‍∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)μJ+k)Pn(J+1,2J)=∑k=1J‍((1-∑l=1J‍pJ+k,J+l)∑n1,…,nJ=0∞‍a(nJ+k+1,nk)Pn(1,J)μJ+kkkkkkkkkkkkkk×(1-∑l=1J‍pJ+k,J+l))kkkkkk×PnJ+1,…,nJ+k+1,…,n2J+∑k=1J‍(nJ+k∧1)αJ+kPnJ+1,…,nJ+k-1,…,n2J+∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)Pn(1,J)μJ+rpJ+r,J+s)kkkkkkkkklk×(nJ+s∧1)PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J+∑k=1J∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)μJ+kpJ+k,J+kPn(J+1,2J).
By the definition of κ(·) in Theorem 1, that is, for each k∈{1,…,J},
(43)κ(nJ+k)=∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J),
then, we can rewrite (42) as follows:
(44)(∑k=1J‍αJ+k+∑k=1J‍κ(nJ+k)μJ+k)Pn(J+1,2J)=∑k=1J‍κ(nJ+k+1)μJ+k(1-∑l=1J‍pJ+k,J+l)PnJ+1,…,nJ+k+1,…,n2J+∑k=1J‍(nJ+k∧1)αJ+kPnJ+1,…,nJ+k-1,…,n2J+∑r,s=1,s≠rJ‍κ(nJ+r+1)μJ+rpJ+r,J+s(nJ+s∧1)kkkkkkkkkk×PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J+∑k=1J‍κ(nJ+k)μJ+kpJ+k,J+kPn(J+1,2J).
From the given distribution in Theorem 1, we have
(45)PnJ+1,…,nJ+k+1,…,n2JPn(J+1,2J)=λJ+kκ(nJ+k+1)μJ+k,PnJ+1,…,nJ+k-1,…,n2JPn(J+1,2J)=κ(nJ+k)μJ+kλJ+k,PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2JPn(J+1,2J)=κ(nJ+s)λJ+rμJ+sκ(nJ+r+1)λJ+sμJ+r.
Then, by (44) to (45), we know that it is sufficient to show
(46)∑k=1J‍αJ+k+∑k=1J‍κ(nJ+k)μJ+k=∑k=1J‍(1-∑l=1J‍pJ+k,J+l)λJ+k+∑k=1J‍κ(nJ+k)αJ+kμJ+kλJ+k+∑r,s=1J‍κ(nJ+s)λJ+rμJ+spJ+r,J+sλJ+s.
Owing to the routing structure and the traffic equation (3), we have the following observations:
(47)∑k=1J‍(1-∑l=1J‍pJ+k,J+l)λJ+k=∑k=1J‍αJ+k,∑r,s=1J‍κ(nJ+s)λJ+rμJ+spJ+r,J+sλJ+s=∑s=1J‍(κ(nJ+s)μJ+sλJ+s)∑r=1J‍λJ+rpJ+r,J+s=∑s=1J‍κ(nJ+s)μJ+s-∑s=1J‍κ(nJ+s)μJ+sαJ+sλJ+s.
Then, substituting (47) into (46), we get the necessary equality.Next, we derive the initial distributionP0s(J+k) corresponding to state nJ+k=0. By the network stability conditions (12) and ρj<1 for each j∈{1,…,J}, it follows from ∑nJ+k=0∞‍PnJ+ks(J+k)=1 that the following initial distribution is well posed:
(48)P0s(J+k)=11+∑nJ+k=1∞∏rJ+k=1nJ+k(λJ+k/κ(rJ+k)μJ+k).
Hence, we complete the proof of Theorem 1 for priority types h∈{1,2}.Step 3.
To finish the induction procedure, in this step, we first suppose that the result described in Theorem1 is true for all classes with priority rank h∈{1,…,I-1}; then, we show that it is true for all classes with priority rank h∈{1,…,I}. By the similar illustration in getting (41), we have the following differential equations:
(49)P˙n((I-1)J+1,IJ)(t)=-∑k=1J‍α(I-1)J+kPn((I-1)J,IJ)(t)-∑k=1J∑n1,…,n(I-1)J=0∞‍a(nIJ+k,∑u=0I-2‍nuJ+k)μ(I-1)J+kkkkkkklk×(1-p(I-1)J+k,(I-1)J+k)kkkkkklk×Pn(1,(I-1)J)(t)Pn((I-1)J+1,IJ)(t)+∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k+1,∑u=0I-2‍nuJ+k)μ(I-1)J+kkkkkkkkk×(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)Pn(1,(I-1)J)(t)kkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+k+1,…,nIJ(t)+∑k=1J‍(n(I-1)J+k∧1)α(I-1)J+kkkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+k-1,…,nIJ(t)+∑r,s=1,s≠rJ∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+r+1,∑u=0I-2‍nuJ+r)kkkkkkkkkkk×μ(I-1)J+rp(I-1)J+r,(I-1)J+skkkkkkkkkkk×Pn(1,(I-1)J)(t)(n(I-1)J+s∧1)kkkkkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,nIJ(t).
Next, we show that the given distribution corresponding the lowest priority type in Theorem 1 is a steady-state solution of those equations described by (49). It suffices to demonstrate that the derivatives in the above equations are all made zero by setting Pn((I-1)J,IJ)(t)=Pn((I-1)J,IJ), that is, to prove that
(50)∑k=1J‍α(I-1)J+kPn((I-1)J,IJ)+∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k,∑u=0I-2‍nuJ+k)kkkkkk×Pn(1,(I-1)J)μ(I-1)J+kPn((I-1)J+1,IJ)=∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k+1,∑u=0I-2‍nuJ+k)Pn(1,(I-1)J)kklkkkk×μ(I-1)J+k(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)kklkkkk×Pn(I-1)J+1,…,n(I-1)J+k+1,…,nIJ+∑k=1J‍(n(I-1)J+k∧1)α(I-1)J+kPn(I-1)J+1,…,n(I-1)J+k-1,…,nIJ+∑r,s=1,s≠rJ∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+r+1,∑u=0I-2‍nuJ+r)kkkkkkkkkkk×Pn(1,(I-1)J)μ(I-1)J+r×p(I-1)J+r,(I-1)J+s(n(I-1)J+s∧1)×Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,nIJ+∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k,∑u=0I-2‍nuJ+k)Pn(1,(I-1)J)kkklkkkkl×μ(I-1)J+kp(I-1)J+k,(I-1)J+kPn((I-1)J+1,IJ).
By the definition of κ(·) in Theorem 1, we have
(51)κ(n(I-1)J+k)=∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k,∑u=0I-2‍nuJ+k)Pn(1,(I-1)J).
Then, we can rewrite (50) as follows:
(52)(∑k=1J‍α(I-1)J+k+∑k=1J‍κ(n(I-1)J+k)μ(I-1)J+k)Pn((I-1)J+1,IJ)=∑k=1J‍κ(n(I-1)J+k+1)μ(I-1)J+k(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)kkkkkkk×Pn(I-1)J+1,…,n(I-1)J+k+1,…,nIJ+∑k=1J‍(n(I-1)J+k∧1)α(I-1)J+kPn(I-1)J+1,…,n(I-1)J+k-1,…,nIJ+∑r,s=1,s≠rJ‍κ(n(I-1)J+k+1)kkkkkkkkkkkk×μ(I-1)J+rp(I-1)J+r,(I-1)J+s(n(I-1)J+s∧1)kkkkkkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,nIJ+∑k=1J‍κ(n(I-1)J+k)μ(I-1)J+kp(I-1)J+k,(I-1)J+kPn((I-1)J+1,IJ).
From the given distribution in Theorem 1, we have
(53)Pn(I-1)J+1,…,n(I-1)J+k+1,…,n2JPn((I-1)J+1,2J)=λ(I-1)J+kκ(n(I-1)J+k+1)μ(I-1)J+k,Pn(I-1)J+1,…,n(I-1)J+k-1,…,n2JPn((I-1)J+1,2J)=κ(n(I-1)J+k)μ(I-1)J+kλ(I-1)J+k,Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,n2JPn((I-1)J+1,2J)=κ(n(I-1)J+s)λ(I-1)J+rμ(I-1)J+sκ(n(I-1)J+r+1)λ(I-1)J+sμ(I-1)J+r.
From (52) to (53), it is sufficient to show
(54)∑k=1J‍α(I-1)J+k+∑k=1J‍κ(n(I-1)J+k)μ(I-1)J+k=∑k=1J‍(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)λ(I-1)J+k+∑k=1J‍κ(n(I-1)J+k)α(I-1)J+kμ(I-1)J+kλ(I-1)J+k+∑r,s=1J‍κ(n(I-1)J+s)λ(I-1)J+rμ(I-1)J+sp(I-1)J+r,(I-1)J+sλ(I-1)J+s.
Owing to the routing structure and the traffic equation (3), we have the following observations:
(55)∑k=1J‍(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)λ(I-1)J+k=∑k=1J‍α(I-1)J+k,∑r,s=1J‍κ(n(I-1)J+s)λ(I-1)J+rμ(I-1)J+sp(I-1)J+r,(I-1)J+sλ(I-1)J+s=∑s=1J‍(κ(n(I-1)J+s)μ(I-1)J+sλ(I-1)J+s)∑r=1J‍λ(I-1)J+rp(I-1)J+r,(I-1)J+s=∑s=1J‍κ(n(I-1)J+s)μ(I-1)J+s-∑s=1J‍κ(n(I-1)J+s)μ(I-1)J+sα(I-1)J+sλ(I-1)J+s.
Then, substituting (55) into (54), we get the necessary equality.Next, we derive the initial distributionP0s((I-1)J+k) corresponding to state n(I-1)J+k=0. By the network stability conditions (12) and ρj<1 for each j∈{1,…,J}, it follows from ∑n(I-1)J+k=0∞‍Pn(I-1)J+ks((I-1)J+k)=1 that the following initial distribution is well posed:
(56)P0s((I-1)J+k)=(1+∑n(I-1)J+k=1∞‍∏r(I-1)J+k=1n(I-1)J+k‍λ(I-1)J+kκ(r(I-1)J+k)μ(I-1)J+k)-1.
Hence, we complete the proof of Theorem 1.


### 5.2. Proof of Proposition2

Net I. We justify the stability condition (12) for Net I by induction in terms of the number of job types; that is, h=1,2,…,I.First, we consider the case thath=2. Note that cj=1 for all j∈{1,2,…,J} and step one in the proof of Theorem 1, we have that, for nJ+k>0 with k∈{1,2,…,J},
(57)κ(nJ+k)=P0s(k)=1-λkμk.
Therefore, condition (12) is true for h=2. Furthermore, we have
(58)P0s(J+k)=1-λJ+k(1-λk/μk)μJ+k.Second, forh≤I-1, we suppose that
(59)κ(n(h-1)J+k)=1-∑i=0h-2‍λiJ+kμiJ+k.
Hence, we have
(60)P0s((h-1)J+k)=1-λ(h-1)J+k(1-∑i=0h-2λiJ+k/μiJ+k)μ(h-1)J+k.
From the induction assumptions (59) and (60), we know that (12) is true for h≤I-1.Finally, we show that (12) holds for h=I. In fact, from the definition of κ(n(I-1)J+k), we know that, for n(I-1)J+k>0,
(61)κ(n(I-1)J+k)=∏i=0I-2‍P0s(iJ+k)=1-∑i=0I-2‍λiJ+kμiJ+k.
Then, we see that (12) is true for h=I. Hence, we complete the proof of Net I.Net II. Note that, for each k∈{1,…,J} and nJ+k≥cs(k), we have
(62)κ(nJ+k)=∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)=∑nk=0cs(k)-1‍a(nJ+k,nk)Pnks(k)=P0s(k)(cs(k)+(cs(k)-1)(λkμk)+⋯+(λk/μk)cs(k)-1(cs(k)-1)!)=(∑nk=0cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!)×(∑nk=0cs(k)-1‍(λk/μk)nknk!+(λk/μk)cs(k)cs(k)!(1-λk/cs(k)μk))-1=(cs(k)(1-λkcs(k)μk)∑nk=0cs(k)-1‍(cs(k)-nk)(λk/μk)nkcs(k)nk!)×(1+∑nk=0cs(k)-2‍(1(nk+1)!-1cs(k)nk!)(λkμk)nk+1)-1=cs(k)(1-λkcs(k)μk).
Furthermore, for 0<nJ+k<cs(k), we have
(63)κ(nJ+k)=∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)=∑nk=0cs(k)-1‍a(nJ+k,nk)Pnks(k)=P0s(k)((λkμk)cs(k)-1nJ+k(1+11!(λkμk)kkkkkkkkkkkkkkk+⋯+1(cs(k)-nJ+k)!(λkμk)cs(k)-nJ+k)kkkkkkkkl+nJ+k-1(cs(k)-nJ+k+1)!(λkμk)cs(k)-nJ+k+1kkkkkkkkl+⋯+1(cs(k)-1)!(λkμk)cs(k)-1)=(nJ+k∑nk=0cs(k)-nJ+k‍(λk/μk)nknk!kkkkkl+∑nk=cs(k)-nJ+k+1cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!)×(∑nk=0cs(k)-1‍(λk/μk)nknk!+(λk/μk)cs(k)cs(k)!(1-λk/cs(k)μk))-1=(+∑nk=cs(k)-nJ+k+1cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!nJ+k)nJ+k(1-λkcs(k)μk)kkkkk×(∑nk=0cs(k)-nJ+k‍(λk/μk)nknk!kkkkkkkkkkk+∑nk=cs(k)-nJ+k+1cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!nJ+k))×(1+∑nk=0cs(k)-2‍(1(nk+1)!-1cs(k)nk!)(λkμk)nk+1)-1≥nJ+k(1-λkcs(k)μk).
Then, we can see that (12) is true. Hence, we complete the proof of Proposition 2.


### 5.3. Illustration of Example3

For each classk∈{1,2,…,J}, the initial distribution P0s(k) can be calculated as follows, since ρs(k)<1:
(64)P0s(k)=1-λk/2μk1+λk/2μk.
Then, for nJ+k=1, we have
(65)κ(nJ+k)=P0s(k)(1+λkμk)=(1-λk/2μk)(1+λk/μk)1+λk/2μk≥1-λk2μk.
Furthermore, for nJ+k≥2, we have
(66)κ(nJ+k)=P0s(k)(2+λkμk)=2(1-λk2μk).
Therefore, since ρs(k)<1, we have
(67)P0s(J+k)=(1+λJ+kκ(1)μJ+kk+λJ+kκ(1)μJ+kλJ+k/μJ+k(1-λJ+k/κ(2)μJ+k)κ(2)μJ+k)-1=(1-λk/2μk-λJ+k/2μJ+k)(1+λk/μk)1+λk/2μk+λJ+k/2μJ+k-(1/2)(λk/μk)2.
Finally, we can compute κ(n2J+k) for n2J+k≥2 as follows:
(68)κ(n2J+k)=2P0s(k)P0s(J+k)+P1s(k)P0s(J+k)+P0s(k)P1s(J+k)=P0s(k)P0s(J+k)(2+λkμk+λJ+kμJ+k)=(2(1-λk2μk-λJ+k2μJ+k)(1+λkμk)(1-λk2μk)kkkkk×(1+λk2μk+λJ+k2μJ+k))×((1+λk2μk+λJ+k2μJ+k-12(λkμk)2)(1+λk2μk))-1=(2(1-λk2μk-λJ+k2μJ+k)kkkkk×(1+λkμk+λJ+k2μJ+k-(λk2μk)2-14(λkμk)3kkkkkkkk+2(1-λk2μk)λkλJ+k4μkμJ+k))×(1+λkμk+λJ+k2μJ+k-(λk2μk)2kkkkkkk-14(λkμk)3+λkλJ+k4μkμJ+k)-1.
Thus, if λk/μk≤1, from the condition ρs(k)<1, we see that
(69)λ2J+kκ(n2J+k)μ2J+k≤λ2J+k2(1-λk/2μk-λJ+k/2μJ+k)μ2J+k<1;
that is, the condition (12) holds.


## 5.1. Proof of Theorem1

For convenience, we introduce some additional notations. LetΔEk(t), ΔSk(t), and ΔAk(t) be defined by
(22)ΔEk(t)≡Ek(t+Δt)-Ek(t),ΔSk(t)≡Sk(t+Δt)-Sk(t),ΔAk(t)≡Ak(t+Δt)-Ak(t),
which denote the cumulative external arrivals to class k, the cumulative number of jobs finished services at class k, and the total arrivals to class k in [t,t+Δt]. Then, we can justify Theorem 1 by induction as in the following several steps.Step 1.
We consider the steady-state distribution for the highest rank classes with each indexk∈{1,2,…,J}. In this case, the type index i=0 is used in Theorem 1. Owing to the preemptive-resume service discipline and the class routing structure, we know that these J classes form a Jackson network. Then, by the theorem in Jackson [24], we have the following product form:
(23)Pn(1,J)=Pn1s(1)Pn2s(2)⋯PnJs(J),
where Pnks(k) denotes the steady-state probability at state nk for class k(∈{1,2,…,J}) at station s(k). More precisely, it is given by
(24)Pnks(k)={P0s(k)(λk/μk)nknk!,if0≤nk<ck,P0s(k)(λk/μk)nknk!(nk)nk-ck,ifnk≥ck
with P0s(k) being determined by equation ∑nk=0∞‍Pnks(k)=1; that is,
(25)P0s(k)=(∑nk=0cs(k)-1‍1nk!(λkμk)nk+(λk/μk)cs(k)cs(k)!(1-λk/cs(k)μk))-1.Step 2.
We derive the steady-state distribution for the highest rank and the second highest rank classes withk∈{1,2,…,2J}. Owing to the preemptive-resume service discipline and the job routing structure, we have
(26)Pn(1,2J)(t)=Pn(1,J)/n(J+1,2J)(t)Pn(J+1,2J)(t)=Pn(1,J)(t)Pn(J+1,2J)(t),
where the second equality follows from the preemptive assumption and Pn(1,J)/n(J+1,2J)(t) is the conditional probability at time t for classes with k∈{1,…,J} at state n(1,J) in terms of classes with k∈{J+1,…,2J} at state n(J+1,2J); that is,
(27)Pn(1,J)/n(J+1,2J)(t)=P{Q(1,J)(t)=n(1,J)∣Q(J+1,2J)(t)=n(J+1,2J)}.
In order to get the steady state distribution for Pn(J+1,2J)(t), we consider each state n(J+1,2J) at time t for the second highest rank class jobs. There are several ways in which the system can reach it. They can be summarized in the following formula:
(28)Pn(J+1,2J)(t+Δt)=(1-∑k=1J‍αJ+kΔtkkkkkl-∑k=1J∑n1,…,nJ=0∞‍a(nJ+k,nk)μJ+k(1-pJ+k,J+k)kkkkkkkkkkl×Pn(1,J)(t)Δt∑k=1J‍)×Pn(J+1,2J)(t)+∑k=1J‍(∑n1,…,nJ=0∞‍a(nJ+k+1,nk)μJ+kkkkkkkkkkkkkkkkkll×(1-∑l=1J‍pJ+k,J+l)Pn(1,J)(t)Δt)kkkkkkk×PnJ+1,…,nJ+k+1,…,n2J(t)+∑k=1J‍(nJ+k∧1)αJ+kΔtPnJ+1,…,nJ+k-1,…,n2J(t)+∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)kkkkkkkkkkkkkkkkkk×μJ+rpJ+r,J+sPn(1,J)(t)∑n1,…,nJ=0∞‍)kkkkkkkkk×(nJ+s∧1)ΔtPnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J(t)+o(Δt),
where a(nJ+k,nk)=(nJ+k∧(cs(k)-nk))∨0 as defined in (10). o(Δt) is infinitesimal in terms of Δt; that is, o(Δt)/Δt→0 as Δt→∞. Equation (28) can be illustrated in the following disjoint events.Event
A. The J-dimensional queue length process Q(J+1,2J)(t) keeps at state n(J+1,2J) unchanging at times t and t+Δt; that is,
(29)Q(J+1,2J)(t)=n(J+1,2J),Q(J+1,2J)(t+Δt)=n(J+1,2J).
This event involves the following two parts.Part One. Suppose ΔSJ+k(t)=ΔAJ+k(t)=0 for all k∈{1,2,…,J}; that is, no external jobs arrive to classes with indices belonging to {J+1,…,2J} during [t,t+Δt] while no jobs finish their services, either because the jobs being served at time t require longer service times than Δt or because the services are blocked or interrupted by higher rank class jobs. Then, the probability for Part One can be expressed in the following product form:
(30)P{∑k=1J‍ΔSJ+k(t)=∑k=1J‍ΔAJ+k(t)=0}=P{∑k=1J‍ΔSJ+k(t)=∑k=1J‍ΔEJ+k(t)=0}=∏k=1J‍P{=∑s=kJ‍ΔEJ+s+1(t)=0ΔSJ+k(t)=ΔEJ+k(t)=0∣kkkkkkkkkk∑s=kJΔSJ+s+1(t)=∑s=kJ‍ΔEJ+s+1(t)=00∑s=kJΔSJ+s+1(t)}=∏k=1J‍P{EJ+k(t)}=∏k=1J‍(1-αJ+kΔt-aJ+k(t)μJ+kΔt+o(Δt))=1-∑k=1J‍αJ+kΔt-∑k=1J‍aJ+k(t)μJ+kΔt+o(Δt),
where aJ+k(t) is defined to be
(31)aJ+k(t)≡(QJ+k(t)∧(cs(k)-Qk(t)))∨0.
The event EJ+k(t) in the third equality of (30) is defined as follows:
(32)EJ+k(t)={∑s=kJ‍ΔSJ+k(t)=ΔEJ+k(t)=0∣kkkk∑s=kJΔSJ+s+1(t)=∑s=kJ‍ΔEJ+s+1(t)=00∑s=kJΔSJ+s+1(t)}.
To explain the fourth equality of (30), we introduce more notations. Let bk(t)Δt denote the probability that ΔAk(t)=1 for k∈{1,…,J}, and let cJ+k(t)Δt be the probability that ∑s=1k-1‍ΔSJ+s(t)=1 for k∈{2,3,…,J}; that is,
(33)bk(t)Δt≡P{ΔAk(t)=1},cJ+k(t)Δt≡P{∑s=1k-1‍ΔSJ+s(t)=1}.
These probabilities can be explicitly expressed in terms of the external arrival rates, service rates, and network states; for example,
(34)bk(t)Δt=(αk+∑l∈{1,…,J}‍plkμl(cs(l)∧Ql(t)))Δt.
Thus, by the independent assumptions on external arrival and service processes among different stations and classes, and for each k∈{2,3,…,J}, we have
(35)P{EJ+k(t)}=∑n=0∞∑m=0∞‍P{EJ+k(t)∣ΔAk(t)=n,∑s=1k-1‍ΔSJ+s(t)=m}kkkkkkkkkkkk×P{ΔAk(t)=n}P{∑s=1k-1‍ΔSJ+s(t)=n}=P{EJ+k(t)∣ΔAk(t)=0,∑s=1k-1‍ΔSJ+s(t)=0}×P{ΔAk(t)=0}P{∑s=1k-1‍ΔSJ+s(t)=0}+P{EJ+k(t)∣ΔAk(t)=1,∑s=1k-1‍ΔSJ+s(t)=0}×P{ΔAk(t)=1}P{∑s=1k-1‍ΔSJ+s(t)=0}+P{EJ+k(t)∣ΔAk(t)=0,∑s=1k-1‍ΔSJ+s(t)=1}×P{ΔAk(t)=0}P{∑s=0k-1‍ΔSJ+s(t)=1}+o(Δt)=(1-aJ+k(t)μJ+kΔt)(1-αJ+kΔt)×(1-bk(t)Δt)(1-cJ+k(t)Δt)+(1-((QJ+k(t)∧(cs(k)-(Qk(t)+1)))∨0)μJ+kΔt)×(1-αJ+kΔt)×(bk(t)Δt-pkkμk(cs(k)∧Qk(t))Δt)×(1-cJ+k(t)Δt)+(1-((QJ+k(t)∧(cs(k)-Qk(t)))∨0)μJ+kΔt)×(1-αJ+kΔt)×pkkμk(cs(k)∧Qk(t))Δt(1-cJ+k(t)Δt)+(1-(((QJ+k(t)+1)∧(cs(k)-Qk(t)))∨0)μJ+kΔt)×(1-αJ+kΔt)(1-bk(t)Δt)cJ+k(t)Δt+o(Δt)=1-αJ+kΔt-aJ+k(t)μJ+kΔt+o(Δt).
For the case k=1, one can use the similar way to check that the result in the above equation is also true.Part Two. Suppose ΔSJ+k(t)=ΔAJ+k(t)≥1 for at least one k∈{1,2,…,j}; that is, the number (at least one) of jobs that finished their services for class J+k in [t,t+Δt] equals that of jobs that arrived at the class. It is easy to see that the probability corresponding part one is given by
(36)P{ΔSJ+k(t)=ΔAJ+k(t)≥1klforatleastonek∈{1,2,…,j}ΔSJ+k(t)}=∑k=1J‍aJ+k(t)μJ+kpJ+k,J+kΔt+o(Δt),
where aJ+k(t) is defined in (31).Then, it follows from (30) and (36) that the probability for Event A is given by
(37)P{Q(J+1,2J)(t)=n(J+1,2J);Q(J+1,2J)(t+Δt)=n(J+1,2J)}=∑n1,…,nJ=0∞‍P{Q(J+1,2J)(t+Δt)=n(J+1,2J)∣kkkkkkkkkkklQ(1,2J)(t)=n(1,2J)Q(1,2J)(t)}kkkkkkkkk×P{Q(1,J)(t)=n(1,J)}kkkkkkkkk×P{Q(J+1,2J)(t)=n(J+1,2J)}+o(Δt)=(1-∑k=1J‍αJ+kΔtkkkk-∑k=1J‍∑n1,…,nJ=0∞‍a(nJ+k,nk)μJ+k(1-pJ+k,J+k)Pn(1,J)(t)Δt)×Pn(J+1,2J)(t)+o(Δt).Event
B. There is a k∈{1,…,J} such that QJ+k(t)=nJ+1+1, QJ+k(t+Δt)=nJ+1 and for all l∈{1,…,J}∖{k}, QJ+l(t)=QJ+l(t+Δt)=nJ+l. Similar to the discussion in Event A, we can obtain the probability for Event B as follows:
(38)∑k=1J‍P{QJ+k(t)=nJ+k+1,QJ+k(t+Δt)=nJ+k,kkkklkQJ+l(t)=QJ+l(t+Δt)=nJ+l,l∈{1,…,J}∖{k}}=∑k=1J‍P{∑l=1,l≠kJ‍ΔSJ+k(t)=1,ΔAJ+k(t)=0,kkkkklkk∑l=1,l≠kJ‍ΔSJ+l(t)=∑l=1,l≠kJ‍ΔAJ+l(t)=0,kkkkkklkkkkQJ+k(t)=nJ+k+1,QJ+l(t)=nJ+lkkkkkklkkkk∀l∈{1,…,J}∖{k}∑l=1,l≠kJ‍}+o(Δt)=∑k=1J‍P{∑l=1,l≠kJ‍ΔSJ+l(t)=ΔAJ+l(t)=0,ΔAJ+k(t)=0,kkkkkkkkkkkQJ+k(t)=nJ+k+1,QJ+l(t)=nJ+l,kkkkkkkkkkkl∈{1,…,J}∖{k}∣ΔSJ+k(t)=1∑l=1,l≠kJ‍}×P{ΔSJ+k(t)=1}+o(Δt)=∑k=1J‍(∑n1,…,nJ=0∞‍a(nJ+k+1,nk)μJ+kkkkkkkkkkkkkkkk×(1-∑l=1J‍pJ+k,J+l)Pn(1,J)(t)Δt)kkkkkkk×PnJ+1,…,nJ+k+1,…,n2J(t)+o(Δt).Event
C. There is a k∈{1,…,J} such that QJ+k(t)=nJ+k-1, QJ+k(t+Δt)=nJ+k and for all l∈{1,…,J}∖{k}, QJ+l(t)=QJ+l(t+Δt)=nJ+l. Then, we can get the probability for Event C as follows:
(39)∑k=1J‍P{QJ+k(t)=nJ+k-1,QJ+k(t+Δt)=nJ+k,kkkkkQJ+l(t)=QJ+l(t+Δt)=nJ+l,l∈{1,…,J}∖{k}}=∑k=1J‍P{∑l=1,l≠kJ‍ΔSJ+k(t)=0,ΔEJ+k(t)=1,kkkkkkk∑l=1,l≠kJ‍ΔSJ+l(t)=∑l=1,l≠kJ‍ΔAJ+l(t)=0,kkkkkkkkkQJ+k(t)=nJ+k-1,kkkkkkkkkQJ+l(t)=nJ+l∀l∈{1,…,J}∖{k}∑l=1,l≠kJ‍}+o(Δt)=∑k=1J‍(nJ+k∧1)αJ+kΔtPnJ+1,…,nJ+k-1,…,n2J(t)+o(Δt).Event
D. There exist r,s∈{1,…,J} such that QJ+r(t)=nJ+r+1, QJ+r(t+Δt)=nJ+s, QJ+s(t)=nJ+s-1, and QJ+s(t+Δt)=nJ+s and for all l∈{1,…,J}∖{r,s}, QJ+l(t)=QJ+l(t+Δt)=nJ+l. Then, we can get the probability for Event D as follows:
(40)∑r,s=1,s≠rJ‍P{QJ+r(t)=nJ+r+1,QJ+r(t+Δt)=nJ+r;kkkkkkkkkQJ+s(t)=nJ+s-1,QJ+s(t+Δt)=nJ+s;kkkkkkkkkQJ+l(t)=QJ+l(t+Δt)=nJ+lkkkkkkkkk∀l∈{1,…,J}∖{r,s}QJ+r(t)=nJ+r}=∑r,s=1,s≠rJ‍P{∑l=1,l≠r,sJ‍ΔSJ+r(t)=1,ΔAJ+r(t)=0;ΔSJ+s(t)=0,kkkkkkkkkkkkkΔAJ+s(t)=1;kkkkkkkkkkkkk∑l=1,l≠r,sJ‍ΔSJ+l(t)=∑l=1,l≠r,sJ‍ΔAJ+l(t)=0;kkkkkkkkkkkkkQJ+r(t)=nJ+r+1,QJ+s(t)=nJ+s-1,kkkkkkkkkkkkkQJ+l(t)=nJ+l∀l∈{1,…,J}∖{r,s}∑l=1,l≠r,sJ‍}+o(Δt)=∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)μJ+rpJ+r,J+sPn(1,J)(t))kkkkkkkkkk×(nJ+s∧1)PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J(t)+o(Δt).Now we go back to discuss (28). After transferring the Pn(J+1,2J)(t) from right to left, dividing by Δt, and taking the limit as Δt tends to zero, we get the following differential equations:
(41)P˙n(J+1,2J)(t)=(-∑k=1J‍αJ+kkkkk-∑k=1J‍∑n1,…,nJ=0∞‍a(nJ+k,nk)μJ+k(1-pJ+k,J+k)Pn(1,J)(t))×Pn(J+1,2J)(t)+∑k=1J‍(∑n1,…,nJ=0∞‍a(nJ+k+1,nk)μJ+kkkkkkkkkkkkkkkk×(1-∑l=1J‍pJ+k,J+l)Pn(1,J)(t))kkkkk×PnJ+1,…,nJ+k+1,…,n2J(t)+∑k=1J‍(nJ+k∧1)αJ+kPnJ+1,…,nJ+k-1,…,n2J(t)+∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)μJ+rpJ+r,J+sPn(1,J)(t))kkkkkkkkkkk×(nJ+s∧1)PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J(t).
Next we show that the given distribution corresponding h=2 in Theorem 1 is a steady-state solution of those equations described by (41). It is enough to demonstrate that the derivatives in the above equations are all made zero by setting Pn(1,2J)(t)=Pn(1,2J), that is, to prove that
(42)(∑k=1J‍αJ+k+∑k=1J‍∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)μJ+k)Pn(J+1,2J)=∑k=1J‍((1-∑l=1J‍pJ+k,J+l)∑n1,…,nJ=0∞‍a(nJ+k+1,nk)Pn(1,J)μJ+kkkkkkkkkkkkkk×(1-∑l=1J‍pJ+k,J+l))kkkkkk×PnJ+1,…,nJ+k+1,…,n2J+∑k=1J‍(nJ+k∧1)αJ+kPnJ+1,…,nJ+k-1,…,n2J+∑r,s=1,s≠rJ‍(∑n1,…,nJ=0∞‍a(nJ+r+1,nr)Pn(1,J)μJ+rpJ+r,J+s)kkkkkkkkklk×(nJ+s∧1)PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J+∑k=1J∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)μJ+kpJ+k,J+kPn(J+1,2J).
By the definition of κ(·) in Theorem 1, that is, for each k∈{1,…,J},
(43)κ(nJ+k)=∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J),
then, we can rewrite (42) as follows:
(44)(∑k=1J‍αJ+k+∑k=1J‍κ(nJ+k)μJ+k)Pn(J+1,2J)=∑k=1J‍κ(nJ+k+1)μJ+k(1-∑l=1J‍pJ+k,J+l)PnJ+1,…,nJ+k+1,…,n2J+∑k=1J‍(nJ+k∧1)αJ+kPnJ+1,…,nJ+k-1,…,n2J+∑r,s=1,s≠rJ‍κ(nJ+r+1)μJ+rpJ+r,J+s(nJ+s∧1)kkkkkkkkkk×PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2J+∑k=1J‍κ(nJ+k)μJ+kpJ+k,J+kPn(J+1,2J).
From the given distribution in Theorem 1, we have
(45)PnJ+1,…,nJ+k+1,…,n2JPn(J+1,2J)=λJ+kκ(nJ+k+1)μJ+k,PnJ+1,…,nJ+k-1,…,n2JPn(J+1,2J)=κ(nJ+k)μJ+kλJ+k,PnJ+1,…,nJ+r+1,…,nJ+s-1,…,n2JPn(J+1,2J)=κ(nJ+s)λJ+rμJ+sκ(nJ+r+1)λJ+sμJ+r.
Then, by (44) to (45), we know that it is sufficient to show
(46)∑k=1J‍αJ+k+∑k=1J‍κ(nJ+k)μJ+k=∑k=1J‍(1-∑l=1J‍pJ+k,J+l)λJ+k+∑k=1J‍κ(nJ+k)αJ+kμJ+kλJ+k+∑r,s=1J‍κ(nJ+s)λJ+rμJ+spJ+r,J+sλJ+s.
Owing to the routing structure and the traffic equation (3), we have the following observations:
(47)∑k=1J‍(1-∑l=1J‍pJ+k,J+l)λJ+k=∑k=1J‍αJ+k,∑r,s=1J‍κ(nJ+s)λJ+rμJ+spJ+r,J+sλJ+s=∑s=1J‍(κ(nJ+s)μJ+sλJ+s)∑r=1J‍λJ+rpJ+r,J+s=∑s=1J‍κ(nJ+s)μJ+s-∑s=1J‍κ(nJ+s)μJ+sαJ+sλJ+s.
Then, substituting (47) into (46), we get the necessary equality.Next, we derive the initial distributionP0s(J+k) corresponding to state nJ+k=0. By the network stability conditions (12) and ρj<1 for each j∈{1,…,J}, it follows from ∑nJ+k=0∞‍PnJ+ks(J+k)=1 that the following initial distribution is well posed:
(48)P0s(J+k)=11+∑nJ+k=1∞∏rJ+k=1nJ+k(λJ+k/κ(rJ+k)μJ+k).
Hence, we complete the proof of Theorem 1 for priority types h∈{1,2}.Step 3.
To finish the induction procedure, in this step, we first suppose that the result described in Theorem1 is true for all classes with priority rank h∈{1,…,I-1}; then, we show that it is true for all classes with priority rank h∈{1,…,I}. By the similar illustration in getting (41), we have the following differential equations:
(49)P˙n((I-1)J+1,IJ)(t)=-∑k=1J‍α(I-1)J+kPn((I-1)J,IJ)(t)-∑k=1J∑n1,…,n(I-1)J=0∞‍a(nIJ+k,∑u=0I-2‍nuJ+k)μ(I-1)J+kkkkkkklk×(1-p(I-1)J+k,(I-1)J+k)kkkkkklk×Pn(1,(I-1)J)(t)Pn((I-1)J+1,IJ)(t)+∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k+1,∑u=0I-2‍nuJ+k)μ(I-1)J+kkkkkkkkk×(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)Pn(1,(I-1)J)(t)kkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+k+1,…,nIJ(t)+∑k=1J‍(n(I-1)J+k∧1)α(I-1)J+kkkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+k-1,…,nIJ(t)+∑r,s=1,s≠rJ∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+r+1,∑u=0I-2‍nuJ+r)kkkkkkkkkkk×μ(I-1)J+rp(I-1)J+r,(I-1)J+skkkkkkkkkkk×Pn(1,(I-1)J)(t)(n(I-1)J+s∧1)kkkkkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,nIJ(t).
Next, we show that the given distribution corresponding the lowest priority type in Theorem 1 is a steady-state solution of those equations described by (49). It suffices to demonstrate that the derivatives in the above equations are all made zero by setting Pn((I-1)J,IJ)(t)=Pn((I-1)J,IJ), that is, to prove that
(50)∑k=1J‍α(I-1)J+kPn((I-1)J,IJ)+∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k,∑u=0I-2‍nuJ+k)kkkkkk×Pn(1,(I-1)J)μ(I-1)J+kPn((I-1)J+1,IJ)=∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k+1,∑u=0I-2‍nuJ+k)Pn(1,(I-1)J)kklkkkk×μ(I-1)J+k(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)kklkkkk×Pn(I-1)J+1,…,n(I-1)J+k+1,…,nIJ+∑k=1J‍(n(I-1)J+k∧1)α(I-1)J+kPn(I-1)J+1,…,n(I-1)J+k-1,…,nIJ+∑r,s=1,s≠rJ∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+r+1,∑u=0I-2‍nuJ+r)kkkkkkkkkkk×Pn(1,(I-1)J)μ(I-1)J+r×p(I-1)J+r,(I-1)J+s(n(I-1)J+s∧1)×Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,nIJ+∑k=1J∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k,∑u=0I-2‍nuJ+k)Pn(1,(I-1)J)kkklkkkkl×μ(I-1)J+kp(I-1)J+k,(I-1)J+kPn((I-1)J+1,IJ).
By the definition of κ(·) in Theorem 1, we have
(51)κ(n(I-1)J+k)=∑n1,…,n(I-1)J=0∞‍a(n(I-1)J+k,∑u=0I-2‍nuJ+k)Pn(1,(I-1)J).
Then, we can rewrite (50) as follows:
(52)(∑k=1J‍α(I-1)J+k+∑k=1J‍κ(n(I-1)J+k)μ(I-1)J+k)Pn((I-1)J+1,IJ)=∑k=1J‍κ(n(I-1)J+k+1)μ(I-1)J+k(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)kkkkkkk×Pn(I-1)J+1,…,n(I-1)J+k+1,…,nIJ+∑k=1J‍(n(I-1)J+k∧1)α(I-1)J+kPn(I-1)J+1,…,n(I-1)J+k-1,…,nIJ+∑r,s=1,s≠rJ‍κ(n(I-1)J+k+1)kkkkkkkkkkkk×μ(I-1)J+rp(I-1)J+r,(I-1)J+s(n(I-1)J+s∧1)kkkkkkkkkkkk×Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,nIJ+∑k=1J‍κ(n(I-1)J+k)μ(I-1)J+kp(I-1)J+k,(I-1)J+kPn((I-1)J+1,IJ).
From the given distribution in Theorem 1, we have
(53)Pn(I-1)J+1,…,n(I-1)J+k+1,…,n2JPn((I-1)J+1,2J)=λ(I-1)J+kκ(n(I-1)J+k+1)μ(I-1)J+k,Pn(I-1)J+1,…,n(I-1)J+k-1,…,n2JPn((I-1)J+1,2J)=κ(n(I-1)J+k)μ(I-1)J+kλ(I-1)J+k,Pn(I-1)J+1,…,n(I-1)J+r+1,…,n(I-1)J+s-1,…,n2JPn((I-1)J+1,2J)=κ(n(I-1)J+s)λ(I-1)J+rμ(I-1)J+sκ(n(I-1)J+r+1)λ(I-1)J+sμ(I-1)J+r.
From (52) to (53), it is sufficient to show
(54)∑k=1J‍α(I-1)J+k+∑k=1J‍κ(n(I-1)J+k)μ(I-1)J+k=∑k=1J‍(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)λ(I-1)J+k+∑k=1J‍κ(n(I-1)J+k)α(I-1)J+kμ(I-1)J+kλ(I-1)J+k+∑r,s=1J‍κ(n(I-1)J+s)λ(I-1)J+rμ(I-1)J+sp(I-1)J+r,(I-1)J+sλ(I-1)J+s.
Owing to the routing structure and the traffic equation (3), we have the following observations:
(55)∑k=1J‍(1-∑l=1J‍p(I-1)J+k,(I-1)J+l)λ(I-1)J+k=∑k=1J‍α(I-1)J+k,∑r,s=1J‍κ(n(I-1)J+s)λ(I-1)J+rμ(I-1)J+sp(I-1)J+r,(I-1)J+sλ(I-1)J+s=∑s=1J‍(κ(n(I-1)J+s)μ(I-1)J+sλ(I-1)J+s)∑r=1J‍λ(I-1)J+rp(I-1)J+r,(I-1)J+s=∑s=1J‍κ(n(I-1)J+s)μ(I-1)J+s-∑s=1J‍κ(n(I-1)J+s)μ(I-1)J+sα(I-1)J+sλ(I-1)J+s.
Then, substituting (55) into (54), we get the necessary equality.Next, we derive the initial distributionP0s((I-1)J+k) corresponding to state n(I-1)J+k=0. By the network stability conditions (12) and ρj<1 for each j∈{1,…,J}, it follows from ∑n(I-1)J+k=0∞‍Pn(I-1)J+ks((I-1)J+k)=1 that the following initial distribution is well posed:
(56)P0s((I-1)J+k)=(1+∑n(I-1)J+k=1∞‍∏r(I-1)J+k=1n(I-1)J+k‍λ(I-1)J+kκ(r(I-1)J+k)μ(I-1)J+k)-1.
Hence, we complete the proof of Theorem 1.


## 5.2. Proof of Proposition2

Net I. We justify the stability condition (12) for Net I by induction in terms of the number of job types; that is, h=1,2,…,I.First, we consider the case thath=2. Note that cj=1 for all j∈{1,2,…,J} and step one in the proof of Theorem 1, we have that, for nJ+k>0 with k∈{1,2,…,J},
(57)κ(nJ+k)=P0s(k)=1-λkμk.
Therefore, condition (12) is true for h=2. Furthermore, we have
(58)P0s(J+k)=1-λJ+k(1-λk/μk)μJ+k.Second, forh≤I-1, we suppose that
(59)κ(n(h-1)J+k)=1-∑i=0h-2‍λiJ+kμiJ+k.
Hence, we have
(60)P0s((h-1)J+k)=1-λ(h-1)J+k(1-∑i=0h-2λiJ+k/μiJ+k)μ(h-1)J+k.
From the induction assumptions (59) and (60), we know that (12) is true for h≤I-1.Finally, we show that (12) holds for h=I. In fact, from the definition of κ(n(I-1)J+k), we know that, for n(I-1)J+k>0,
(61)κ(n(I-1)J+k)=∏i=0I-2‍P0s(iJ+k)=1-∑i=0I-2‍λiJ+kμiJ+k.
Then, we see that (12) is true for h=I. Hence, we complete the proof of Net I.Net II. Note that, for each k∈{1,…,J} and nJ+k≥cs(k), we have
(62)κ(nJ+k)=∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)=∑nk=0cs(k)-1‍a(nJ+k,nk)Pnks(k)=P0s(k)(cs(k)+(cs(k)-1)(λkμk)+⋯+(λk/μk)cs(k)-1(cs(k)-1)!)=(∑nk=0cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!)×(∑nk=0cs(k)-1‍(λk/μk)nknk!+(λk/μk)cs(k)cs(k)!(1-λk/cs(k)μk))-1=(cs(k)(1-λkcs(k)μk)∑nk=0cs(k)-1‍(cs(k)-nk)(λk/μk)nkcs(k)nk!)×(1+∑nk=0cs(k)-2‍(1(nk+1)!-1cs(k)nk!)(λkμk)nk+1)-1=cs(k)(1-λkcs(k)μk).
Furthermore, for 0<nJ+k<cs(k), we have
(63)κ(nJ+k)=∑n1,…,nJ=0∞‍a(nJ+k,nk)Pn(1,J)=∑nk=0cs(k)-1‍a(nJ+k,nk)Pnks(k)=P0s(k)((λkμk)cs(k)-1nJ+k(1+11!(λkμk)kkkkkkkkkkkkkkk+⋯+1(cs(k)-nJ+k)!(λkμk)cs(k)-nJ+k)kkkkkkkkl+nJ+k-1(cs(k)-nJ+k+1)!(λkμk)cs(k)-nJ+k+1kkkkkkkkl+⋯+1(cs(k)-1)!(λkμk)cs(k)-1)=(nJ+k∑nk=0cs(k)-nJ+k‍(λk/μk)nknk!kkkkkl+∑nk=cs(k)-nJ+k+1cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!)×(∑nk=0cs(k)-1‍(λk/μk)nknk!+(λk/μk)cs(k)cs(k)!(1-λk/cs(k)μk))-1=(+∑nk=cs(k)-nJ+k+1cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!nJ+k)nJ+k(1-λkcs(k)μk)kkkkk×(∑nk=0cs(k)-nJ+k‍(λk/μk)nknk!kkkkkkkkkkk+∑nk=cs(k)-nJ+k+1cs(k)-1‍(cs(k)-nk)(λk/μk)nknk!nJ+k))×(1+∑nk=0cs(k)-2‍(1(nk+1)!-1cs(k)nk!)(λkμk)nk+1)-1≥nJ+k(1-λkcs(k)μk).
Then, we can see that (12) is true. Hence, we complete the proof of Proposition 2.


## 5.3. Illustration of Example3

For each classk∈{1,2,…,J}, the initial distribution P0s(k) can be calculated as follows, since ρs(k)<1:
(64)P0s(k)=1-λk/2μk1+λk/2μk.
Then, for nJ+k=1, we have
(65)κ(nJ+k)=P0s(k)(1+λkμk)=(1-λk/2μk)(1+λk/μk)1+λk/2μk≥1-λk2μk.
Furthermore, for nJ+k≥2, we have
(66)κ(nJ+k)=P0s(k)(2+λkμk)=2(1-λk2μk).
Therefore, since ρs(k)<1, we have
(67)P0s(J+k)=(1+λJ+kκ(1)μJ+kk+λJ+kκ(1)μJ+kλJ+k/μJ+k(1-λJ+k/κ(2)μJ+k)κ(2)μJ+k)-1=(1-λk/2μk-λJ+k/2μJ+k)(1+λk/μk)1+λk/2μk+λJ+k/2μJ+k-(1/2)(λk/μk)2.
Finally, we can compute κ(n2J+k) for n2J+k≥2 as follows:
(68)κ(n2J+k)=2P0s(k)P0s(J+k)+P1s(k)P0s(J+k)+P0s(k)P1s(J+k)=P0s(k)P0s(J+k)(2+λkμk+λJ+kμJ+k)=(2(1-λk2μk-λJ+k2μJ+k)(1+λkμk)(1-λk2μk)kkkkk×(1+λk2μk+λJ+k2μJ+k))×((1+λk2μk+λJ+k2μJ+k-12(λkμk)2)(1+λk2μk))-1=(2(1-λk2μk-λJ+k2μJ+k)kkkkk×(1+λkμk+λJ+k2μJ+k-(λk2μk)2-14(λkμk)3kkkkkkkk+2(1-λk2μk)λkλJ+k4μkμJ+k))×(1+λkμk+λJ+k2μJ+k-(λk2μk)2kkkkkkk-14(λkμk)3+λkλJ+k4μkμJ+k)-1.
Thus, if λk/μk≤1, from the condition ρs(k)<1, we see that
(69)λ2J+kκ(n2J+k)μ2J+k≤λ2J+k2(1-λk/2μk-λJ+k/2μJ+k)μ2J+k<1;
that is, the condition (12) holds.


## 6. Conclusion

The research conducted in the paper is to iteratively derive the product-form solutions of stationary distributions for a particular type of preemptive priority multiclass queueing networks with multiserver stations. This type of queueing systems can typically be used to model the stochastic dynamics of some large scale backbone networks with multiprocessor shared-memory switches or local (edge) cloud computing centers with parallel-server pools. The queueing networks are Markovian with exponential interarrival and service time distributions. The obtained iterative solutions can be used to conduct performance analysis or as comparison criteria for approximation and simulation studies. Numerical comparisons with existing Brownian approximating model related to general interarrival and service times are provided to show the effectiveness of our current designed algorithm and our former derived BAM. Furthermore, based on the iterative solutions, we also give some analysis concerning network stability for some cases of these queueing systems, which provides some insight for more general study.



---

*Source: 767651-2014-10-13.xml*