筆記 - 網路效能模擬與分析

發表於 更新於 分類於 Disqus：

  此為「網路效能分析與模擬」學習筆記。主要為discrete-event simulation的筆記，參考書目為A. M. Law的Simulation Modeling & Analysis, 4ed。

Next-event time advance

  Next-event time advance的模擬方法，當有事件發生時simulation clock才會前進到事件發生的時間點，因此實際模擬時間可以大幅縮短。


  這裡以single-server queueing system為例，提出模擬需要的變數：

  • $t_i$ = time of arrival of ith customer
  • $A_i$ = $t_i - t_{i -1}$ = inter-arrival time
  • $S_i$ = service time of ith customer
  • $D_i$ = delay in queue of ith customer
  • $c_i$ = time of departure of ith customer
  • $e_k$ = time of occurance of kth event

  其中event可以是顧客的抵達、顧客的離開(結束服務後)，因此 $e_i$ 就是simulation clock的時間點。這些變數大多是隨機變數，因此模擬時需要使用到random number generator。


  下圖說明next-event time-advance模擬方法的流程圖。注意模擬流程中使用event-scheduling。模擬會在滿足某個條件的時候停止，比如說當queue裡的人數超過某個數字，模擬就會停止，或是某個顧客的delay time超過一定值，模擬就停止。


圖片來源：自己，改自參考書目。

  系統中的entities：servers and customers，server的attributes：server status(idle or busy)、queue length；customers的attributes：抵達時間、離開時間、排隊時間等。

Single-server queue model

  假設inter-arrival time $A_i$ 為i.i.d.的隨機變數，service time $S_i$ 也是i.i.d.隨機變數，且 $A_i$ 和 $S_i$ 互相獨立。如果server正在服務前一位顧客，則抵達的顧客必須排隊，遵循FIFO的原則，先來到先服務。


  這個系統有以下效能指標：

  • Average delay in queue $d(n)$

      基本上是一個隨機變數，因為每次模擬得到的average(over customers)都不一樣。其估計量(estimate)為

    $$
    \hat{d}(n) = \frac{\sum_{i = 1}^{n}{D_i}}{n}
    $$

  • Average number of customers in queue $q(n)$

      排隊人數(或者說queue length)的期望值，也是一個隨機變數。其估計量為

    $$
    \begin{align}
    \hat{q}(n) &= \frac{\sum_{i = 0}^{\infty}{i T_i}}{T(n)}\\[10pt]
    &= \frac{\int_{0}^{T(n)}{Q(t) \rm{d} t}}{T(n)}
    \end{align}
    $$

      其中 $T(n) = T_0 + T_1 + \cdots$。$Q(t)$ 為時間t的queue length。

  • Utilization $u(n)$

      使用率，指server忙碌的時間佔所有時間($T(n)$)的比例。定義busy function為

    $$
    B(t) =
    \begin{cases}
    1,\quad\text{if server busy at time t}\\
    0,\quad\text{if server not busy at time t}
    \end{cases}
    $$

      則使用率的估計量為

    $$
    \hat{u}(n) = \frac{\int_{0}^{T(n)}{B(t) \rm{d} t}}{T(n)}
    $$

  下圖顯示single-server queue model的模擬模型，為系統初始時的snapshot。其中event list的 A 表示arrival，D 表示departure；# delayed 表示離開queue的人數，A under Q(t) 表示Q(t)函數下的面積。


圖片來源：自己，改自參考書目。

  要注意 times of arrival 記錄的不是所有顧客的arrival time，而是目前在queue的顧客的arrival time。模擬的停止條件可以是，當N個人離開queue──也就是 # delayed = N 時，模擬就停止。


  Tie-breaking：同時發生departure和arrival，哪一個要先處理？這個現象叫做time ties，需要採取tie-breaking rule來決定哪一個先處理。


使用event graph的例子。

Inventory system

  要解決的問題為：決定接下來 $n$ 個月每個月要訂多少貨。Time between demand (TBD) 為i.i.d.指數隨機變數，每次的需求量也是一個離散隨機變數。用 $\text{K} + i\text{Z}$ 表示訂貨成本，$\text{K}$ 為overhead，$Z$ 為訂貨量，$i$ 即為單位訂貨量成本。Lead time也model為一個隨機變數。而最關鍵的訂貨策略為一stationary policy $(s, S)$，表示如下：

$$
Z =
\begin{cases}
\begin{align}
\text{S} - I,\quad &\text{if}\ I < s\\[10pt]
0,\quad &\text{if}\ I \geq s
\end{align}
\end{cases}
$$

  如果inventory level小於到來的demand，就會有backlog(積欠)，inventory level必須減去積欠的量，所以會是負值。之後訂貨到貨時，先扣掉backlog的量，剩餘的量才會加入inventory level。


  此外，用 $I^{+}(t) = \max{\{I(t), 0\}}$ 表示inventory裡實際存有的量，$I^{-}(t) = \max{\{-I(t), 0\}}$ 表示backlog的量。假設每個月的holding cost為每單位貨量 $h$ 元，則每個月的平均holding cost為 $h\overline{I^{+}}$，其中

$$
\overline{I^{+}} = \frac{\int_{0}^{n}{I^{+}(t) dt}}{n}
$$

  為inventory的每月平均實際存量。同樣的，也可以計算每月平均積欠量：

$$
\overline{I^{-}} = \frac{\int_{0}^{n}{I^{-}(t) dt}}{n}
$$

  然後得到每月平均backlog cost (shortage cost)為 $\pi\overline{I^{-}}$。Event graph、狀態變數的部分，請從C程式來推敲、對照。

Modeling (more) complex systems

  基本的做法是：把系統分割成幾個子系統，再來model子系統間的互動。


  以tandem queue為例：當後一級的buffer size滿了，前一級的departure就無法排程，此為blocking。除了原本的 $M/M/1$ model外，還要新增event、狀態變數來處理blocking。


  再一例，network example，如下圖。圖中的路徑(route) $L_i$ 可以視為 $M/M/1$ 裡的server，N2這個節點可以視為end-service的event。


Open network model。圖片取自上課講義。

  本學期project主題，簡化的video capture server模型就是使用到tandem queue的例子，只不過因為後一級的server buffer容量為無限大，所以不需要模擬blocking的情形。

補充：模擬常用到的隨機分布

  首先是連續隨機變數。

  • Exponential：應該很熟悉了，可以用來描述interarrival time和time to failure。
  • Gamma：如果有 $m$ 個i.i.d.的exponential隨機變數 $X_i$，則 $X_1 + X_2 + \cdots + X_m$ 是一個gamma($m$, $\beta$)隨機變數。可以用來描述machine repair time或task completion time。
  • Weibull：Rayleigh分布為其特例，可以用來描述task completion time、time to failure。
  • Normal：用來描述許多隨機事件的總和(CLT)。
  • Log-normal：用來描述許多隨機變數的乘機(CLT)。
  • Beta；常用來描述隨機的比例，比如說半成品中瑕疵品的比例，或是用在計畫完成時間的modeling，如PERT network的應用。

  接著是離散隨機變數。

  • Bernoulli：用來描述「成功或失敗」的隨機實驗結果(Bernoulli trial)。
  • Binomial：用來描述 $N$ 個獨立的Bernoulli trial中成功的次數。
  • Geometric：用來描述連續 $N$ 個獨立Bernoulli trial結果為失敗。
  • Poisson：用來描述固定一段時間內某事件的發生次數，如inventory的需求物件數。

  最後談談empirical distribution。如果data無法擬合到任何理論上的機率分布(e.g. 資料量太少)，可能就要用empirical distribution來描述。假設我有 $N$ 筆資料，且為遞增排序，$X_1 \leq X_2 \leq \cdots \leq X_N$，則累積分布函數(cdf)可以近似為

$$
F(x) =
\begin{cases}
\begin{align}
&0 &{\rm if}\ x < X_1\\[10pt]
\frac{i - 1}{N - 1} + &\frac{x - X_i}{(N - 1)(X_{i + 1} - X_i)} &{\rm if}\ X_i \leq x < X_{i + 1}\\[10pt]
&1 &{\rm if}\ X_N \leq x
\end{align}
\end{cases}
$$

  如果各個資料點的數值不得而知、只有整理好的histogram，也可以用類似的方法(interpolating)來建立empirical distribution。假設我有 $N$ 筆資料，分類到 $[a_0, a_1),\ [a_1, a_2),\ \cdots,\ [a_{k - 1}, a_k)$ 共 $k$ 個bins (intervals)，每一個bin包含 $N_i$ 個資料點。定義 $G(a_j) = (N_1 + N_2 + \cdots + N_j)/N$，則cdf可以近似為

$$
G(x) =
\begin{cases}
\begin{align}
&0 &{\rm if}\ x < a_0\\[10pt]
G(a_{j - 1}) + &\frac{x - a_{j - 1}}{a_j - a_{j - 1}}[G(a_j) - G(a_{j - 1})] &{\rm if}\ a_{j - 1} \leq x < a_j\\[10pt]
&1 &{\rm if}\ a_k \leq x
\end{align}
\end{cases}
$$

  有了cdf，就可以透過uniform distribution來生成empirical distribution的隨機變數值(inverse transform)。

補充：如何生成隨機變數值

  通常我們可以生成uniform分布的隨機變數值，但是要如何生成任一機率分布的隨機變數值？這節就要來回答這個問題。


  舉例來說，標準常態分布 $N(\mu, \sigma^2)$ 的隨機變數值怎麼產生？最早的Box-Muller方法，假設 $U_1$ 和 $U_2$ 為 $U(0, 1)$ 分布，則

$$
\begin{align}
&X_1 = \sqrt{-2 \ln{U_1}}\cos{2\pi U_2}\\[10pt]
&X_2 = \sqrt{-2 \ln{U_1}}\sin{2\pi U_2}
\end{align}
$$

  不過，這種方法要求 $U_1$ 和 $U_2$ 必須為 i.i.d.，而且三角函數的計算也使得這種方法運算速度較慢。下面介紹兩種較常使用的方法。

Inverse Transform

  假設一個機率分布其cdf $F(x)$ 為strictly increasing (因此有反函數)，則可以透過uniform distribution來生成 $X$ 的隨機變數值：

  1. 生成 $U \sim U(0, 1)$ 的隨機變數值 $u$
  2. $x = F^{-1}(u)$ 即為所求

  我們可以驗證得到的 $X$ 確實為目標機率分布：

$$
P[X \leq x] = P[F^{-1}(U) \leq x] = P[U \leq F(x)] = F(x)
$$

  離散隨機變數值也可以用這個方法來生成，只不過inverse mapping的部分就比較麻煩。更一般性的作法為：

$$
X = \min{\{x:\ F(x) \geq U \}}
$$

  Inverse transform的缺點是必須知道目標機率分布的cdf。

Acceptance-Rejection Method

  假設所求的機率分布pdf為 $f(x)$。令 $t(x)$ 為majorizing function，$t(x) \geq f(x)\ \forall x$。定義 $r(x) = t(x)/c$，其中 $c = \int_{-\infty}^{\infty}{t(x) dx}$，因此 $r(x)$ 是一個合法的pdf。則生成 $X$ 隨機變數值的方法為：

  1. 生成 $Y \sim r$ 的隨機變數值 $y$
  2. 生成 $U \sim U(0, 1)$ 的隨機變數值 $u$，$U$ 獨立於 $Y$
  3. 若 $u \leq f(y)/t(y)$，則 $x = y$；否則，回到第一步再試一次。

  這應該是一種蒙地卡羅的方法。接著我們來驗證得到的 $X$ 確實為目標機率分布。定義 $A$ 為「接受生成的 $Y$ 值」的事件，$A$ 發生時 $X = Y$，因此

$$
P[X \leq x] = P[Y \leq x\ |\ A] = \frac{P[Y \leq x, A]}{P[A]}
$$

  不過

$$
P[A\ |\ Y = y] = P[U \leq \frac{f(y)}{t(y)}] = \frac{f(y)}{t(y)}
$$

  因此 $P[Y \leq x\ |\ A]$ 的分子部分為

$$
\begin{align}
P[Y \leq x, A] &= \int_{-\infty}^{x}{P[Y \leq x, A\ |\ Y = y]r(y) dy}\\[10pt]
&= \int_{-\infty}^{x}{P[A\ |\ Y = y]\frac{t(y)}{c} dy}\\[10pt]
&= \frac{1}{c}\int_{-\infty}^{x}{f(y) dy}
\end{align}
$$

  分母的部分則是 $P[A] = 1/c$。因此

$$
P[X \leq x] = \frac{P[Y \leq x, A]}{P[A]} = \frac{\frac{1}{c}\int_{-\infty}^{x}{f(y) dy}}{1/c} = \int_{-\infty}^{x}{f(y) dy}
$$

  也就是 $X$ 的cdf，故得證。

Internet of Things (IoT)

  老師說他想講IoT，雖然和網路效能分析好像沒什麼關係。這部分請參考review paper [1]。

Summary

  修課心得：老師人好，也學到一些離散事件模擬的概念。至於隨機亂數產生器和IoT的部分，我沒有太認真聽，就不寫這部分的筆記。總之還是有學到東西的。

Reference

[1] C. Milarokostas et. al, “A Comprehensive Study on LPWANs With a Focus on the Potential of LoRa/LoRaWAN Systems,” in IEEE Communications Surveys & Tutorials, vol. 25, no. 1, pp. 825-867, Firstquarter 2023.