# Queueing Theory - 3.2 Single-Server Queues (M/M/1)

contributed by <[`kaeteyaruyo`](https://github.com/kaeteyaruyo)>

###### tags: `Queueing Theory`

> Figure 3.3 Rate transition diagram for the M/M/1 queue.

* M/M/1 是一種規律的排隊模型，具有以下性質
* M: Interarrival time 呈現指數分佈，到達率是定值 $\lambda$，其 PDF 為 $a(t) = \lambda e^{-\lambda t}$
* M: Service time 呈現指數分佈，服務率是定值 $\mu$，其 PDF 為 $b(t) = \mu e^{-\mu t}$
* 1: Server 數只有 1 個
* Fig. 3.3 是 M/M/1 的 Rate transition diagram，跟 birth-death process 的圖長的很像，只差在**到達率和服務率是不隨系統人數變化的定值**
* 每一個狀態代表**系統當中**（不是隊伍中）尚存有的人數
* 顧客到達可以視為一個 birth 事件，且 $\lambda_n = \lambda$
* 顧客離開可以視為一個 death 事件，且 $\mu_n = \mu$
* M/M/1 模型的 global balance equation 為
$$
\begin{cases}
(\lambda + \mu) p_n = \mu p_{n+1} + \lambda p_{n-1} \quad (n \geq 1) \\
\lambda p_0 = \mu p_1
\end{cases} \tag{3.6}
$$
其中 steady state 的 $p_n$ 的通式可以用三種方法來解，分別是
* Iterative method
* Generating funtion
* Operator D

## 3.2.1 Solving for $\{p_n\}$ Using an Iterative Method

* 直接套入 [(3.3) 式](https://hackmd.io/qI9i83Z3TxuReILYuCQfQQ) 的結論，會得知
$$
p_n = p_0 \cdot \prod_{i=1}^{n} \frac{\lambda}{\mu} = p_0 \Big(\frac{\lambda}{\mu}\Big)^{n}
$$
* 加上 Boundary condition $\sum_{\forall n} p_n = 1$，可以算出 $p_0$ 就是
$$
\begin{align}
& 1 = \sum_{n=0}^{\infty} p_n = \sum_{n=0}^{\infty} p_0 \Big(\frac{\lambda}{\mu}\Big)^{n} = p_0 \sum_{n=0}^{\infty} \rho^{n} \\
& \Rightarrow p_0 = \Big(\sum_{n=0}^{\infty} \rho^n\Big)^{-1}
\end{align}
$$
其中 $\rho = \frac{\lambda}{\mu}$，如同我們在 Chapter 1.4 時定義的，代表 traffic intensity 或是 utilization rate
* 上式中 $\sum_{n=0}^{\infty} \rho^n$ 的部份是一個公比是 $\rho$ 的無窮等比級數和，所以可以直接代入無窮等比級數和的公式，求得
$$
\sum_{n=0}^{\infty} \rho^n = \frac{1}{1 - \rho} \quad (\rho < 1) \\
\Rightarrow p_0 = 1 - \rho
$$
* $\rho$ 可以簡單想成這個 M/M/1 系統**忙碌的程度**。例如：若 $\rho = 0.7$，這代表所有的時間當中，大概有七成的時間這個系統中的 server 是處於忙碌的狀態
* 若 $\rho < 1$，代表這個系統可以穩定運作。**只有在這樣的情況， M/M/1 系統才會有穩定態**（另一個解釋：上面的等比級數和只有在公比小於 1 的時候才會收斂）
* 若 $\rho = 1$，代表這個系統已經滿載，此時系統會處於一個穩定增長的狀態，因為服務速度無法大於顧客抵達的速度，所以每當一個新的使用者進來，server 就愈難將這個新的顧客消化掉。這樣的系統不會收斂，沒有穩定態
* 若 $\rho > 1$，代表系統已經超過負荷，此時隊伍會無限制的增長，這樣的系統不會收斂，沒有穩定態
* 將 $p_0$ 代回 $(3.3)$ 式，即可獲得 $p_n$ 的通式
$$
p_n = \rho^np_0 = (1 - \rho) \rho^n \quad (\rho < 1)
$$

## 3.2.2 Solving for $\{p_n\}$ Using Generating Functions

* 定義**生成函數 (Generating funtion)** 為一個參數為複數的函數，其格式如
$$
P(z) = \sum_{n=0}^{\infty} z^n p_n, \quad |z| \leq 1
$$
* 將 $p_n$ 各項乘上 $z^n$ （他是第幾項就乘上 $z$ 的幾次方）後加總，就是生成函數
* 其參數 $z$ 只能在複數平面上的單位圓內，如此一來這個函數才會收斂
* 有時候直接找 $p_n$ 不好找，我們就先把他包裝成生成函數，然後試圖去**尋找這個函數的每個 $z^n$ 項前方的係數**，就可以得知 $p_n$ 的值了（只在乎係數，不在乎函數值）
* 先將 $(3.6)$ 式除上 $\mu$，化簡等式中的符號
$$
\begin{cases}
(\rho + 1) p_n = \rho p_{n-1} + p_{n+1} \quad (n \geq 1) \\
\rho p_0 = p_1
\end{cases}
$$
* 然後將兩式左右兩邊都乘上 $z^n$ 後加總，得到
$$
(\rho + 1) \sum_{n=1}^{\infty} z^n p_n + \rho p_0 = \rho \sum_{n=1}^{\infty} z^n p_{n-1} + \sum_{n=0}^{\infty} z^n p_{n+1}
$$
注意最後一項的加總下標是從 0 開始，因為下式等號右手邊的 $p_1$ 被歸在裏面（因為他沒有乘上 $\rho$）
* 我們想要將上式改以生成函數的多項式來表示，因此要儘量把 $\sum z^n p_n$ 的部份都用 $P(z)$ 代換掉
* 觀察生成函數的定義式，有兩點要注意
* $z^n$ 的次方數要和 $p_n$ 的底標一致
* 加總的下標**從 0 開始**（第一項一定要是 $z^0 p_0$）
* 如果式子中的級數和不滿足上述兩個條件，就要把他喬成正確的格式，才可以代換
* 有了以上概念，就可以魔改上面那條式子變成
$$
\begin{align}
& (\rho + 1)\sum_{n=0}^{\infty} z^n p_n - p_0 = \rho z \sum_{n=1}^{\infty} z^{n-1} p_{n-1} + \frac{1}{z} \sum_{n=0}^{\infty} z^{n+1} p_{n+1} \\
\Rightarrow & (\rho + 1)P(z) - p_0 = \rho z P(z) + \frac{1}{z} \big[P(z) - p_0\big]
\end{align}
$$
* 經過移項整理和因式分解後，可以得到 $P(z)$ 為
$$
[1 - \rho z - z + \rho z^2]P(z) = (1 - z) p_0 \\
\Rightarrow P(z) = \frac{p_0}{1 - \rho z}
$$
但 $p_0$ 還不知道是多少，所以還不算是通式
* 加上 boundary condition $\sum_{n=0}^{\infty} p_n = 1$ 可以求出 $p_0$。剛好把生成函數的 $z$ 代入 1，就是在求 $p_n$ 的加總，因此我們可以求出
$$
P(1) = \sum_{n=0}^{\infty} p_n 1^n = \sum_{n=0}^{\infty} p_n = 1 = \frac{p_0}{1 - \rho} \\
\Rightarrow p_0 = 1 - \rho, P(z) = \frac{1-\rho}{1-\rho z} \quad (\rho < 1, |z| \leq 1)
$$
* 但事實上我們要的不是 $P(z)$ 的通式，我們要的式 $P(z) = \sum_{n=0}^{\infty} z^n p_n$ 的 $p_n$。也就是說，我們必須想辦法把上面這個 $P(z)$ 的通式改寫成 $z^n$ 的級數和，然後觀察他們的係數才是答案
* 如果我們把上式中分子的 $1 - \rho$ 遮起來不看，會發現等式右手邊會是一個公比是 $\rho z$ 的等比級數和，因此我們可以把 $P(z)$ 改寫成
$$
\begin{align}
P(z) &= (1 - \rho) \Big(\frac{1}{1 - \rho z}\Big) \\
&= (1 - \rho) (1 + \rho z + (\rho z)^2 + (\rho z)^3 + ...) \\
&= \sum_{n=0}^{\infty} (1 - \rho) \rho^n z^n
\end{align}
$$
* 改寫至此，會發現這樣的格式已經吻合了我們一開始定義的生成函數的格式，因此可以觀察 $z^n$ 的係數得知
$$
p_n = (1 - \rho) \rho^n \quad (\rho = \frac{\lambda}{\mu} < 1)
$$
* 對於 M/M/1 這種很規律的模型來說，用這個方法實在是複雜到有點白痴。但是以後我們可能會遇到 $p_n$ 不那麼規律的模型，到時候就換成 iterative method 很難用了，所以才需要先介紹這種方法

## 3.2.3 Solving for $\{p_n\}$ Using Operators

* [差分方程式](https://zh.wikipedia.org/wiki/%E9%81%9E%E8%BF%B4%E9%97%9C%E4%BF%82%E5%BC%8F)是一種透過遞迴方式定義一個數列的方程式，數列的一般項 $a_n$ 被表示成他的前幾項的關係式
* 差分方程式的**階數**代表遞迴關係式中，和 $a_n$ 差最遠的那一項距離 $a_n$ 有多遠。例如：二階差分方程式就是指 $a_n$ 可以用 $a_{n-1}$ 和 $a_{n-2}$ 的關係式描述
* 若 $a_n$ 被定義成他前幾項的**線性組合**，那我們就稱這樣的差分方程式為線性差分方程式
* 依照以上定義，一個線性二階差分方程式的通式可以寫成
$$
a_n = Aa_{n-1} + Ba_{n-2}
$$
* 針對一個數列 $\{a_n\}$，我們可以定義一個 operator $D$，當他作用在數列某一項時，會得到他的下一項，也就是
$$
Da_n \equiv a_{n+1} \quad \text{(for all $n$)}
$$
* 對數列的某一項重複使用 operator $D$ $m$ 次，會得到他的後面第 $m$ 項，也就是
$$
D^ma_n \equiv a_{n+m} \quad \text{(for all $n$ and $m$)}
$$
* 因此，一個線性差分方程式就可以改以 operator $D$ 的方程式來表示。以二階線性方程式來舉例，他可以被改寫如下
$$
\begin{align}
& C_2 a_{n+2} + C_1 a_{n+1} + C_0 a_n = 0 \\
\Rightarrow & C_2 D^2 a_n + C_1 D a_n + C_0 a_n = 0 \\
\Rightarrow & (C_2 D^2 + C_1 D + C_0) a_n = 0 \tag{3.19}
\end{align}
$$
其中 $C_n$ 是各項的係數
* $(3.19)$ 式中， $a_n$ 的係數是一個 $D$ 的二次式，此二次式等於 0 的方程式即為該差分方程式的**特徵多項式**，也就是
$$
(C_2 D^2 + C_1 D + C_0) = 0
$$
* 此特徵多項式的根可以決定 $\{a_n\}$ 的通式。也就是不使用遞迴定義，可以直接用多項式的根來表達 $a_n$，其通式為 $a_n = \sum_{i=1}^k d_i r_i^n$。以上述二次式來說，通常會有兩個相異的解 $r_1$ 和 $r_2$，則 $a_n$ 就可以表示為
$$
a_n = d_1 r_1^n + d_2 r_2^n
$$
其中 $d_1$ 和 $d_2$ 是任意的常數
* 觀察 $p_n$ 在 $n \geq 1$ 的通式，將其移項整理並進行 index shifting 後，會得到
$$
\mu p_{n+2} = (\lambda + \mu) p_{n+1} - \lambda p_n, \; n \geq 0
$$
正好形如一個線性二階差分方程式。因此我們可以使用解差分方程式的手法來求出 $p_n$ 的通解
* 利用 operator $D$ 來求解 $p_n$ 的通解
$$
\begin{align}
& \mu p_{n+2} = (\lambda + \mu) p_{n+1} - \lambda p_n, \; n \geq 0 \\
\Rightarrow & \mu D^2 p_n - (\lambda + \mu) D p_n + \lambda p_n = 0, \; n \geq 0 \\
\Rightarrow & [\mu D^2 - (\lambda + \mu) D + \lambda] p_n = 0, \; n \geq 0 \\
\Rightarrow & (D - 1)(\mu D - \lambda) p_n = 0 \\
\Rightarrow & p_n = d_1 (1)^n + d_2 \big(\frac{\lambda}{\mu}\big)^n = d_1 + d_2 \rho \\
\Rightarrow & p_n = d_1 (1)^n + d_2 \big(\frac{\lambda}{\mu}\big)^n = d_1 + d_2 \rho \\
\end{align}
$$
* 解到這裡還有 $d_1$ 和$d_2$ 這兩個未知數，可以使用 $p_1 = \rho p_0$ 這個 boundary condition 推得
$$
p_1 = d_1 + d_2 \rho = \rho p_0, \; \therefore d_1 = 0, d_2 = p_0 \\
\Rightarrow p_n = p_0 \rho^n
$$
* 上式中 $p_0$ 依然是未知數，可以再透過 $\sum p_n = 1$ 的 boundary condition 求得
$$
\sum_{n=0}^{\infty}p_n = p_0 + \rho p_0 + \rho^2 p_0 + ... = \Big(\frac{1}{1 - \rho}\Big) p_0 = 1 \\
\Rightarrow p_0 = 1- \rho \quad (\rho < 1)
$$
* 綜合上述，終於可以求得 $p_n$ 的一般式為
$$
p_n = (1 - \rho)\rho^n
$$

## 3.2.4 Measures of Effectiveness

* 在解出 M/M/1 系統的 steady-state distribution 之後，我們就可以利用他來計算系統的各種效能指標。一般來說我們會在意的效能指標有兩種：和數人頭有關的，和時間有關的。數人頭有關的比較簡單，通常只是在算級數和；和時間有關的就比較難，通常會牽扯到積分式
* Mean system size: 平均來說系統中有多少人（包含隊伍和 server 中的人）。簡單來說就是系統人數的期望值，也就是系統處於 $n$ 個人數的狀態的機率乘上 $n$ 個人數
$$
L = \sum_{n=0}^{\infty} n p_n = \sum_{n=0}^{\infty} n \rho^n (1 - \rho) = (1 - \rho) \sum_{n=0}^{\infty} n \rho^n = (1 - \rho) \cdot \frac{\rho}{(1 - \rho)^2} = \frac{\rho}{1 - \rho}
$$
* Mean queue size: 平均來說隊伍中有多少人。因為 M/M/1 系統只有一個 server，所以隊伍人數永遠都會是系統人數減 1，因此就是要計算系統處於 $n$ 個人數的狀態的機率乘上 $n-1$ 個人數
$$
L_q = \sum_{n=1}^{\infty} (n-1) p_n = \sum_{n=0}^{\infty} np_n - \sum_{n=1}^{\infty} p_n = L - (1 - p_0) = L - \rho = \frac{\rho^2}{1 - \rho}
$$
* Mean queue size for nonempty queues: 在隊伍非空的前提下（也就是系統中至少有兩個人的前提下），平均隊伍中有多少人
* 必須先定義隊伍非空前提下的系統人數條件機率 $p_n^{'}$
$$
\begin{align}
p_n^{'} &= Pr\{n \text{ events in system} | n \geq 2\} \\
&= \frac{Pr\{n \text{ events in system and } n \geq 2\}}{Pr\{n \geq 2\}} \\
&= \frac{p_n}{1 - p_0 - p_1} = \frac{p_n}{\rho^2}
\end{align}
$$
* 隊伍非空的前提下，平均隊伍長度就是
$$
L_q^{'} = E[N_q | N_q \neq 0] = \sum_{n=2}^{\infty} (n-1)p_n^{'} = \frac{1}{1 - \rho}
$$

## 3.2.5 Waiting time distribution

* ~~要計算跟時間有關的效能指標實在太難了，難到值得我們把他獨立出來變另一小節~~
* 在計算跟時間有關的效能指標的時候，往往需要先求出該時間效能指標的機率密度函數，才可以藉由機率密度函數去求取平均值（期望值）。而由於時間是一個連續的變數，所以要求取機率密度函數時，往往會牽扯到積分式
* Waiting time: 顧客等在隊伍當中的平均時間長度，需要先求出在隊伍中等待的時間 $W_q$ 的機率密度函數才有辦法求取
* 令 $T_q$ 為一個隨機變數，代表一個顧客待在隊伍裡面的等待時間
* 當某個顧客到達隊伍時，發現系統裡已經有 $n$ 個人的機率為 $q_n \equiv Pr\{n \text{ in system at an arrival}\}$
* 相較之下， $p_n$ 是不限定任何的時間點，觀察到現在系統當中有 $n$ 個人的機率（$p_n$ 和 $q_n$ 的關係有點類似 CTMC 和他的 embedding DTMC 之間的關係）
* 在一個 M/M/1 的排隊系統中，$q_n = p_n$。因為 M/M/1 的抵達過程是 poisson process，具有 memoryless property，所以不管什麼時候去觀測系統人數，都跟當下是什麼時間點沒有關係
* 在滿足以下條件的時候，$q_n$ 都會等於 $p_n$
1. 隊伍的空間無限大，不會有任何顧客想進去進不去
2. 系統的抵達率必須是常數，不是隨時間變化的函數
* 因為 M/M/1 $q_n = p_n$，所以等等在計算 $W_q(t)$ 的時候，有看到 $q_n$ 都可以改寫成 $p_n$
* 系統的排隊時間分佈 $W_q(t)$ 就會是 $T_q$ 的 CDF。為了求出 $W_q(t)$，我們需要分成兩個時間區間討論
* 在 $t = 0$ 的時候，一定是不用等的
$$
\begin{align}
W_q(0) &= Pr\{T_q \leq 0\} = Pr\{T_q = 0\} \\
&= Pr\{\text{System is empty at an arrival}\} \\
&= q_0 = p_0 = 1 - \rho
\end{align}
$$
* 在 $t > 0$ 的時候，可以拆成一定要等和一定不用等的兩種狀況相加
$$
\begin{align}
W_q(t) &= Pr\{T_q \leq t\} \\
&= \sum_{n=1}^{\infty}Pr\{T_q \leq t | \text{arrival found n in system}\} \times Pr\{\text{arrival found n in system}\} + Pr\{T_q = 0\}\\
&= \sum_{n=1}^{\infty}Pr\{T_q \leq t | \text{arrival found n in system}\} \times q_n + W_q(0)\\
&= \sum_{n=1}^{\infty}Pr\{\text{time for n completions } \leq t | \text{arrival found n in system}\} \times p_n + W_q(0) \\
&= \sum_{n = 1}^{\infty} (1-\rho) \rho^n \int_0^t \frac{\mu(\mu x)^{n-1}}{(n-1)!} e^{-\mu x} dx + (1-\rho) \\
&= (1-\rho)\rho \int_0^t \mu e^{-\mu x} \sum_{n = 1}^{\infty} \frac{(\rho \mu x)^{n-1}}{(n-1)!} dx + (1-\rho) \\
&= (1-\rho) \rho \int_0^t \mu e^{-\mu x(1-\rho)} dx + (1 - \rho) \\
&= \rho[1 - e^{-\mu (1-\rho)t}] + (1-\rho) \\
&= 1 - \rho e^{-\mu (1-\rho)t}
\end{align}
$$
* 一定要等的狀況下所需的等待時間不超過 $t$ 的機率，就是「顧客在抵達系統時發現已經有 $n$ 個人在裏面，而且等待時間還是不超過 $t$」的機率，也就是「當某個顧客抵達系統時發現已經有 $n$ 個人在系統中了的前提下，等待時間依然沒有超過 $t$ 的條件機率」，乘上「顧客抵達系統時發現已經有 $n$ 個人在系統中」的機率，而後者就是 $q_n$
* 一定不用等的情況就跟 $W_q(0)$ 是一樣的
* 「有 $n$ 個人排在我前面的所需等待時間」，可以等價想成「前面 $n$ 個人的服務時間加總」，因為在一個先進先出的系統當中，我一定得等到我前面的人都被服務完了，才會停止等待
* 在 M/M/1 系統當中，因為每個顧客的服務時間都呈現指數分佈，所以「前 $n$ 個人的服務時間加總」，就是在加總 $n$ 個指數分佈的隨機變數，也就是 Erlang-type n 的分佈
* $\sum_{n=1}^{\infty} \frac{(\rho \mu x)^{n-1}}{(n-1)!}$ 是 $e^{\rho\mu x}$ 的泰勒展開式
* $\mu (1-\rho) e^{-\mu (1-\rho) x}$ 就是參數為 $\mu (1-\rho)$ 的指數分佈的 PDF，他的積分就是該指數分佈的 CDF，也就是 $1 - e^{-\mu (1-\rho)t}$
* 現在有 $W_q(t)$ 了，他的平均值 $W_q$ 就是平均等待時間，同時也是在算 $T_q$ 的期望值
$$
\begin{align}
W_q &= E[T_q] = \int_0^{\infty} t dW_q(t) \\
&= 0 \cdot (1 - \rho) + \int_0^{\infty} t dW_q(t) \\
&= \int_0^{\infty} t [\rho \mu (1 - \rho) e^{-\mu (1 - \rho)t}] dt \\
&= \rho \int_0^{\infty} t [\mu (1 - \rho) e^{-\mu (1 - \rho)t}] dt \\
&= \rho \cdot \frac{1}{\mu (1 - \rho)} = \frac{\rho}{\mu (1 - \rho)}
\end{align}
$$
* System time: 顧客等在系統當中待的平均時間長度，需要先算出 $W(t)$ 才能求得。方法跟計算 $W_q$ 差不多，唯一的差別在於 $W_q$ 有可能會是 0，但 $W$ 不可能是 0，所以不需要分成兩種請況計算
* 定義 $T$ 為一個隨機變數，代表一個顧客待在系統裡面的時間
* $W(t)$ 就是 $T$ 的 CDF，計算方式同 $W_q$，但不需要加上不用等的那一項
$$
\begin{align}
W(t) &= Pr\{T \leq t\}
= \sum_{n=0}^{\infty} Pr\{T \leq t | \text{arrival found n in system} \} \times p_n \\
&= \sum_{n=0}^{\infty} (1 - \rho) \rho^n \int_0^t \frac{\mu(\mu x)^n}{n!} e^{-\mu x} dx \\
&= (1-\rho) \int_0^t \mu e^{-\mu x} \sum_{n=0}^{\infty} \frac{(\rho \mu x)^n}{n!}dx \\
&= \int_0^t \mu (1 - \rho) e^{-\mu (1-\rho)x} dx = 1 - e^{-\mu (1-\rho)t}
\end{align}
$$
* 觀察 $W(t)$ 的長相，會發現他就直接是一個指數分佈了，所以計算他的平均值不需要用到積分，可以直接透過指數分佈的性質知道他的平均值 $W$ 是
$$
W = \frac{1}{\mu - \lambda}
$$
* 有時候我們不想要知道這些時間相關的效能指標的完整分佈函數，我們只是想要知道均值，此時可以利用 [Little's Law](https://hackmd.io/kfLYKMjuTY2VaXDuTbsw8Q?view#14-Little%E2%80%99s-Law)，透過系統人數 $L$ ($L_q$) 和抵達率 $\lambda$ 來快速求取 $W$ ($W_q$) 的平均值