I have been working so long to figure out the variance of the queue size L for M/GI/1 queue. From many resources, I can find the Expected Queue Size L, but for some reason I cannot find any textbook/articles/papers deriving the Variance of L.
I know the expected queue size for M/G/1 is as follows: $$L_q = \lambda \left( \bar{p}\cdot \frac{\rho}{1-\rho}\cdot \frac{1+\frac{\text{Var}(p)}{\bar{p}^2}}{2}\right)$$
where $\lambda$ is an arrival rate, $\rho$ is a service utilization rate ($=\lambda / \mu$), and $\mu=1/\bar{p}$ is a service rate. ($\bar{p}$ is a mean service time.)
This can be derived by Khintchine-Pollaczek Theorem and Little's Law. Khintchine-Pollaczek Theorem is the following: $$E(W_q)=\bar{p}\cdot \frac{\rho}{1-\rho}\cdot \frac{1+C^2(p)}{2}$$
where
$$C^2(p) = \frac{\text{Var}(p)}{\bar{p}^2}$$.
Could anyone tell me how I can derive the variance of $L_q$ from here? Or, if anyone could give me an equation of $Var(L_q)$, that would help me a lot.
Thank you so much for your help.
Because the number of customers in a $M/G/1$ queue at a departure is equal to the number of customers that arrived during the departing customer's sojourn time, we can write the probability generating function for the number of customers in terms of the Laplace transform of the sojourn time T: \begin{align} G_N(z) &= \sum_{k=0}^\infty z^k \int_0^\infty \frac{(\lambda x)^k}{k!}e^{-\lambda x}f_T(x)\ \mathsf dx\\ &= \int_0^\infty \left(\sum_{k=0}^\infty \frac{(\lambda z x)^k}{k!}\right)e^{-\lambda x}f_T(x)\ \mathsf dx\\ &= \int_0^\infty e^{-\lambda(1-z)x}f_T(x)\ \mathsf dx\\ &= L_T(\lambda(1-z). \end{align} It follows from basic properties of probability generating functions and Laplace transforms that $$G_N^{(k)}(1^{-}) = \mathbb E\left[\prod_{i=0}^{k-1} (N-i)! \right] = (-1)^k L_T^{(k)}(0)\lambda^k = \lambda^k \mathbb E[T^k].$$ Note that $k=1$ corresponds to Little's law, i.e. $\mathbb E[N] = \lambda \mathbb E[T].$ Takács (1962) gives the following recurrence for the moments of the waiting time: $$ \mathbb E[W^k] = \frac\lambda{1-\rho}\sum_{i=1}^k\binom ki \frac{\mathbb E[S^{i+1}]}{i+1}\mathbb E[W^{k-i}], $$ and as the sojourn time $T=W+S$ is the sum of the waiting time and service time (which are independent), it follows that $$ \mathbb E[T^k] = \sum_{i=0}^k \binom ki \mathbb E[W^i]\mathbb E [S^{k-i}]. $$ As $\mathbb E[N] =\lambda \mathbb E[T]$ and $\mathbb E[N(N-1))] = \lambda^2\mathbb E[T^2]$, we may compute the variance of the number of customers in the system: $$ \mathsf{Var}(N) = \mathbb E[N(N-1)]+\mathbb E[N]-\mathbb E[N^2] = \frac{\lambda\mathbb E[S^3]}{3(1-\rho)} + \left(\frac{\lambda\mathbb E[S^2]}{2(1-\rho)}\right)^2 + \frac{\lambda(3-2\rho)\mathbb E[S^2]}{2(1-\rho)} + \rho(1-\rho). $$ Finally, from $$\mathbb E[Q^2] = \sum_{k=1}^\infty (k-1)^2\pi_k = \sum_{k=1}^\infty k^2\pi_k -2\sum_{k=1}^\infty k\pi_k + \sum_{k=1}^\infty \pi_k = \mathbb E[N^2]-2\mathbb E[N] + \rho $$ we may compute the variance of the queue length by $$ \mathsf{Var}(Q) = \frac{\lambda\mathbb E[S^3]}{3(1-\rho)} + \left(\frac{\lambda\mathbb E[S^2]}{2(1-\rho)}\right)^2 + \frac{\lambda\mathbb E[S^2]}{2(1-\rho)}. $$