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

Question

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$.

Answer 1

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.

Answer 2

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.

Answer 3

$\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.