# Is there a general method for solving this type of recurrence?

Edit: Here is the original problem; it is possible that my recurrence for the stationary distribution $\pi$ is incorrect.

Consider a single server queue where customers arrive according to a Poisson process with intensity $\lambda$ and request i.i.d. $\mathsf{Exp}(\mu)$ service times. The server is subject to failures and repairs. The lifetime of a working server is an $\mathsf{Exp}(\theta)$ random variable, while the repair time is an $\mathsf{Exp}(\alpha)$ random variable. Successive lifetimes and repair times are independent, and are independent of the number of customers in the queue. When the server fails, all the customers in the queue are forced to leave, and while the server is under repair no new customers are allowed to join.

Edit: I have revised the recurrence.

In a problem on queueing theory I've derived the following recurrence: \begin{align} \pi_1 &=\left(\frac{\lambda+\theta}\mu\right)\pi_0 - \frac{\alpha\theta}{\mu(\alpha+\theta)}\\ \pi_{n+1} &= \left(1+\frac{\lambda+\theta}\mu\right)\pi_n - \frac\lambda\mu\pi_{n-1},\ n\geqslant1. \end{align} where $\lambda$, $\mu$, $\theta$, and $\alpha$ are positive constants and $$\sum_{i=0}^\infty \pi_i = \frac\alpha{\alpha+\theta}.$$

After a lot of tedious algebra, I found that $$\scriptsize\pi_n = \left(\frac{\alpha \theta \left(2 \theta (\lambda +\mu )+(\lambda -\mu )^2\right) \left(\theta +\lambda +\mu+ \sqrt{\theta ^2+2 \theta (\lambda +\mu )+(\lambda -\mu )^2}\right)^n}{(\alpha +\theta ) (2 \mu )^n \sqrt{\theta ^2+2 \theta (\lambda +\mu )+(\lambda -\mu )^2}}\right)(1+\pi_0)$$ for $n\geqslant 1$. To save space, let $$\mathcal C:=\sqrt{\theta ^2+2 \theta (\lambda +\mu )+(\lambda -\mu )^2}.$$

Summing over $n$ and solving for $\pi_0$, I found $$\pi_0 =\frac{\alpha \mu \left(2 \theta (\lambda +\mu )+(\lambda -\mu )^2\right) \left(\lambda -\mu-\theta-\mathcal C \right)}{2 \theta (\alpha +\theta ) \mathcal C},$$

and so $$\pi_n=\left(\frac{ \left(2 \theta (\lambda +\mu )+(\lambda -\mu )^2\right) \left(\lambda -\mu-\theta-\mathcal C \right)+2 \theta (\alpha +\theta ) \mathcal C }{2(\alpha +\theta )^2\mathcal C^2\left(\alpha^2\mu \left(2 \theta (\lambda +\mu )+(\lambda -\mu )^2\right)\right)^{-1}} \right)\left(\frac{\lambda+\mu+\theta+\mathcal C }{2\mu}\right)^n.$$

If you see any errors let me know...

I'm also wondering what conditions on $\lambda,\mu,\theta$, and $\alpha$ are necessary for $\sum_{i=0}^\infty \pi_i$ to converge. For context, this is a $M/M/1$ queue with arrival rate $\lambda$, service rate $\mu$, but with an added state $D$ with transitions of rate $\theta$ from each state $n$ to $D$ and a transition of rate $\alpha$ from $D$ to $0$.

• Pushing all $\pi_k$ terms to the LHS, you've got a denumerably-infinite system of linear equations. Perhaps it'd be useful to view it as $A\Pi = B$ where $A,\Pi,B$ are an infinite matrix and infinite column vectors respectively. – Semiclassical May 29 '16 at 1:16
• The typo was actually on my part---mea culpa! But the additional context is helpful. – Semiclassical May 29 '16 at 4:27
• I note that you've changed the summation to have upper limit $i=n$ instead of $i=n-1$. If this was intentional, then the answers both I and @FelixMarin provided are no longer quite right (Mine, for instance, will have $\frac{\theta}{\mu}x P(x)$ instead of $\frac{\theta}{\lambda}x P(x)$.) – Semiclassical May 29 '16 at 14:17
• Yes, that is because we defined $\varphi_n:=\sum_{i=0}^n \pi_i$. – Math1000 May 29 '16 at 14:19
• Can you explain how you derived the recurrence? It seems that in modelling the problem using a Markov chain, we would need to have a special state for when the server is out of repair, distinct from the ordinary state with an empty queue. How is this reflected in your $\pi_n$ notation? I am assuming $\pi_n$ is supposed to represent the probability of being in state $n$, under the equilibrium distribution. – Brent Kerby May 29 '16 at 14:21

Hint: By defining $\phi_n = \sum_{i=0}^n \pi_i$, we may express $\pi_{n+1}$ and $\phi_{n+1}$ as linear combinations of $\pi_n$ and $\phi_n$ plus constants, giving a first-order linear matrix difference equation.

EDIT: The equation we end up with is of the form $$x_{(n+1)} = Ax_{(n)}+b$$ where the vectors $x_{(n)}=\begin{pmatrix}\pi_n \\ \phi_n\end{pmatrix}\in\mathbb R^2$ are to be solved for, and $A$ is a known $2\times 2$ matrix and $b\in\mathbb R^2$ a known vector. If $A$ can be diagonalized, then, writing $A=SDS^{-1}$ and setting $y_{(n)}=S^{-1}x_{(n)}$, the system becomes $$y_{(n+1)} = Dy_{(n)}+S^{-1}b$$ In this case, the system decomposes into univariate difference equations, which can be solved separately, and then one uses $x_{(n)} = Sy_{(n)}$ to solve for the original unknowns.

• This is an impressively simple approach. Well done! Though, should the definition for $\phi_n$ start at $i=1$ or $i=0$? – Semiclassical May 29 '16 at 4:41
• Good catch! I just edited it to fix that. – Brent Kerby May 29 '16 at 4:56
• @BrentKerby I'm not familiar with multivariate recurrences. Could you elaborate on how to solve for $\pi_n,\varphi_n$? – Math1000 May 29 '16 at 6:41
• Actually, I'm a little surprised that the summation ends at $i=n$ rather than $i=n-1$. In the former case, $\phi_n$ depends on $\pi_n$ whereas in the latter you can write $\phi_n=\phi_{n-1}+\pi_n$. – Semiclassical May 29 '16 at 14:33
• @Math1000, I've added an explanation. I may have caused some confusion by calling it a multivariate difference equation, a term which is used to describe a discrete analogue of a PDE; I should have called it a matrix difference equation. – Brent Kerby May 29 '16 at 14:33

$\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{\mathrm{i}} \newcommand{\iff}{\Leftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\, #2 \,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ The question: $\ds{\pi_{n+1} = \overbrace{\frac{\lambda+\mu\theta}{\mu}}^{\ds{\equiv a}}\ \pi_n\ +\ \overbrace{\frac\theta\mu}^{\ds{\equiv b}}\sum_{i=0}^{n-1}\pi_i\ -\ \overbrace{\frac{\alpha\theta}{\mu(\alpha+\theta)}}^{\ds{\equiv c}} \quad\imp\quad\pi_{n+1} = a\pi_n + b\sum_{i=0}^{n-1}\pi_i - c\quad}$ and $\ds{\sum_{i = 0}^{\infty}\pi_{i}\ =\ \underbrace{{\alpha \over \alpha + \theta}}_{\ds{\equiv\ d}}\ =\ {c \over b}}$

With $\verts{z} < 1$: \begin{align} \sum_{n = 0}^{\infty}\pi_{n + 1}z^{n} & = a\sum_{n = 0}^{\infty}\pi_nz^{n} + b\sum_{n = 0}^{\infty}z^{n}\sum_{i=0}^{n - 1}\pi_i - c\sum_{n = 0}^{\infty}z^{n} \\[3mm]\imp {1 \over z}\pars{\sum_{n = 0}^{\infty}\pi_{n}z^{n} - \pi_{0}} & = a\sum_{n = 0}^{\infty}\pi_nz^{n} + b\sum_{i = 0}^{\infty}\pi_{i}\sum_{n = 1 + i}^{\infty}z^{n} - {c \over 1 - z} \\[3mm]\imp \pars{{1 \over z} - a}\sum_{n = 0}^{\infty}\pi_{n}z^{n} & = {\pi_{0} \over z} + b\sum_{i = 0}^{\infty}\pi_{i}{z^{i + 1} \over 1 - z} - {c \over 1 - z} \\[3mm]\imp \pars{{1 \over z} - a - b\,{z \over 1 - z}}\sum_{n = 0}^{\infty}\pi_{n}z^{n} & = {\pi_{0} \over z} - {c \over 1 - z} \end{align}

Then, $$\sum_{i = 0}^{\infty}\pi_{i}z^{i} = {\pars{\pi_{0} + c}z - \pi_{0} \over \pars{b - a}z^{2} + \pars{a + 1}z - 1}$$

In order to get the set $\braces{\pi_{i}}$, expand th right hand side in powers of $z$. Maybe, some other conditions on the magnitud of $z$ will be required along the way.

Note: The present form of the answer reflects an earlier version of the problem. I plan to modify it to reflect the changes, hopefully sooner rather than later...

Let $P(x)=\sum_{n=0}^\infty \pi_n x^n$. The first few lines of this recurrence are \begin{align} \pi_1+\frac{1}{\mu}\left(\frac{\alpha\theta}{\alpha+\theta}\right) &= \left(\frac{\lambda+\mu\theta}{\mu}\right)\pi_0,\\ \pi_2+\frac{1}{\mu}\left(\frac{\alpha\theta}{\alpha+\theta}\right) &= \left(\frac{\lambda+\mu\theta}{\mu}\right)\pi_1+\frac{\theta}{\mu}\pi_0,\\ \pi_3+\frac{1}{\mu}\left(\frac{\alpha\theta}{\alpha+\theta}\right) &= \left(\frac{\lambda+\mu\theta}{\mu}\right)\pi_2+\frac{\theta}{\mu}\pi_1+\frac{\theta}{\mu}\pi_0,\\ \end{align} and so on. Multiplying each line by powers of $x$ (starting with $x^1$) and summing yields

$$P(x)-\pi_0 +\frac{1}{\mu}\left(\frac{\alpha\theta}{\alpha+\theta}\right)(x+x^2+x^3+\cdots) = \left(\frac{\lambda+\mu\theta}{\mu}\right)xP(x)+\frac{\theta}{\lambda}\left(x^2+x^3+\cdots\right)P(x).$$ We clean this up by multiplying both sides by $(1-x)$, yielding $$(1-x)(P(x)-\pi_0)+\frac{1}{\mu}\left(\frac{\alpha\theta}{\alpha+\theta}\right)x=\left(\frac{\lambda+\mu\theta}{\mu}\right)x(1-x)P(x)+\frac{\theta}{\mu}x^2 P(x).$$ As a check, for $x=1$ this implies $$\frac{1}{\mu}\left(\frac{\alpha\theta}{\alpha+\theta}\right)=\frac{\theta}{\mu} P(1)\implies P(1)=\sum_{n=0}^\infty \pi_n =\frac{\alpha}{\alpha+\theta}$$ which is the desired normalization condition. What remains is to solve for $P(x)$ and then expand to obtain the coefficients $\{\pi_n\}$, a task I leave to the interested reader.

• I find that $$P(x) = \frac{\pi _0 \mu (\alpha +\theta )-x \left(\alpha \left(\theta +\pi _0 \mu \right)+\pi _0 \theta \mu \right)}{(\alpha +\theta ) \left(\mu +x^2 (\theta (\mu -1)+\lambda )-x (\theta \mu +\lambda +\mu )\right)}.$$ Needless to say, $\texttt{SeriesCoefficient}$ in Mathematica yields a truly horrifying result :( – Math1000 May 29 '16 at 7:03
• What I find more bothersome than the horrible-looking result is that it seems underdetermined, since the higher coefficients all depend on $\pi_0$. Is there a condition that's missing? @Math1000 – Semiclassical May 29 '16 at 13:48
• Here $\pi$ is a probability distribution, so $\pi_D + \sum_{i=0}^\infty \pi_i=1$. So it suffices to solve for $\pi_n$ in terms of $\pi_0$. – Math1000 May 29 '16 at 13:52