Difference Equations Approach for Multi-Server Queueing Models with Removable Servers

Abstract

We consider an extended form of the M^X/M/c queue with two types of server groups: Static as well as dynamic (which turn on/off in a state-dependent manner) servers. The two server groups may have homogenous or non-homogenous service rates. The model is further extended to feature setup and delayed-off times, finite capacity, and k staffing levels. This class of queues is solved via the difference equations approach, which addresses narratives in the literature and achieves higher numerical efficiency than the direct method. While the model of this queueing system is not new, the methodology for solving it is. Comparisons between our model and classic queues are provided followed by concluding remarks, including a summary of key observations.

Appendices

Appendix 1: Proving the existence of roots

In this appendix we prove that (8) has r roots inside the unit circle |z| = 1. We first multiply both sides of the root equation (8) by z^r yielding

$$ {z}^r=\frac{1}{\lambda + c\mu +l{\mu}_1}\left[\lambda \sum \limits_{h=1}^r{b}_h{z}^{r-h}+\left( c\mu +l{\mu}_1\right){z}^{r+1}\right] $$

Let f(z) = z^r and $g(z)=-\frac{1}{\lambda + c\mu +l{\mu}_1}\left[\lambda \sum \limits_{h=1}^r{b}_h{z}^{r-h}+\left( c\mu +l{\mu}_1\right){z}^{r+1}\right]$ such that f(z) + g(z) = 0. Consider the magnitudes of f(z) and g(z) on the contour |z| = 1 − τ where τ is positive and sufficiently small. This gives

$$ \left|f(z)\right|={\left(1-\tau \right)}^r=1-\tau r+o\left(\tau \right) $$

and

$$ \left|g(z)\right|\le \frac{1}{\lambda + c\mu +l{\mu}_1}\left[\lambda \sum \limits_{h=1}^r{b}_h{\left|z\right|}^{r-h}+\left( c\mu +l{\mu}_1\right){\left|z\right|}^{r+1}\right] $$

Letting |z| = 1 − τ in the right-hand side of the above expression leads to the following:

$$ {\displaystyle \begin{array}{l}\left|g(z)\right|\le \frac{1}{\lambda + c\mu +l{\mu}_1}\left[\lambda \sum \limits_{h=1}^r{b}_h\left(1-\left(r-h\right)\tau \right)+\left( c\mu +l{\mu}_1\right)\left(1-\left(r+1\right)\tau \right)+o\left(\tau \right)\right]\\ {}\le \frac{1}{\lambda + c\mu +l{\mu}_1}\left[\lambda + c\mu +l{\mu}_1-\left(\lambda + c\mu +l{\mu}_1\right) r\tau +\lambda E\left[X\right]\tau -\left( c\mu +l{\mu}_1\right)\tau +o\left(\tau \right)\right]\end{array}} $$

Using the definition of ρ from Section 2, the above expression can be rearranged to give

$$ \left|g(z)\right|\le 1- r\tau -\frac{\left( c\mu +l{\mu}_1\right)\left(1-\rho \right)\tau }{\lambda + c\mu +l{\mu}_1}+o\left(\tau \right) $$

The fact that ρ < 1 implies that |g(z)| < |f(z)| on |z| = 1 − τ. Since f(z) and g(z) satisfy the conditions of Rouché’s theorem it follows that (8) has r roots inside the unit circle.

Appendix 2: Reduction of N _R

In this Section we demonstrate that (9) is also true for n ≤ i ≤ n + r − 1. The benefit of doing so is in the analytical reduction of N_R by r which subsequently enables even further reduction of N_R (the effect of such a reduction is demonstrated in Table 1). We begin this procedure by letting i = n + 2r − 1 in the balance equation (6) and expressing probabilities using (9) where applicable:

$$ \left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{n+2r-1}=\lambda \sum \limits_{h=1}^{r-1}{b}_h\sum \limits_{j=1}^r{C}_j{z}_j^{n+2r-1-h}+\lambda {b}_r{\pi}_{n+r-1,1}+\left( c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{n+2r} $$

This can be rearranged:

$$ \left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{n+2r-1}\left[1-\frac{1}{\lambda + c\mu +l{\mu}_1}\left\{\lambda \sum \limits_{h=1}^r{b}_h{z}_j^{-h}-\left( c\mu +l{\mu}_1\right){z}_j\right\}+\frac{\lambda {b}_r{z}_j^{-r}}{\lambda + c\mu +l{\mu}_1}\right]=\lambda {b}_r{\pi}_{n+r-1,1} $$

Applying (8) to the above expression and given that λ and b_r are both strictly positive, we have

$$ {\pi}_{n+r-1,1}=\sum \limits_{h=1}^r{C}_h{z}_h^{n+r-1} $$

By letting i = n + 2r − 2, n + 2r − 3, …, n + r + 1, n + r, we have the following result:

$$ {\pi}_{i,1}=\sum \limits_{h=1}^r{C}_h{z}_h^i,\left(n\le i\le n+r-1\right) $$

(27)

By deriving (27), we have reduced N_R by r, it went from 2n − m + 2r − 1 to 2n − m + r − 1.

1.1 Appendix 2.1: Further reduction of N _R

Further reduction of N_R is desired as it enables efficient numerical computations. To perform such a reduction we must distinguish and treat each of the following two cases separately: Case 1 occurs when r ≥ n and Case 2 occurs when r < n.

1.1.1 Appendix 2.1.1: Case 1: r ≥ n

In this case, an incoming batch of size h, (1 ≤ h ≤ r) could be equal to or larger than n so that the l dynamic servers are turned on immediately upon arrival of the batch. Using (27) we concluded that there are N_R = 2n − m + r − 1 unknowns: {π_i,0, 0 ≤ i ≤ n − 1}, {π_i,1, m + 1 ≤ i ≤ n − 1}, and C_h, (1 ≤ h ≤ r). To further reduce N_R we let i = n + r − 1 in balance equation (5) and express π_{i, 1} with (27) where applicable. Doing so gives

$$ \left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{n+r-1}=\lambda \left({b}_r{\pi}_{n-1,0}+\sum \limits_{h=1}^r{b}_h{\pi}_{n+r-1-h,1}\right)+\left( c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{n+r} $$

The above expression is then rearranged to yield

$$ \left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{n+r-1}\left[1-\frac{1}{\lambda + c\mu +l{\mu}_1}\left\{\lambda \sum \limits_{h=1}^r{b}_h{z}_j^{-h}-\left( c\mu +l{\mu}_1\right){z}_j\right\}+\frac{\lambda }{\lambda + c\mu +l{\mu}_1}{b}_r{z}_j^{-r}\right]=\lambda {b}_r\left({\pi}_{n-1,0}+{\pi}_{n-1,1}\right) $$

Applying (8) to the above expression and given that λ and b_r are both strictly positive, we have

$$ \sum \limits_{j=1}^r{C}_j{z}_j^{n-1}={\pi}_{n-1,0}+{\pi}_{n-1,1} $$

We let i = n + r − 2, n + r − 3, …, m + r + 2, m + r + 1 in balance equation (5) and prove that

$$ \sum \limits_{j=1}^r{C}_j{z}_j^i={\pi}_{i,0}+{\pi}_{i,1},\left(m+1\le i\le n-1\right) $$

(28)

We proceed further for the remaining values of i (i.e. i = m + r, m + r − 1, …, r + 1, r). Let i = m + r in balance equation (5) and express π_{i, 0} + π_{i, 1} with (28) where applicable. Doing so gives

$$ \left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{m+r}=\lambda \left(\sum \limits_{h=m+r-n+1}^r{b}_h{\pi}_{m+r-h,0}+\sum \limits_{h=1}^{r-1}{b}_h{\pi}_{m+r-h,1}\right)+\left( c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{m+r+1} $$

which can be rearranged to

$$ {\displaystyle \begin{array}{l}\left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{m+r}\\ {}=\lambda \left\{\sum \limits_{h=1}^{m+r-n}{b}_h{\pi}_{m+r-h,1}+\sum \limits_{h=m+r-n+1}^{r-1}{b}_h\left({\pi}_{m+r-h,0}+{\pi}_{m+r-h,1}\right)+{b}_r{\pi}_{m,0}\right\}+\left( c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{m+r+1}\end{array}} $$

or

$$ \left(\lambda + c\mu +l{\mu}_1\right)\sum \limits_{j=1}^r{C}_j{z}_j^{m+r}\left[1-\frac{1}{\lambda + c\mu +l{\mu}_1}\left\{\lambda \sum \limits_{h=1}^r{b}_h{z}_j^{-h}+\left( c\mu +l{\mu}_1\right){z}_j\right\}+\frac{b_r{z}_j^{-r}}{\lambda + c\mu +l{\mu}_1}\right]=\lambda {b}_r{\pi}_{m,0} $$

Applying (8) to the above expression and given that λ and b_r are both strictly positive, we have

$$ {\pi}_{m,0}=\sum \limits_{j=1}^r{C}_j{z}_j^m $$

By letting i = m + r − 1, m + r − 2, …, r + 1, r in balance equation (5), it can be shown that

$$ \sum \limits_{j=1}^r{C}_j{z}_j^i={\pi}_{i,0},\left(0\le i\le m\right) $$

(29)

Therefore when r ≥ n, by deriving expression (28) and (29) we have further reduced N_R by n so that it is reduced from 2n − m + r − 1 to n − m + r − 1. The needed N_R equations can be generated from the balance equations such that {π_{i, s}, i ≥ 0, s = 0, 1} can be explicitly expressed as

$$ {\pi}_{i,s}=\left\{\begin{array}{c}\sum \limits_{l=1}^r{C}_l{z}_l^i,\left(0\le i\le m,s=0\right)\kern10.25em \\ {} already\ determined,\left(m+1\le i\le n-1,s=0\right)\\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i-{\pi}_{i,0},\left(m+1\le i\le n-1,s=1\right)\kern3.5em \\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i,\left(i\ge n,s=1\right)\kern12.75em \end{array}\right. $$

(30)

where the ‘already determined’ probabilities are those that are simultaneously found along with the C_h’s when solving the N_R equations.

1.1.2 Appendix 2.1.2: Case 2: r < n

In this case, we assume that an incoming batch of size h, (1 ≤ h ≤ r) will prompt the l dynamic servers to turn on immediately upon arrival of the batch. With (27) found we have concluded that there are N_R = 2n − m + r − 1 unknowns: {π_i,0, 0 ≤ i ≤ n − 1}, {π_i,1, m + 1 ≤ i ≤ n − 1}, and C_h, (1 ≤ h ≤ r). In reducing N_R for Case 2 we must further consider two subcases: n − r ≤ m and n − r > m. As a remark, readers will later see that both of these subcases lead to the reduction of N_R by r. However, such a separation needs to be made as the expressions for π_i,s for each subcase are different.

1.1.3 Appendix 2.1.2.1: Subcase 1: n − r ≤ m

The procedure to compute {π_i,0, 0 ≤ i ≤ n − 1}, {π_i,1, m + 1 ≤ i ≤ n − 1}, and C_h, (1 ≤ h ≤ r) when n − r ≤ m follows the same procedure as provided in Appendix 2.1.1 up to the derivation of (28). However, after (28), instead of letting i = m + r, m + r − 1, …, r + 1, r in the balance equation (5), we let i = m + r, m + r − 1, …, n + 1, n as r < n. Doing so leads to

$$ \sum \limits_{j=1}^r{C}_j{z}_j^i={\pi}_{i,0},\left(n-r\le i\le m\right) $$

(31)

Therefore when n − r ≤ m, by deriving expression (31) we have further reduced N_R by r, from 2n − m + r − 1 to 2n − m − 1. The needed N_R equations can be generated from the balance equations such that {π_i,s, i ≥ 0, s = 0, 1} can be explicitly expressed as

$$ {\pi}_{i,s}=\left\{\begin{array}{c} already\ determined,\left(0\le i\le n-r-1,s=0\right)\ \\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i,\left(n-r\le i\le m,s=0\right)\kern8.5em \\ {} already\ determined,\left(m+1\le i\le n-1,s=0\right)\\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i-{\pi}_{i,0},\left(m+1\le i\le n-1,s=1\right)\kern3.5em \\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i,\left(i\ge n,s=1\right)\kern12.75em \end{array}\right. $$

(32)

where the ‘already determined’ probabilities are those that are simultaneously found along with the C_h’s when solving the N_R equations.

1.1.4 Appendix 2.1.2.2: Subcase 2: n − r > m

The procedure to compute {π_i,0, 0 ≤ i ≤ n − 1}, {π_i,1, m + 1 ≤ i ≤ n − 1}, and C_h, (1 ≤ h ≤ r) when n − r > m is slightly different than the procedure provided in Appendix 2.1.2.1. Instead of letting i = n + r − 1, n + r − 2, …, m + r + 2, m + r + 1 in the balance equation (5) in Appendix 2.1.2.1, we let i = n + r − 1, n + r − 2, …, n + 1, n as n − r > m. Doing so leads to

$$ \sum \limits_{j=1}^r{C}_j{z}_j^i={\pi}_{i,0}+{\pi}_{i,1},\left(n-r\le i\le n-1\right) $$

(33)

Therefore when n − r > m, by deriving expression (33) we have further reduced N_R by r (as done in Appendix 2.1.2.1). The needed N_R equations can be generated from the balance equations such that {π_{i, s}, i ≥ 0, s = 0, 1} can be explicitly expressed as

$$ {\pi}_{i,s}=\left\{\begin{array}{c} already\ determined,\left(0\le i\le n-r-1,s=0,1\right)\\ {} already\ determined,\left(n-r\le i\le n-1,s=0\right)\kern0.75em \\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i-{\pi}_{i,0},\left(n-r\le i\le n-1,s=1\right)\kern4.25em \\ {}\sum \limits_{l=1}^r{C}_l{z}_l^i,\left(i\ge n,s=1\right)\kern13.25em \end{array}\right. $$

(34)

where the ‘already determined’ probabilities are those that are simultaneously found along with the C_h’s when solving the N_R equations.

Appendix 3: Balance equations for the extension to k staffing levels

The transition dynamics of the M^X/M/c + l/(m, n) queue with k staffing levels are provided. While the balance equations (1) and (2) from the baseline model remain unchanged, the rest of the balance equations are modified to the following:

$$ \left(\lambda + c\mu \right){\pi}_{i,0}=\lambda \sum \limits_{h=1}^{\mathit{\min}\left(i,r\right)}{b}_h{\pi}_{i-h,0}+ c\mu {\pi}_{i+1,0},\left(c\le i\le {m}_1-1\right) $$

(35)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^{s-1}{l}_j{\mu}_j\right){\pi}_{m_s,s-1}=\lambda \sum \limits_{h={m}_s-{n}_1+1}^{\min \left(r,{m}_s\right)}{b}_h{\pi}_{m_s-h,0}+\lambda \sum \limits_{j=1}^{s-2}\sum \limits_{h={m}_s-{n}_{j+1}+1}^{\min \left(r,{m}_s-{m}_j-1\right)}{b}_h{\pi}_{m_s-h,j}+\lambda \sum \limits_{h=1}^{\min \left(r,{m}_s-{m}_{s-1}-1\right)}{b}_h{\pi}_{m_s-h,s-1}+\left( c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{m_s+1,s}+\left( c\mu +\sum \limits_{j=1}^{s-1}{l}_j{\mu}_j\right){\pi}_{m_s+1,s-1},\left(1\le s\le k\right) $$

(36)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^{s-1}{l}_j{\mu}_j\right){\pi}_{i,s-1}=\lambda \sum \limits_{h=i-{n}_1+1}^{\min \left(r,i\right)}{b}_h{\pi}_{i-h,0}+\lambda \sum \limits_{j=1}^{s-2}\sum \limits_{h=i-{n}_{j+1}+1}^{\min \left(r,i-{m}_j-1\right)}{b}_h{\pi}_{i-h,j}+\lambda \sum \limits_{h=1}^{\min \left(r,i-{m}_{s-1}-1\right)}{b}_h{\pi}_{i-h,s-1}+\left( c\mu +\sum \limits_{j=1}^{s-1}{l}_j{\mu}_j\right){\pi}_{i+1,s-1},\left({m}_s+1\le i\le {n}_s-2,1\le s\le k\right) $$

(37)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^{s-1}{l}_j{\mu}_j\right){\pi}_{n_s-1,s-1}=\lambda \sum \limits_{h={n}_s-{n}_1}^{\min \left(r,{n}_s-1\right)}{b}_h{\pi}_{n_s-1-h,0}+\lambda \sum \limits_{j=1}^{s-2}\sum \limits_{h={n}_s-{n}_{j+1}}^{\min \left(r,{n}_s-{m}_j-2\right)}{b}_h{\pi}_{n_s-1-h,j}+\lambda \sum \limits_{h=1}^{\min \left(r,{n}_s-{m}_{s-1}-2\right)}{b}_h{\pi}_{n_s-1-h,s-1},\left(1\le s\le k\right) $$

(38)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{m_s+1,s}=\left( c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{m_s+2,s},\left(1\le s\le k\right) $$

(39)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{i,s}=\lambda \sum \limits_{h=1}^{\min \left(r,i-{m}_s-1\right)}{b}_h{\pi}_{i-h,s}+\left( c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{i+1,s},\left({m}_s+2\le i\le {n}_s-1,1\le s\le k\right) $$

(40)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{i,s}=\lambda \sum \limits_{h=i-{n}_1+1}^{\min \left(r,i\right)}{b}_h{\pi}_{i-h,0}+\lambda \sum \limits_{j=1}^{s-1}\sum \limits_{h=i-{n}_{j+1}+1}^{\min \left(r,i-{m}_j-1\right)}{b}_h{\pi}_{i-h,j}+\lambda \sum \limits_{h=1}^{\min \left(r,i-{m}_s-1\right)}{b}_h{\pi}_{i-h,s}+\left( c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{i+1,s},\left({n}_s\le i\le {n}_s+r-1,1\le s\le k\right) $$

(41)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{i,s}=\lambda \sum \limits_{h=1}^{\min \left(r,i-{n}_s\right)}{b}_h{\pi}_{i-h,s}+\left( c\mu +\sum \limits_{j=1}^s{l}_j{\mu}_j\right){\pi}_{i+1,s},\left({n}_s+r\le i\le {m}_{s+1}-1,1\le s\le k-1\right) $$

(42)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^k{l}_j{\mu}_j\right){\pi}_{i,k}=\lambda \sum \limits_{h=1}^{\min \left(r,i-{n}_k\right)}{b}_h{\pi}_{i-h,k}+\left( c\mu +\sum \limits_{j=1}^k{l}_j{\mu}_j\right){\pi}_{i+1,k},\left({n}_k+r\le i\le {n}_k+2r-1\right) $$

(43)

$$ \left(\lambda + c\mu +\sum \limits_{j=1}^k{l}_j{\mu}_j\right){\pi}_{i,k}=\lambda \sum \limits_{h=1}^{\min \left(r,i-{n}_k-r\right)}{b}_h{\pi}_{i-h,k}+\left( c\mu +\sum \limits_{j=1}^k{l}_j{\mu}_j\right){\pi}_{i+1,k},\left(i\ge {n}_k+2r\right) $$

(44)

Appendix 4: Extension to the M ^X/M/c + l/(m, n)/K queue

The baseline model can be extended to feature a finite capacity such that the total number of jobs held by the system is finite. Therefore, the M^X/M/c + l/(m, n) queue with finite capacity can house up to K, (1 ≤ K <  + ∞) jobs in the system where K includes the jobs in queue as well as those being served by both the static and dynamic servers (if any). Therefore we have the M^X/M/c + l/(m, n)/K queue with the joint steady-state distribution {π_{i, s}, 0 ≤ i ≤ K, s = 0, 1} and the normalizing condition

$$ \sum \limits_{i=0}^{n-1}{\pi}_{i,0}+\sum \limits_{i=m+1}^K{\pi}_{i,1}=1 $$

(45)

With the introduction of finite capacity, an incoming batch can be rejected if its size (h) exceeds the available space (K − i). When h > K − i the model M^X/M/c + l/(m, n)/K is subject to one of the following two rejection policies: partial rejection of a batch occurs when out of h jobs the K − i jobs are admitted into the system and the remaining (h − K + i) jobs are rejected. Total rejection of a batch occurs when, given the same condition, the entire batch is rejected. The balance equations that describe the system dynamics of the M^X/M/c + l/(m, n)/K queue can be derived by modifying the balance equation (7) from Section 2.1 to incorporate each rejection policy. These are provided in the following two sections.

1.1 Appendix 4.1: M ^X/M/c + l/(m, n)/K queue with partial rejection

$$ \left(\lambda + c\mu +l{\mu}_1\right){\pi}_{i,1}=\lambda \sum \limits_{h=1}^{\mathit{\min}\left(i-n-r,r\right)}{b}_h{\pi}_{i-h,1}+\left( c\mu +l{\mu}_1\right){\pi}_{i+1,1},\left(n+2r\le i\le K-1\right) $$

(46)

$$ \left( c\mu +l{\mu}_1\right){\pi}_{K,1}=\lambda \sum \limits_{j=1}^r\sum \limits_{h=j}^r{b}_h{\pi}_{K-j,1} $$

(47)

1.2 Appendix 4.2: M ^X/M/c + l/(m, n)/K queue with total rejection

$$ \left(\lambda + c\mu +l{\mu}_1\right){\pi}_{i,1}=\lambda \sum \limits_{h=1}^{\mathit{\min}\left(i-n-r,r\right)}{b}_h{\pi}_{i-h,1}+\left( c\mu +l{\mu}_1\right){\pi}_{i+1,1},\left(n+2r\le i\le K-r-1\right) $$

(48)

$$ \left(\lambda \sum \limits_{h=1}^{K-i}{b}_h+ c\mu +l{\mu}_1\right){\pi}_{i,1}=\lambda \sum \limits_{h=1}^{\mathit{\min}\left(i-n-r,r\right)}{b}_h{\pi}_{i-h,1}+\left( c\mu +l{\mu}_1\right){\pi}_{i+1,1},\left(K-r\le i\le K-1\right) $$

(49)

$$ \left( c\mu +l{\mu}_1\right){\pi}_{K,1}=\lambda \sum \limits_{h=1}^r{b}_h{\pi}_{K-h,1} $$

(50)

The above balance equations can be solved via the difference equations approach as demonstrated in earlier sections of this paper.

Appendix 5: Properties of difference equations

The difference equations approach we introduced in solving the baseline model and its extensions is based on interpreting the model’s balance equations as difference equations. By doing so, we can express the solution in terms of roots by leveraging the well-established properties of linear difference equations. As discussed in Chaudhry and Templeton (1983), an equation of the type

$$ {a}_0{f}_{x+n}+{a}_1{f}_{x+n-1}+\dots +{a}_{n-1}{f}_{x+1}+{a}_n{f}_x={b}_x,\left(x=1,2,\dots \right) $$

where the a_i are known constants, f_i are unknown functions to be determined, and b_x is a given function of x, is called a nonhomogeneous linear difference equation of order n. If b_x = 0, for all x, then it is called a homogenous linear difference equation with constant coefficients. A general solution to the above nonhomogeneous equation consists of two parts:

• 1. A linear combination of all solutions to the homogeneous equation; and

• 2. A particular solution to the nonhomogeneous equation.

The solution to the homogeneous part of the equation proceeds along the following lines. Letting f_x = Cz^x in the homogeneous equation leads to

$$ {a}_0C{z}^{x+n}+{a}_1C{z}^{x+n-1}+\dots +{a}_{n-1}C{z}^{x+1}+{a}_nC{z}^x=0 $$

and

$$ {a}_0{z}^n+{a}_1{z}^{n-1}+\dots +{a}_{n-1}z+{a}_n=0 $$

The last equation in z, being an n-th degree equation, gives n roots (real or complex, distinct or coincident). As a consequence, assuming that the roots are distinct, the general solution of the homogeneous part is written as

$$ {f}_x=\sum \limits_{j=1}^n{C}_j{z}_j^x $$