Closed form for probabilities from generating function of queue length in M/D/1 queue. Let $\pi$ be the stationary distribution for a M/D/1 queue with unit service time and $\lambda<1$. From the Pollaczek-Khinchin transform we know that the generating function for $\pi$ is
$$
G(z) := \sum_{n=0}^\infty \pi_n z^n = \frac{(1-\lambda)(1-z)e^{-\lambda(1-z)}}{e^{-\lambda(1-z)}-z}$$
I'd like to find a closed form for $\pi_n$. So far I have
\begin{align}
G(z) &= \frac{(1-\lambda)(1-z)e^{-\lambda(1-z)}}{e^{-\lambda(1-z)}-z}\\
&= \frac{(1-\lambda)(1-z)}{1-ze^{\lambda(1-z)}}\\
&= (1-\lambda)(1-z)\sum_{n=0}^\infty e^{\lambda(1-z)n}z^n
\end{align}
Section 3.5.2 in these notes, suggests that a "lengthy" calculation yields
$$\pi_n = (1-\lambda)\sum_{i=0}^n e^{i\lambda}(-1)^{n-i}\frac{(i\lambda + n-i)(i\lambda)^{n-i-1}}{(n-i)!}. $$
I tried expanding $e^{\lambda(1-z)n}$ in the series $\sum_{k=0}^\infty\frac{(\lambda(1-z)n)^k}{k!}$ and invoking Fubini but ended up with the term
$$\sum_{n=0}^\infty n^k z^n $$
which doesn't appear to have a nice closed form. Any suggestions on how to proceed?
 A: Starting from $$G(z) = \frac{(1-\lambda)(1-z)e^{-\lambda(1-z)}}{e^{-\lambda(1-z)}-z},$$ we have
\begin{align}
G(z) &= \frac{(1-\lambda)(1-z)}{1-ze^{\lambda(z-1)}}\\
&= (1-\lambda)(1-z)\sum_{j=0}^\infty \left(ze^{\lambda(z-1)}\right)^j\\
&= (1-\lambda)(1-z)\sum_{j=0}^\infty e^{j\lambda(z-1)}z^j\\
&= (1-\lambda)(1-z)\sum_{j=0}^\infty e^{j\lambda} \sum_{k=0}^\infty \frac{(-j\lambda z)^k}{k!}z^j\\
&= (1-\lambda)(1-z)\sum_{j=0}^\infty \sum_{k=0}^\infty e^{j\lambda} (-1)^k \frac{(j\lambda)^k}{k!}z^j\\
&= (1-\lambda)(1-z)\sum_{j=0}^\infty \sum_{n=j}^\infty e^{j\lambda} (-1)^{n-j} \frac{(j\lambda)^{n-j}}{(n-j)!}z^n\\
&=(1-\lambda)(1-z)\sum_{n=0}^\infty\sum_{j=0}^\infty e^{j\lambda} (-1)^{n-j} \frac{(j\lambda)^{n-j}}{(n-j)!}z^n\\
&(1-\lambda)\left[ \sum_{n=0}^\infty\sum_{j=0}^n e^{j\lambda} (-1)^{n-j}\frac{(j\lambda)^{n-j}}{(n-j)!}z^n -\sum_{n=1}^\infty\sum_{j=0}^{n-1} e^{j\lambda}(-1)^{n-j-1}\frac{(j\lambda)^{n-j-1}}{(n-j-1)!}z^n \right].
\end{align}
Now, $\pi_0 = 1-\lambda$ and $\pi_1 = (1-\lambda)(e^\lambda-1)$, so for $n\geqslant2$ we have
\begin{align}
\pi_n &= (1-\lambda)\sum_{j=0}^n e^{j\lambda}(-1)^{n-j}\frac{ j\lambda(j\lambda )^{n-j-1}}{(n-j)!} + \sum_{j=0}^n e^{j\lambda} (-1)^{n-j} \frac{(j\lambda)^{n-j-1}(n-j)}{(n-j)!}\\
&= (1-\lambda)\sum_{j=0}^n e^{j\lambda}(-1)^{n-j}\frac{(j\lambda+n-j)(j\lambda)^{n-j-1}}{(n-j)!}.
\end{align}
