Asymptotic product-form steady-state for generalized Jackson networks in multi-scale heavy traffic

Jackson (1957) introduced a class of first-come-first-serve (FCFS) queueing networks that are known today as open Jackson networks. The defining characteristics of a Jackson network are (a) all customers visiting a service station are homogeneous in terms of the service time distribution and the routing probabilities, and (b) all interarrival and service time distributions, referred to as the primitive distributions, are exponential. For a Jackson network with $J$ stations, the queue length vector process is a continuous time Markov chain (CTMC). Jackson’s pioneering contribution is that when the traffic intensity at each station is less than one, the CTMC has a unique stationary distribution that is of product form, meaning that the steady-state queue lengths at the various stations in the network are independent. This product-form result makes the computation of steady-state performance measures highly tractable.

When the primitive distributions are general, the corresponding network is known as a generalized Jackson network (GJN). For a GJN, the product-form result no longer holds in general. Creating and justifying approximations for GJNs have been an active research area for more than 60 years. In this paper, we prove the following asymptotic product-form result for a GJN: for a GJN with $J$ stations, satisfying the multi-scale heavy traffic condition

$\displaystyle 1\ll\frac{1}{1-\rho_{1}}\ll\frac{1}{1-\rho_{2}}\ll\cdots\ll\frac{1}{1-\rho_{J}},$

(1.1)

the following approximation holds

$\displaystyle\Big{(}(1-\rho_{1})Z_{1},\ldots,(1-\rho_{J})Z_{J}\Big{)}\approx\left(\rho_{1}d_{1}E_{1},\ldots,\rho_{J}d_{J}E_{J}\right),$

(1.2)

where $E_{1},\ldots,E_{J}$ are i.i.d exponential random variables with mean $1$, and $(d_{1},\ldots,d_{J})>0$ is some deterministic vector. In (1.2), $\rho_{j}<1$ is the traffic intensity at station $j$ defined in (2.4), $Z=(Z_{1},\ldots,Z_{J})$ is the vector of steady-state queue length. The right side of (1.2) indicates a product-form distribution. Regarding the multi-scale heavy traffic condition (1.1), $1/(1-\rho_{j})$ will be called the heavy traffic normalization factor for station $j$ and $a\ll b$ is interpreted as “$b$ being much larger than $a$”. Condition (1.1) says that the heavy traffic normalization factors for different stations are all large, but of widely separated magnitudes. A precise definition of the multi-scale heavy traffic condition will be defined in (3.5) in Section 3, and a precise version of (1.2) in terms of convergence in distribution will be stated in Theorem 3.1.

Our result is the first of its kind. It has many potential ramifications. An obvious benefit of our result is that the right side of (1.2) offers a scalable method to approximately evaluate the steady-state performance of a GJN. In Section 4, we demonstrate that approximations based on (1.2) can be accurate sometimes even when the load-separation condition in (1.1) is violated. Our asymptotic product-form result has been extended to multiclass queueing networks operating under static buffer priority service policies in Dai and Huo (2024). An active research direction is to explore whether the asymptotic product-form phenomenon holds for a much broader class of stochastic systems including stochastic processing networks as defined in Dai and Harrison (2020). These exciting developments have the potential to radically speed up the performance analysis and optimization of real world complex systems.

When a GJN is in conventional heavy traffic, namely all the stations have roughly equal heavy traffic normalization factors, Reiman (1984) and Johnson (1983) prove that the scaled queue length vector process converges to a multi-dimensional reflecting Brownian motion (RBM). In their conventional heavy traffic limit theorem, one common scaling is applied (based on the heavy traffic normalization factor for an arbitrarily chosen station) across all the different stations. Their RBMs are first introduced in Harrison and Reiman (1981), and the stationary distributions of these RBMs were shown to exist in Harrison and Williams (1987a) when $\rho_{j}<1$ for each station $j$. These pioneering works have propelled the development of Brownian models for queueing network performance analysis and controls in the last forty years; see, for example, Harrison (1988), Harrison and Nguyen (1993), Chen and Yao (2001), Kang et al. (2009), Ata et al. (2024).

The stationary distribution of a multi-dimensional RBM is typically not of product-form except under a “skew-symmetry” condition characterized in Harrison and Williams (1987b); see, e.g., Harrison (1978). Numerical algorithms have been developed for computing the stationary distributions of RBMs in low dimensions; see Dai and Harrison (1992) and Shen et al. (2002). It is likely that these algorithms, even when convergent, do not scale well when the number $J$ of stations gets large. Moreover, a recent simulation-based algorithm proposed in Blanchet et al. (2021) has a complexity that scales linearly in dimension $J$ when the RBM data satisfies (a) a uniform (in $J$) contraction condition and (b) a uniform (in $J$) stability condition. The requirement of uniformity in $J$ is a strong condition when $J$ is large. In other words, the complexity guarantee of their work may not cover many real-world high-dimensional network examples. On the other hand, the approximation based on (1.2) can be easily computed even when $J$ is (very) large. The approximation (1.2) can be either used directly to approximate the original networks or it can be used at a “pre-conditioner” for numerical algorithms, such as that of Dai and Harrison (1992). Specifically, their algorithm computes a density relative to a reference measure which must be chosen well in order to guarantee convergence. The approximation (1.2) is a possible choice for that reference measure; its effectiveness will be explored elsewhere.

Our proof of Theorem 3.1 employs the BAR-approach recently developed in Braverman et al. (2017, 2024) for continuous-time queueing networks with general primitive distributions. The BAR-approach was developed initially for providing an alternative to the limit interchange procedure in Gamarnik and Zeevi (2006), Budhiraja and Lee (2009), Gurvich (2014), and Ye and Yao (2018), justifying the convergence of the prelimit stationary distribution to the RBM stationary distribution under conventional heavy traffic. Our proof of Theorem 3.1 relies heavily on a steady-state moment bound for the scaled queue length vector. This bound is proved in a companion paper Guang et al. (2023).

The parameter $d_{j}$ in (1.2) has an explicit formula (3.9) in terms of the first two moments of the primitive distributions and the routing probabilities. The distribution of $d_{j}E_{j}$ corresponds to the stationary distribution of a one-dimensional RBM approximating a single station queue. The “single station” can be constructed precisely and has a good intuitive interpretation as described immediately below the statement of Theorem 3.1, consistent with the bottleneck analysis advanced in Chen and Mandelbaum (1991). Our work contributes to a line of research pioneered by Whitt (1983) on decomposition approaches for the performance analysis of GJNs, and partly justifies the sequential bottleneck decomposition (SBD) heuristic proposed in Dai et al. (1994). For a recent survey of decomposition approaches including “robust queueing approximations”, readers are referred to Whitt and You (2022). We leave it as a future work to compare the approximation (1.2) with the extensive numerical work in Whitt and You (2022).

We note that when the normalization factors $(1-\rho_{j})^{-1}$ are widely separated, this implies that in heavy traffic, the diffusion temporal scalings $(1-\rho_{j})^{-2}$ that dictate the time scales governing the dynamics at each station are also widely separated. It is this “time scale separation” that leads to the independence across stations. A precise statement of this mathematical result and its proof has been developed in Chen et al. (2025).

We first define an open generalized Jackson network, following closely Section 2.1 of Braverman et al. (2017) both in terms of terminology and notation. There are $J$ single-server stations in the network. Each station has an infinite-size buffer that holds jobs waiting for service. Each station potentially has an external job arrival process. An arriving job is processed at various stations in the network, one station at a time, before eventually exiting the network. When a job arrives at a station and finds the server busy processing another job, the arriving job waits in the buffer until its turn for service. After being processed at one station, the job either is routed to another station for further processing or exits the network, independent of all previous history. Jobs at each station, either waiting or in service, are homogeneous in terms of service times and routing probabilities. Jobs at each station are processed following the first-come-first-serve (FCFS) discipline.

Let ${\cal J}=\{1,\ldots,J\}$ be the set of stations. The following notational system follows Braverman et al. (2024) closely. For each station $j\in\cal{J}$, we are given two nonnegative real numbers $\alpha_{j}\in\mathbb{R}_{+}$ and $\mu_{j}>0$, two sequences of i.i.d. random variables $T_{e,j}=\{T_{e,j}(i),i\in\mathbb{N}\}$ and $T_{s,j}=\{T_{s,j}(i),i\in\mathbb{N}\}$, and one sequence of i.i.d. random vectors $\Phi_{j}=\{\Phi_{j}(i),i\in\mathbb{N}\}$, where $\mathbb{N}$ denotes the set of natural numbers. We allow the possibility that some stations may not have external arrivals by assuming $\alpha_{j}=0$. If such a station $j$ exists, all terms associated with external arrivals at station $j$ can be removed without ambiguity. To maintain clean notation, we assume each station $j\in\cal{J}$ has external arrival. We assume $T_{e,1}$, …, $T_{e,J}$, $T_{s,1}$, …, $T_{s,J}$, $\Phi_{1}$, $\ldots$, $\Phi_{J}$ are independent. Following Braverman et al. (2024), we assume $\mathbb{E}[T_{e,j}(1)]=1$, $\mathbb{E}[T_{s,j}(1)]=1$ such that $T_{e,j}(1)$ and $T_{s,j}(1)$ are unitized. For $j\in\cal{J}$, we use $T_{e,j}(i)/\alpha_{j}$ to represent the interarrival time between the $(i-1)$th and $i$th external arriving jobs at station $j$, and $T_{s,j}(i)/\mu_{j}$ to denote the $i$th job service time at station $j$. Therefore, $\alpha_{j}$ is interpreted as the external arrival rate and $\mu_{j}$ as the service rate at station $j$. The random vector $\Phi_{j}(1)$ takes vector $e^{(k)}$ with probability $P_{jk}$ for $k\in\mathcal{J}$ and takes vector $e^{(0)}$ with probability $P_{j0}\equiv 1-\sum_{k\in\mathcal{J}}P_{jk}$, where $e^{(k)}$ is the $J$-vector with component $k$ being 1 and all other components being $0$, and $e^{(0)}$ is the $J$-vector of zeros. When $\Phi_{j}(i)=e^{(k)}$, the $i$th departing job from station $j$ goes next to station $k\in\mathcal{J}$. When $\Phi_{j}(i)=e^{(0)}$, the job exits the network. We assume the $J\times J$ routing matrix $P=(P_{ij})$ is transient, which is equivalent to requiring that $(I-P)$ is invertible. This assumption ensures that each arriving job will eventually exit the network, and thus the network is known as an open generalized Jackson network.

Denote by $c^{2}_{e,j}=\operatorname{Var}\left(T_{e,j}(1)\right)$ and $c^{2}_{s,j}=\operatorname{Var}\left(T_{s,j}(1)\right)$. It is known that $c^{2}_{e,j}$ and $c^{2}_{s,j}$ are the squared coefficients of variation for interarrival time and service time, respectively.

We consider a family of generalized Jackson networks indexed by $r\in(0,1)$. We denote by $\mu^{(r)}_{j}$ the service rate at station $j$ in the $r$th network. To keep the presentation clean, we assume the service rate vector $\mu^{(r)}$ is the only network parameter that depends on $r\in(0,1)$, and that the service time $T_{s,j}/\mu^{(r)}_{j}$ depends on $r$ solely through the service rate. The arrival rate $\alpha$, the unitized interarrival and service times, and routing vectors do not depend on $r$. Consequently, the nominal total arrival rate $\lambda$ does not depend on $r$.

Condition (3.5) is equivalent to

$\displaystyle 1-\rho^{(r)}_{j}=b_{j}r^{j}/\mu^{(r)}_{j},\quad j\in\mathcal{J}.$

Condition (3.5) implies $\mu^{(r)}_{j}\to\lambda_{j}$ and

$\displaystyle 1/(1-\rho_{1}^{(r)})=\frac{\mu^{(r)}_{1}}{b_{1}}\frac{1}{r}\to\infty,\quad\frac{1/(1-\rho_{j+1}^{(r)})}{1/(1-\rho_{j}^{(r)})}=\frac{b_{j}\mu^{(r)}_{j+1}}{b_{j+1}\mu^{(r)}_{j}}\frac{1}{r}\to\infty,\quad j=1,\ldots,J-1.$

as $r\downarrow 0$, which is consistent with (1.1). We call condition (3.5) the multi-scale heavy traffic condition. Condition (3.5) is one way to achieve “widely separated traffic normalization factors”. Clearly, it is not the only way. We choose to use the form in condition (3.5) to make our proofs clean. We leave it to future research to understand the most general form of load-separation under which Theorem 3.1 still holds. In the rest of this paper, we choose $b_{j}=1$ for $j\in\mathcal{J}$. When $b_{j}\neq 1$, all the proofs remain essentially the same.

The $(J+1+\delta_{0})$th moment assumption is used to derive in Section 7 the asymptotic basic adjoint relationship (7.75) and Taylor expansions (7.72) and (7.73). These three equations are critical in the proof of Theorem 3.1.

The following assumption is mild. Dai (1995) proves that the following assumption holds when $\rho_{j}^{(r)}<1$ for each station $j$, and the interarrival times are “unbounded and spread-out”.

Under Assumption 3.3, we use

$\displaystyle X^{(r)}=\big{(}Z^{(r)},R_{e}^{(r)},R_{s}^{(r)}\big{)}$

to denote the steady-state random vector that has the stationary distribution. It is well known that

$\displaystyle\mathbb{P}\{Z^{(r)}_{j}=0\}=1-\rho^{(r)}_{j}=r^{j}/\mu^{(r)}_{j},\quad j\in\mathcal{J}.$

(3.7)

See, for example, (A.11) and (A.13) of Braverman et al. (2017) for a proof. Therefore, under multi-scale heavy traffic condition (3.5), each server approaches $100\%$ utilization.

The following is the main result of this paper.

As discussed in the introduction, the pre-limit stationary distribution (the left side of 3.8) is not of product form, nor is stationary distribution of conventional heavy traffic limit Reiman (1984). It is somewhat surprising that the multi-scale limit $E_{1},\ldots,E_{J}$ are independent. Theorem 3.1 suggests an obvious approximation to the stationary distribution of the queue length vector process for a generalized Jackson network with widely separated values of $(\mu_{j}-\lambda_{j})^{-1}$, $j\in\mathcal{J}$. Our approximation is

$\displaystyle\big{(}Z_{1}(\infty),Z_{2}(\infty),\ldots,Z_{J}(\infty)\big{)}\stackrel{{\scriptstyle d}}{{\approx}}\left(\frac{\rho_{1}}{1-\rho_{1}}d_{1}E_{1},\frac{\rho_{2}}{1-\rho_{2}}d_{2}E_{2},\ldots,\frac{\rho_{J}}{1-\rho_{J}}d_{J}E_{J}\right),$

(3.11)

where $(E_{1},\ldots,E_{J})$ is described in Theorem 3.1, and $\stackrel{{\scriptstyle d}}{{\approx}}$ denotes “has approximately the same distribution as” (and carries no rigorous meaning beyond that imparted by Theorem 3.1).

We now define the $J\times J$ matrix $w=(w_{ij})$ used in (3.9) and (3.10). For each $j\in\mathcal{J}$,

		$\displaystyle(w_{1j},\ldots,w_{j-1,j})^{\prime}=(I-P_{j-1})^{-1}P_{[1:j-1],j},$		(3.12)
		$\displaystyle(w_{jj},\ldots,w_{J,j})^{\prime}=P_{[j:J],j}+P_{[j:J],[1:j-1]}(I-P_{j-1})^{-1}P_{[1:j-1],j},$		(3.13)

where prime denotes the transpose, and, for $A,B\subseteq\mathcal{J}$, $P_{A,B}$ is the submatrix of $P$ whose entries are $(P_{k\ell})$ with $k\in A$ and $\ell\in B$ with the following conventions: for $\ell\leq k$,

$\displaystyle[\ell:k]=\{\ell,\ell+1,\ldots,k\},\quad P_{A,k}=P_{A,\{k\}},\quad P_{\ell,B}=P_{\{\ell\},B},\quad P_{k}=P_{[1:k],[1:k]},$

and any expression involving $P_{A,B}$ is interpreted to be zero when either $A$ or $B$ is the empty set. The expressions in (3.12) and (3.13) may look complicated, but $(w_{ij})$ has the following probabilistic interpretation. Recall that the routing matrix $P$ can be embedded within a transition matrix of a DTMC on state space $\{0\}\cup\mathcal{J}$ with $0$ being the absorbing state. For each $i\in\mathcal{J}$ and $j\in\mathcal{J}$, $w_{ij}$ is the probability that starting from state $i$, the DTMC will eventually visit state $j$ while avoiding visiting states $\{0,j+1,\ldots,J\}$. See more discussion of $(w_{ij})$ in Lemma A.1 in Section 5.

The variance parameter $\sigma_{j}^{2}$ in (3.10) has the following appealing interpretation. Consider a renewal arrival process $\{E(t),t\geq 0\}$ with arrival rate $\nu$ and squared coefficient of variation $c^{2}$ of the interarrival time. Each arrival has probability $p_{j}$ of going to station $j$ instantly, independent of all previous history. Let $\{A_{j}(t),t\geq 0\}$ be the corresponding arrival process to station $j$. It is known that

$\displaystyle\lim_{t\to\infty}\frac{\text{Var}(A_{j}(t))}{t}=\bigl{(}c^{2}p_{j}^{2}+p_{j}(1-p_{j})\bigr{)}\nu.$

(3.14)

Thus, each of the first $J$ terms in the right side of (3.10) corresponds to the variance parameter in (3.14) of $J$ arrival sources to station $j$. The first $j-1$ terms correspond to arrival sources from lightly-loaded stations, the $j$th term corresponds to the external arrival process to station $j$, and the last $J-j$ terms correspond to arrival sources from service completions at heavily-loaded stations. For a lightly-loaded station $\ell<j$, its external arrival process serves as the renewal arrival process, each arrival having probability $w_{\ell j}$ to join station $j$. For a heavily-loaded station $k>j$, its service process serves as the renewal arrival process (as if the server is always busy), each arrival having probability $w_{kj}$ to join station $j$. The resulting formula $d_{j}$ in (3.9) is consistent with the diffusion approximation of a $\sum_{i}G_{i}/GI/1$ queue with $J$ external arrival processes and immediate feedback probability $w_{jj}$ after the service at station $j$.

Theorem 3.1 will be proved in Section 5. The proof relies on the BAR in Section 6 and the asymptotic BAR in Section 7. In addition, the proof relies critically on a uniform moment bound that is proved in the companion paper (Guang et al., 2023, Theorem 2). For easy reference, we state that result here.

The moment bound (3.15) holds of course when interarrival and service time distributions are exponential because $Z_{j}^{(r)}$ is geometrically distributed with mean $\rho^{(r)}_{j}/(1-\rho^{(r)}_{j})$. In the exponential case, Lemma 23 of Xu (2022) provides an alternative proof via Lyapunov functions.

In this section, we illustrate the value of approximation (3.11) that is rooted in Theorem 3.1. In Section 4.1, we propose an approximation formula for the cumulative distribution function (cdf) of the queue length vector. In Section 4.2, we perform a simulation study on a two-station queueing network to demonstrate the accuracy of the approximation. We experiment with two sets of network parameters. In one set, we vary the level of load separation between stations, and in the other, we vary the variability of primitive distributions. In all cases, our proposed approximation performs well in multiple performance metrics, even in the case when the load separation condition is grossly violated. The latter phenomenon deserves a future study.

In this section, we present a proof of Theorem 3.1 under an additional assumption (see Assumption 5.1 below). The proof relies on a transform version of the basic adjoint relationship (BAR), which will be developed in Section 6. Recall that $X^{(r)}=(Z^{(r)},R_{e}^{(r)},R^{(r)}_{s})$ denotes the steady state of the Markov process that governs the dynamics of the GJN. Define the moment generating function (MGF), or Laplace transform, $\phi^{(r)}$ of $Z^{(r)}$ as

$\displaystyle\phi^{(r)}(\theta)=\mathbb{E}\Big{[}\exp{\Big{(}\sum_{j\in\mathcal{J}}\theta_{j}Z_{j}^{(r)}\Big{)}}\Big{]},\quad\theta\in\mathbb{R}^{J}_{-}\equiv\{y\in\mathbb{R}^{J},y\leq 0\}.$

Then $\phi^{(r)}(r\eta_{1},\ldots,r^{J}\eta_{J})$, with $\eta=(\eta_{1},\ldots,\eta_{J})^{\prime}\in\mathbb{R}^{J}_{-}$, is the MGF of $(rZ^{(r)}_{1},\ldots,r^{J}Z^{(r)}_{J})$. It is well known that proving (3.8) is equivalent to proving

$\displaystyle\lim_{r\to\infty}\phi^{(r)}(r\eta_{1},\ldots,r^{J}\eta_{J})=\prod_{j\in\mathcal{J}}\frac{1}{1-d_{j}\eta_{j}},\quad\eta=(\eta_{1},\ldots,\eta_{J})^{\prime}\in\mathbb{R}^{J}_{-}.$

(5.20)

To prove (5.20), it suffices to prove the following lemma. Equation (5.20) then follows from 5.21 by induction.

We will present two proofs for Proposition 5.1. The first proof is given in this section under the following bounded-support assumption.

Assumption 5.1 is stronger than the moment assumption, Assumption 3.2. The rest of this section is devoted to the proof of Proposition 5.1 under Assumption 5.1. This proof relies on a transform version of the basic adjoint relationship (BAR) (5.27), which will be proved in Section 6. In Section 7, we will present the second proof of Proposition 5.1 without Assumption 5.1.

In the rest of this section, instead of working with $\phi^{(r)}(\theta)$ directly, we will work with

$\displaystyle{\psi}^{(r)}(\theta)=\mathbb{E}[{f}_{\theta}(X^{(r)})],\quad\text{ and }\quad{\psi}^{(r)}_{i}(\theta)=\mathbb{E}[f_{\theta}(X^{(r)})\mid Z^{(r)}_{i}=0]\quad i\in\mathcal{J},$

(5.23)

for each $\theta\in\mathbb{R}^{J}_{-}$, where for $x=(z,u,v)\in\mathbb{Z}^{J}_{+}\times\mathbb{R}^{J}_{+}\times\mathbb{R}^{J}_{+}$,

$\displaystyle{f}_{\theta}(x)=\exp\Big{(}\langle\theta,z\rangle-\sum_{i}\alpha_{i}\gamma_{i}(\theta_{i})u_{i}-\sum_{i}\mu^{(r)}_{i}\xi_{i}(\theta)v_{i}\Big{)}.$

(5.24)

In (5.24), $\gamma_{i}(\theta)$ and $\xi_{i}(\theta)$ are defined via

	$\displaystyle e^{\theta_{i}}\mathbb{E}[e^{-\gamma_{i}(\theta_{i})T_{e,i}(1)}]=1,$		(5.25)
	$\displaystyle\bigl{(}P_{i0}e^{-\theta_{i}}+\sum_{j\in\mathcal{J}}P_{ij}e^{\theta_{j}-\theta_{i}}\bigr{)}\mathbb{E}[e^{-\xi_{i}(\theta)T_{s,i}(1)}]=1.$		(5.26)

The advantage of working with ${\psi}^{(r)}$ is that it, together with $({\psi}^{(r)}_{1},\ldots,{\psi}^{(r)}_{J})$, satisfies the following transform version of BAR.

The fact that $\gamma_{i}(\theta_{i})$ and $\xi_{i}(\theta)$ are well defined for all $\theta\in\mathbb{R}^{J}$ follows from Lemma 4.1 of Braverman et al. (2017) and condition (5.22). The transform-BAR (5.27) follows from BAR (3.20) of Braverman et al. (2017) by taking $f(x)=f_{\theta}(x)$ in (5.24). The proof of Lemma 5.1 will be recapitulated in Section 6.

When the primitive distributions are exponential, $\gamma_{i}(\theta_{i})$ and $\xi_{i}(\theta)$ have the following explicit expressions: for each $\theta\in\mathbb{R}^{J}$ and $i\in\mathcal{J}$

	$\displaystyle\gamma_{i}(\theta_{i})=e^{\theta_{i}}-1,$		(5.29)
	$\displaystyle\xi_{i}(\theta)=P_{i0}\big{(}e^{-\theta_{i}}-1\big{)}+\sum_{j\in\cal{J}}P_{ij}\big{(}e^{\theta_{j}-\theta_{i}}-1\big{)}.$		(5.30)

Although $\gamma_{i}(\theta_{i})$ and $\xi_{i}(\theta)$ are defined implicitly in general through (5.25) and (5.26), the following Taylor expansions are what we need to access these implicit functions. As $\theta\to 0$,

	$\displaystyle\gamma_{i}(\theta_{i})=\theta_{i}+\frac{1}{2}\gamma^{*}_{i}(\theta_{i})+o(\theta_{i}^{2}),$		(5.31)
	$\displaystyle\xi_{i}(\theta)=\bar{\xi}_{i}(\theta)+\frac{1}{2}\xi^{*}_{i}(\theta)+o(\lvert\theta\rvert^{2}),$		(5.32)

where $\lvert\theta\rvert=\sum_{i}\lvert\theta_{i}\rvert$, and

	$\displaystyle\gamma_{i}^{*}(\theta_{i})=c_{e,i}^{2}\theta_{i}^{2},$		(5.33)
	$\displaystyle\bar{\xi}_{i}(\theta)=\sum_{j\in\mathcal{J}}P_{ij}\theta_{j}-\theta_{i},$		(5.34)
	$\displaystyle\xi^{*}_{i}(\theta)=c_{s,i}^{2}\Big{(}\sum_{j\in\mathcal{J}}P_{ij}\theta_{j}-\theta_{i}\Big{)}^{2}+\sum_{j\in\mathcal{J}}P_{ij}\theta_{j}^{2}-\Big{(}\sum_{j\in\mathcal{J}}P_{ij}\theta_{j}\Big{)}^{2}.$		(5.35)

In (5.31) and (5.32), we adopt the standard notation that $f(r)=o(g(r))$ means $f(r)/g(r)\to 0$ as $r\to 0$, and in (5.33) and (5.35), $c^{2}_{e,i}$ and $c^{2}_{s,i}$ are the variances of the unitized random variables $T_{e,i}$ and $T_{s,i}$, respectively. The $c^{2}_{e,i}$ and $c^{2}_{s,i}$ are also equal to the squared coefficient of variations (SCVs) of the corresponding interarrival and service time distributions. A direct Taylor expansion of the functions in (5.29) and (5.30) are consistent with the expansions in (5.31) and (5.32) with $c^{2}_{e,i}=1$ and $c^{2}_{s,i}=1$.