具有工作故障的流模型近似解分析

Asymptotic Solution of the Fluid Model with Working Breakdown Policy

摘要: 本文分析了具有工作故障策略的流模型。当流模型处于正常工作状态时，服务台以较高速率为顾客提供服务，若发生故障，服务台将进入工作故障状态，同时维修立刻开始，期间服务台服务速率降低。当维修完成时，流模型将进入正常工作状态。通过归一化技术(Uniformization Technique)与递推公式，得到了流模型库存量尾分布函数的近似表达式，并得到了稳态下库存量各阶矩的表达式。最后通过数值例子分析了系统性能指标随参数的变化趋势。

Abstract: In this paper, we analyze the fluid model with working breakdown policy. During the normal work-ing period, the service serves customers at a high rate. When breakdown happens, server will turn into the partial working period and continue to provide service to arriving customers, but with a lower rate. When the repair is completed, the fluid model will enter the normal working period. Through the recursive formula and Uniformization Technique, we obtain the recursive expression of the tail joint distribution function and the expressions of each moment of the buffer content in stability condition are obtained. Finally, the variation trend of system performance index with pa-rameters is analyzed by numerical examples.

1. 引言

流模型作为传统离散时间排队模型的推广，已经广泛应用于现实生产生活中，国内外众多学者已经对其进行了广泛而深入的研究。Mitra [1] 利用流模型对制造系统中出现的排队现象进行建模，在提出了一种新的优化方案的同时利用可逆马尔可夫过程的特征值提出了适用范围更广的流模型。Bekker与Mandjes [2] 利用流模型对计算机通信系统进行建模，得出了系统的尾部渐近性能指。Gautam [3] 系统地总结了流模型的理论与其应用。在研究由具有有限状态不可约马尔可夫链驱动的流模型时Nabli [4] 提出了一种用于逼近流模型库存量尾分布函数的地推算法。利用该算法，Nabli [5] [6] [7] 研究了许多由马尔可夫链驱动的流模型，并对流模型的各阶矩进行了分析。

完全可靠的服务台在现实生产生活中几乎不存在，因此，在过去几十年中已有众多学者对这类具有不可靠服务台的排队系统进行了研究。早期对具有工作故障策略的排队系统的研究可查阅Kalidass与Kastrui [8] 及其参考文献。随后众多学者研究了具有各种特征的不可靠排队系统 [9] [10] [11] 。但是注意到很少有学者研究由具有工作故障策略的排队系统驱动的流模型。孙红霜与叶晴晴 [12] 由具有工作故障策略的M/M/1重试排队系统驱动的流模型，通过矩阵几何解法与Laplace-Stieltjes变换得到了该流模型稳态下的平均库存量。

本文受工作故障策略与流模型的启发，研究了一个具有工作故障的流模型。在此模型中，流模型的服务台拥有两种状态，分别为工作故障状态与正常工作状态。当流模型处于正常工作状态时，服务台以较高速率为顾客提供服务且有概率发生故障；若发生故障，服务台将进入工作故障状态，同时维修立刻开始，期间服务台服务速率减低。当维修完成时，流模型将进入正常工作状态。

本文剩余部分结构如下：第2节中详细描述了具有工作的流模型，后构建了该模型库存量的尾分布函数；第3节利用归一化技术与递推公式，得到了该流模型库存量尾分布函数的近似表达式，同时给出了在给定误差限下所需截断步长的计算方法；第4节给出了该流模型库存量各阶矩的计算公式；第5节通过数值例子，分析了系统性能指系统参数的变化趋势。

2. 流模型描述

考虑具有工作故障的流模型，流模型可以假设如下：

1) 顾客到达的时间间隔服从参数为 $\lambda$ 的指数分布，即流模型的输入速率为 $\lambda$ 。

2) 服务台在运行过程中会出现故障，且在正常工作状态下服务台寿命服从参数为 $\alpha$ 的指数分布，即平均寿命时间为 $1/\alpha$ ；修理时间服从参数为 $\beta$ 的指数分布，即平均修理时间为 $1/\beta$ 。

3) 在正常工作状态下，服务台服务速率为 $\mu$ ，即正常工作状态下流模型的输出速率为 $\mu$ 。当故障发生时，服务器将进入工作故障状态且继续为到达的客户提供服务，但服务速率降低为 $\eta \left(\eta <\mu \right)$ ，即部分工作期间流模型的输出速率为 $\eta$ 。

4) 假设顾客到达的时间间隔，正常工作期间的服务速率，部分工作期间的服务速率，故障发生时间与维修时间相互独立。此外，假设服务顺序为先进先出(FIFO)。

显然，该流模型的状态可以用马尔可夫链 $X=\left(X\left(t\right)\right)$ 来表示，其中 $X\left(t\right)$ 为时刻t流模型所处的状态，且

$X\left(t\right)=\{\begin{array}{c}0,流模型处于工作故障状态\\ 1,流模型处于正常工作状态\end{array}$

其状态空间为 $S=\left(0,1\right)$ ，对应的无穷小生成元为：

$A=\left(\begin{array}{cc}-\beta & \beta \\ \alpha & -\alpha \end{array}\right)$ (1)

考虑具有两阶段故障的流模型的库存量。假设流模型的输入速率可以用顾客到达速率表示，输出速率可以用服务速率表示。令 $Q\left(t\right)$ 表示时刻t的库存量，且为非负随机变量。设缓冲器的净输入速率(输入速率–输出速率)为二维随机程 $\left\{\left(Q\left(t\right),X\left(t\right)\right),t\ge 0\right\}$ 的函数，且满足

$\frac{\text{d}Q\left(t\right)}{\text{d}t}=\{\begin{array}{l}\lambda -\eta ,X\left(t\right)=0,Q\left(t\right)\ge 0\\ 0,X\left(t\right)=1,Q\left(t\right)=0\\ \lambda -\mu ,X\left(t\right)=1,Q\left(t\right)>0\end{array}$ (2)

记净输入速率为 $d_i$ ， $i=0,1$ 且 $d_0=\lambda -\eta$ ， $d_1=\lambda -\mu$ 。式(2)给出了该流模型处于不同状态下库存量的变化趋势。具体来说，当流模型处于工作故障状态时，库存量随速率 $\left|d_0\right|$ 线性增加；当流模型处于正常工作状态，且库存量为0时，库存量将保持不变；当流模型处于正常工作状态，且库存量大于0时，库存量随速率 $\left|d_1\right|$ 线性增加。此外，定义依赖于净输入速率的状态空间

$S^{+}=\left\{i\in S|d_i>0\right\},S^{-}=\left\{i\in S|d_i\le 0\right\}$ (3)

即： $S^{+}=\left\{0\right\}$ ， $S^{-}=\left\{1\right\}$ 。

定义 $\pi =\left({\pi }_0,{\pi }_1\right)$ 为 $X=\left(X\left(t\right)\right)$ 的稳态概率。行向量 $\pi$ 满足

$\pi A=0,\pi 1^{\text{T}}=1$ (4)

其中 $0$ 与 $1$ 分别是由0和1构成的二维行向量，符号T为转置算子。该流模型的稳态条件为

$\rho =\lambda /\mu <1,d_0{\pi }_0+d_1{\pi }_1<0$ (5)

因此，过程 $\left\{\left(Q\left(t\right),X\left(t\right)\right),t\ge 0\right\}$ 是一个二维马尔可夫过程，定义该马尔可夫过程在时刻t的尾联合分布函为

$F_{i,t}\left(x\right)=P\left\{Q\left(t\right)>x,X\left(t\right)=i\right\}$ (6)

当稳态条件(5)成立时，定义过程 $\left\{\left(Q\left(t\right),X\left(t\right)\right),t\ge 0\right\}$ 的极限分布为 $\left\{Q_{\infty },X_{\infty }\right\}$ 。令

$F_i\left(x\right)=\underset{t\to \infty }{\lim }P\left\{Q\left(t\right)>x,X\left(t\right)=i\right\}=P\left\{Q_{\infty }>x,X_{\infty }=i\right\},i=0,1$ (7)

为了表示方便，引入矩阵

$F\left(x\right)=\left(F_0\left(x\right),F_1\left(x\right)\right)$ (8)

$D=diag\left(d_0,d_1\right)$ (9)

由文献 [13] ，可知，行向量 $F\left(x\right)$ 满足下列微分方程

$F^{\prime }\left(x\right)D=F\left(x\right)A$ (10)

其中 $F^{\prime }\left(x\right)$ 为 $F\left(x\right)$ 在x处的导数。为了求解尾分布函数 $F_i\left(x\right)$ ，需要求解微分方程(10)，同时确定初始条件 $F\left(0\right)$ 。当 $F_i\left(x\right)$ 确定时，流模型库存量在稳态条件下的尾分布函数可由全概率公式计算得到

$P\left\{Q_{\infty }>x\right\}=F_0\left(x\right)+F_1\left(x\right)=F\left(x\right)1^{\text{T}}$ (11)

同时，该流模型的平衡方程(输入速率 = 输出速率)可表示为

$\lambda P\left\{Q_{\infty }=0,X_{\infty }=1\right\}+\mu P\left\{Q_{\infty }>0,X_{\infty }=1\right\}+\eta P\left\{Q_{\infty }\ge 0,X_{\infty }=0\right\}=\lambda$ (12)

3. 模型近似解

在本节中，将使用单值化技术求解该流模型尾分布函数的近似解。令 $a=\max \left(-a_{ii}\right)$ ，其中 $a_{ii}$ 为矩阵 $A$ 的第i个对角线元素。同时令 $d=\min \left\{d_i|d_i>0\right\}$ ，即 $d=d_0$ 。随后定义

$P=\frac{A}{a}+I$ (13)

其中： $I$ 是一个二维单位矩阵。

显然，矩阵 $P$ 是随机的，即矩阵 $P$ 的所有元素都是非负的且每一行和为1。此外，矩阵 $P$ 满足

$\pi P=\pi$ (14)

由式(2)与稳态条件(5)可知 $d_0>0$ ， $d_1<0$ ，故矩阵 $D$ 是可逆的，由式(10)，得

$F\left(x\right)=F\left(0\right)\exp \left(xA{D}^{-1}\right)$ (15)

随后，需要确定初始条件 $F\left(0\right)$ ，由式(3)，可知当 $X\left(t\right)=0$ 时，库存量 $Q\left(t\right)$ 以概率1满足 $Q\left(t\right)>0$ ，故 $F_1\left(0\right)={\pi }_1$ ，此外由平衡方程(12)，得 $F_2\left(0\right)=-d_0{\pi }_0/d_2$ 。在下面的定理中，本文将给出一个算法来计算向量 $F\left(x\right)$ 。

定理1 当稳态条件(5)成立时，对于任意 $x\ge 0$ ，尾联合分布函数 $F_i\left(x\right)$ 满足

$F_i\left(x\right)={\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}{b}_i\left(k\right)}$ (16)

其中系数 $b_i\left(k\right)$ 满足：

当 $k=0$ 时

$b\left(0\right)=F(\; 0\; )$

其中

$b\left(0\right)=(b_0\left(0\right),b_1(\; 0\; ))$

当 $k\ge 1$ 时

$\{\begin{array}{l}{b}_0\left(k\right)=\left(1-\frac{d}{d_0}\frac{\beta }{a}\right)b_0\left(k-1\right)+\frac{d}{d_0}\frac{\alpha }{a}{b}_1\left(k-1\right)\\ b_1\left(k\right)=-\frac{d_0{b}_0\left(k\right)}{d_1}\end{array}$

此外， $b_i\left(k\right)$ 是非负的。

证明矩阵 $A{D}^{-1}$ 可以写成分块矩阵的形式

$A{D}^{-1}=-\frac{a}{d}I+\frac{a}{d}G$ (17)

其中

$G=\left(P-I\right)d{D}^{-1}+I$

故

$F\left(x\right)=F\left(0\right)\exp \left(xA{D}^{-1}\right)={\text{e}}^{-\left(ax/d\right)}F\left(0\right)\exp \left(\frac{ax}{d}G\right)={\text{e}}^{-\left(ax/d\right)}{\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}\frac{{\left(ax/d\right)}^k}{k!}b\left(k\right)}$ (18)

其中 $b\left(k\right)=F\left(0\right)G^k$ ，同时得

$\{\begin{array}{l}b\left(0\right)=F\left(0\right)\\ b\left(k\right)=b\left(k-1\right)G\end{array}$

即

$\{\begin{array}{l}{b}_i\left(0\right)=F_i\left(0\right),i\ge 0\\ b_0\left(k\right)=\left(1-\frac{d}{d_0}\frac{\beta }{a}\right)b_0\left(k-1\right)+\frac{d}{d_0}\frac{\alpha }{a}{b}_1\left(k-1\right)\\ b_1\left(k\right)=\frac{d}{d_1}\frac{\beta }{a}{b}_0\left(k-1\right)+\left(1-\frac{d}{d_1}\frac{\alpha }{a}\right)b_1\left(k-1\right)\end{array}$

注意到， $d/d_1<0$ 而 $1-\alpha d/a{d}_1>1$ ，在进行计算时可能会发生正负抵消而使得算法不稳定。为了避免这一问题，给出另一个与 $i=1$ 有关的关系式，在这个式子中将只出现正数。

定义行向量 $d=\left(d_0,d_1\right)$ ，得

$\left[\left(P-I\right)d{D}^{-1}+I\right]d^{\text{T}}=d\left(P-I\right)1^{\text{T}}+d^{\text{T}}=d^{\text{T}}$ (19)

当 $k\in \mathbb{N}$ 时， $b\left(k\right)d^{\text{T}}=b\left(k-1\right)G{d}^{\text{T}}=b\left(k-1\right)d^{\text{T}}$ 。考虑到 $b\left(0\right)=F\left(0\right)$ ，得

$b\left(k\right)d^{\text{T}}=b\left(0\right)d^{\text{T}}=F\left(0\right)d^{\text{T}}=d_0{F}_0\left(0\right)+d_1{F}_1\left(0\right)=d_0{\pi }_0+d_1\left(-\frac{d_0{\pi }_0}{d_1}\right)=0$ (20)

由此，得到 $d_0{b}_0\left(k\right)+d_1{b}_1\left(k\right)=0$ ，即 $b_1\left(k\right)=-d_0{b}_0\left(k\right)/d_1$ 。

由参数d与 $\alpha$ 的定义，可知 $\alpha /a$ 与 $\beta /a$ 为以1为界的正数且 $d/d_0=1$ ，此外 $1-d\beta /d_0\alpha$ 与 $1-d/d_1$ 是非负的。由于 $b_i\left(0\right)=F_i\left(0\right)$ 且 $F_i\left(0\right)$ ， $i=0,1$ 是非负的，使用归纳法即可证明对于所有 $k\ge 0$ ，有 $b_i\left(k\right)\ge 0$ ， $i=0,1$ 。证毕。

由式(11)，得

$P\left\{Q_{\infty }>0\right\}=b\left(0\right)1^{\text{T}}=\left(d_0-d_1\right){\pi }_0\left(-d_1\right)$ (21)

在下一个定理中，本文将给出数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 的单调性和收敛性。

定理2当 $i\ge 1$ 时， ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 是一个以 ${\pi }_i$ 为界的递减数列，且满足

$\underset{k\to +\infty }{\lim }{b}_i\left(k\right)=0$ (22)

证明首先，通过归纳法证明数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 的单调性。

当 $i=0$ 时

$\begin{array}{c}{b}_0\left(1\right)=\left(1-\frac{d}{d_0}\frac{\beta }{a}\right)b_0\left(0\right)+\frac{d}{d_0}\frac{\alpha }{a}{b}_1\left(0\right)\\ =\left(1-\frac{d}{d_0}\right)b_0\left(0\right)+\frac{d}{d_0}\left(1-\frac{\beta }{a}\right)b_0\left(0\right)+\frac{d}{d_0}\frac{\alpha }{a}{b}_1\left(0\right)\end{array}$ (23)

由于 $b\left(0\right)=F\left(0\right)$ ，得

$b_0\left(1\right)=\left(1-\frac{d}{d_0}\right){\pi }_0+\frac{d}{d_0}\left[\left(1-\frac{\beta }{a}\right){\pi }_0+\frac{\alpha }{a}\left({\pi }_1-P\left\{Q_{\infty }=0\right\}\right)\right]$ (24)

考虑到 $\pi P=\pi$ ，即 ${\pi }_0\left[1-\beta /a\right]+{\pi }_1\alpha /a={\pi }_0$ 得

$\begin{array}{c}{b}_0\left(1\right)=\left(1-\frac{d}{d_0}\right){\pi }_0+\frac{d}{d_0}\left[{\pi }_0-\frac{\alpha }{a}P\left\{Q_{\infty }=0\right\}\right]\\ ={\pi }_0-\frac{d}{d_0}\frac{\alpha }{a}P\left\{Q_{\infty }=0\right\}\\ \le {\pi }_0=b_0\left(0\right)\end{array}$ (25)

当 $i=1$ 时

$\begin{array}{c}{b}_1\left(1\right)=\frac{d}{d_1}\frac{\beta }{a}{b}_0\left(0\right)+\left(1-\frac{d}{d_1}\frac{\alpha }{a}\right)b_1\left(0\right)\\ =\left(1-\frac{d}{d_1}\right)b_1\left(0\right)+\frac{d}{d_1}\frac{\beta }{a}{b}_0\left(0\right)+\frac{d}{d_1}\left(1-\frac{\alpha }{a}\right)b_1\left(0\right)\\ =\left(1-\frac{d}{d_1}\right)b_1\left(0\right)+\frac{d}{d_1}\left[\frac{\beta }{a}{b}_0\left(0\right)+\left(1-\frac{\alpha }{a}\right)b_1\left(0\right)\right]\\ =\left(1-\frac{d}{d_1}\right)\left[{\pi }_1-P\left\{Q_{\infty }=0\right\}\right]b_1\left(0\right)+\frac{d}{d_1}\left[{\pi }_1-P\left\{Q_{\infty }=0\right\}\left(1-\frac{\alpha }{a}\right)\right]\\ ={\pi }_1-P\left\{Q_{\infty }=0\right\}+\frac{d}{d_1}\frac{\alpha }{a}P\left\{Q_{\infty }=0\right\}\end{array}$ (26)

由 $b_1\left(0\right)=F_1\left(0\right)={\pi }_1-P\left\{Q_{\infty }=0\right\}$ 且 $d/d_2<0$ ，得

$b_1\left(1\right)=b_1\left(0\right)+\frac{d}{d_1}\frac{\alpha }{a}P\left\{Q=0\right\}\frac{\alpha }{a}\le b_1\left(0\right)$ (27)

此外，由定理1，得

$\begin{array}{l}{b}_0\left(k+1\right)-b_0\left(k\right)=\left(1-\frac{d}{d_0}\frac{\beta }{a}\right)\left[b_0\left(k\right)-b_0\left(k-1\right)\right]+\frac{d}{d_0}\frac{\alpha }{a}\left[b_1\left(k\right)-b_1\left(k-1\right)\right]\\ b_1\left(k+1\right)-b_1\left(k\right)=-\frac{d_0}{d_1}\left[b_0\left(k+1\right)-b_0\left(k\right)\right]\end{array}$

由参数的定义，可知 $d/d_0=1$ ，且 $\alpha /a\in \left(0,1\right]$ ， $\beta /a\in \left(0,1\right]$ 。由此，可以通过归纳法证明当 $i=0,1$ 时，数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 是单调递减的。接下来，证明数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 的收敛性。由于数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 是非负的单调递减序列，其极限一定存在。令 $l_i$ 为数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 的极限。由式 ，当 $x\to +\infty$ 时，得

$F_i\left(x\right)={\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}{b}_i\left(k\right)}={\displaystyle \underset{k=0}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}{l}_i}$ (28)

由于 $x\to +\infty$ 时，尾分布函数 $F_i\left(x\right)=\underset{t\to \infty }{\lim }P\left\{Q\left(t\right)>x,X\left(t\right)=i\right\}=0$ ，故极限 $l_i$ 必须为0。证毕。

为了计算定理1中给出的无穷级数，本文给出了一个算法，该算法可以保证结果误差不超过给定的误差限 $\epsilon$ 。对于一个给定的正实数x，首先定义泊松级数的截断步骤为 $R\left(x\right)$ ，满足

$R\left(x\right)=\min \left\{n\in N|{\displaystyle \underset{i=0}{\overset{n}{\sum }}{\text{e}}^{-x}\left(\frac{x^i}{i!}\right)}\ge 1-\epsilon \right\}$ (29)

具体的数值计算可参阅文献 [14] ，由于数列 ${\left(b_i\left(k\right)\right)}_{k\in \mathbb{N}}$ 的单调性与收敛性，定义另一个截断步骤

$n_0=\min \left\{R\left(ax/d\right),n_{\epsilon }\right\}$ (30)

其中

$n_{\epsilon }=\min \left\{k\in N|b_0\left(k+1\right)+b_1\left(k+1\right)\le \epsilon \right\}$

由定理1和式 ，得

$P\left\{Q_{\infty }>x\right\}={\displaystyle \underset{k=0}{\overset{n_0}{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}b\left(k\right)1^{\text{T}}}+e\left(n_0\right)$ (31)

同时，误差项 $e\left(n_0\right)$ 满足

$e\left(n_0\right)={\displaystyle \underset{k=n_0+1}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}b\left(k\right)1^{\text{T}}}$ (32)

且，当 $R\left(ax/d\right)<n_{\epsilon }$ 时，由 $b\left(k\right)1^{\text{T}}\le \pi 1^{\text{T}}=1$ ，

$e\left(n_0\right)\le {\displaystyle \underset{k=R\left(ax/d\right)+1}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}}\le 1-{\displaystyle \underset{k=0}{\overset{R\left(ax/d\right)}{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}}\le \epsilon$

当 $R\left(ax/d\right)\ge n_{\epsilon }$ 时，由 $k>n_{\epsilon }$ 时 $b\left(k\right)1^{\text{T}}\le \epsilon$ ，

$e\left(n_0\right)\le {\displaystyle \underset{k=n_{\epsilon }+1}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}b\left(k\right)1^{\text{T}}}\le \epsilon {\displaystyle \underset{k=n_{\epsilon }+1}{\overset{+\infty }{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}}\le \epsilon$

从数值计算的角度来看，若想要计算库存量的尾分布，并使得计算误差在给定的误差限 $\epsilon$ 内，只需要对式 从0到 $n_0$ 求和即可。

4. 库存量的矩

记 $m_i^{\left(n\right)}$ 为 $Q\left(t\right)$ 在状态 $X_{\infty }=i$ 下的n阶矩。由于 $Q\left(t\right)$ 是非负的， $m_i^{\left(n\right)}$ 的表达式为

$m_i^{\left(n\right)}={\displaystyle {\int }_0^{\infty }n{x}^{n-1}P\left(Q_{\infty }>x,X_{\infty }=i\right)\text{d}x}={\displaystyle {\int }_0^{\infty }n{x}^{n-1}{F}_i\left(x\right)\text{d}x}$ (33)

故，流模型稳态条件下库存量 $Q\left(t\right)$ 的n阶矩为

$\begin{array}{c}E\left(Q^{\left(n\right)}\right)={\displaystyle {\int }_0^{+\infty }n{x}^{n-1}P\left(Q>x\right)\text{d}x}\\ =m_0^{\left(n\right)}+m_1^{\left(n\right)}\\ ={\displaystyle {\int }_0^{+\infty }n{x}^{n-1}{F}_0\left(x\right)\text{d}x}+{\displaystyle {\int }_0^{+\infty }n{x}^{n-1}{F}_1\left(x\right)\text{d}x}\\ ={\displaystyle \underset{i=0}{\overset{1}{\sum }}{\displaystyle {\int }_0^{+\infty }n{x}^{n-1}{\displaystyle \underset{k=0}{\overset{n_0}{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}{b}_i\left(k\right)\text{d}x}}}+{\displaystyle {\int }_0^{+\infty }n{x}^{n-1}e\left(n_0\right)\text{d}x}\end{array}$ (34)

由此，流模型稳态条件下的平均库存量为

$\begin{array}{c}E\left(Q\right)=E\left(Q^{\left(1\right)}\right)={\displaystyle {\int }_0^{+\infty }P\left(Q>x\right)\text{d}x}\\ ={\displaystyle {\int }_0^{+\infty }{F}_0\left(x\right)\text{d}x}+{\displaystyle {\int }_0^{+\infty }{F}_1\left(x\right)\text{d}x}+{\displaystyle {\int }_0^{+\infty }{F}_2\left(x\right)\text{d}x}\\ ={\displaystyle \underset{i=0}{\overset{1}{\sum }}{\displaystyle {\int }_0^{+\infty }{\displaystyle \underset{k=0}{\overset{n_0}{\sum }}{\text{e}}^{-\left(ax/d\right)}\frac{{\left(ax/d\right)}^k}{k!}{b}_i\left(k\right)\text{d}x}}}+{\displaystyle {\int }_0^{+\infty }e\left(n_0\right)\text{d}x}\\ =m_0^{\left(1\right)}+m_1^{\left(1\right)}\end{array}$ (35)

5. 数值例子

在本节中，进行了数值实验，来说明系统参数对系统性能指标的影响。

当 $\alpha =1$ ， $\beta =1.25$ ， $\eta =2$ 时，图1展示了当正常服务速率 $\mu$ 为8，9，10时平均库存量 $E\left(Q\right)$ 随顾客到达速率的变化趋势。可以发现，随着顾客到达速率的增大，系统平均库存量不断增加，但在相同达到率下，服务速率 $\mu$ 越大，平均库存量越小。

当 $\alpha =1$ ， $\beta =1.25$ ， $\lambda =2.5$ 时，图2展示了当服务速率 $\eta$ 为1.6，1.7，1.8时平均库存量 $E\left(Q\right)$ 随服务速率 $\mu$ 的变化趋势。可以发现，随着服务速率 $\mu$ 的增大，系统平均库存量不断减小并趋向于0，同时 $\eta$ 越大，对应的平均库存量越小。

当 $\alpha =1$ ， $\beta =1.25$ ， $\lambda =5$ ， $\mu =15$ 时，图3展示了当服务速率 $\eta$ 为4.0，4.2，4.4时库存量尾分布函数 $P\left\{Q_{\infty }>x\right\}$ 的变化趋势。可以发现，随着x的增大，库存量尾分布函数 $P\left\{Q_{\infty }>x\right\}$ 不断减小并趋向0，同时 $\eta$ 越大， $P\left\{Q_{\infty }>x\right\}$ 下降的速率越快。

Figure 1. $E\left(Q\right)$ versus $\mu$ and $\lambda$

图1. $E\left(Q\right)$ 随 $\mu$ 与 $\lambda$ 的变化曲线

Figure 2. $E\left(Q\right)$ versus $\mu$ and $\eta$

图2. $E\left(Q\right)$ 随 $\mu$ 与 $\eta$ 的变化曲线

Figure 3. $P\left\{Q_{\infty }>x\right\}$ versus $\eta$ and x

图3. $P\left\{Q_{\infty }>x\right\}$ 随 $\eta$ 与x的变化曲线

Table 1. Comparison of different truncation steps

表1. 不同截断步骤

当 $\alpha =1$ ， $\beta =1.25$ ， $\lambda =2.5$ ， $\mu =6$ ， $\eta =1$ 时，表1中展示了当给定误差限 $\epsilon =1\times {10}^{-9}$ 时，对于不同x的取值下获得的截断步骤 $n_{\epsilon }$ 与 $R\left(ax/d\right)$ 。不难发现，当系统参数固定时，截断步骤 $R\left(ax/d\right)$ 随着x的增加而增加，但是截断步骤 $n_{\epsilon }$ 是固定的。因此，当x较小时，可以选择 $R\left(ax/d\right)$ 作为计算尾分布函数截断步骤，但随着的x的增大，应该选择 $n_{\epsilon }$ 作为截断步骤。

基金项目

国家自然科学基金青年项目(11901307)。