具有两类故障的M/M/1 + 1排队系统的均衡分析

Equilibrium Analysis for an M/M/1 + 1 Queueing System with Two Types of Failures

摘要: 本文研究了完全可视环境下具有两类故障的M/M/1 + 1排队系统，其中主服务台采用N策略与工作休假机制，并设有备用服务台。主服务台可能发生两类故障，备用服务台在主服务台故障时辅助服务，且在服务期间也可能发生不完全故障。通过建立二维马尔可夫链模型，推导出系统最优进队阈值与稳态概率分布的显式表达式，进一步计算系统的平均队长、顾客平均逗留时间及社会收益等性能指标。通过数值模拟探讨服务率，故障率等关键参数对系统与社会收益的影响。仿真实验验证了解析结果的正确性。研究结果为高可靠性、可调节服务系统的运营优化与资源配置提供了理论依据与决策支持。

Abstract: This paper studies an M/M/1 + 1 queueing system with two types of failures in a fully observable environment. The primary server operates under an N policy with a working vacation mechanism, and a standby server is employed. The primary server is susceptible to two types of failures. The standby server assists the primary server during its failures and may also experience partial failures. By establishing a two-dimensional Markov chain model, explicit expressions for the optimal joining threshold and the steady-state probability distribution of the system are derived. Furthermore, performance metrics such as the average queue length, average customer sojourn time, and social welfare are calculated. Numerical simulations are conducted to investigate the impact of key parameters, including service rates and failure rates, on system behavior and social welfare. Simulation experiments validate the correctness of the analytical results. The findings provide theoretical foundations and decision support for the operational optimization and resource allocation of highly reliable and adjustable service systems.

1. 引言

服务台故障在实际运营中难以完全避免，随着设备的长期运行与老化，服务台可能出现不同形式的故障，进而影响系统的服务能力与顾客体验。Avi和Naor [1]最早对可能发生故障的排队系统展开研究；Economou和Spyridoula [2]以带故障与维修期的M/M/1系统为对象，分别在完全可视与几乎可视情形下讨论顾客的均衡策略。Zhang等[3]研究了带故障的M/M/1重试排队系统在完全可视和几乎不可视下的均衡策略和社会最优策略。Chang等[4]分析了具有顾客反馈与不耐烦行为的不可靠服务台重试排队系统，并比较了截断经典重试策略与恒定重试策略的性能差异。Yu等[5]进一步考虑部分故障与完全故障并存的M/M/1排队系统，给出了完全可视下的阈值策略以及完全不可视下的均衡策略。黎锁平等[6]研究了带启动时间和可修服务台的M/M/1/N单重工作休假排队系统，系统中服务台在休假时低速服务可以发生故障且故障后立即维修，运用矩阵几何法求稳态概率向量，进而求解系统方差、可用度、吞吐率等指标，还通过数值与敏感性分析验证方法的有效性。Lv [7]研究了由1个不可靠服务台与1个维修率变化的服务人员组成的多机可修复系统。唐韵和刘力维[8]考虑了带N策略、不可靠服务台和恒定重试率的M/M/1重试排队系统，得不同服务台状态下顾客均衡到达率，借助粒子群算法求得社会最优到达率，并做系统性能敏感性分析。

在许多应用场景中，服务台启动与保持运行会产生显著成本，因此引入按需启动的控制策略具有现实意义，N策略是其中一类典型方法：当系统内顾客数达到阈值N时，服务台启动，否则保持停歇以节约成本。Yadin和Naor [9]首次将N策略引入排队系统。杨顺利和田乃硕[10]研究了带有N策略工作休假的M/M/1排队系统，提出休假期间服务员以较低速率服务顾客的机制。Guo和Li [11]在完全可视与完全不可视两种信息结构下，研究了带N策略的M/M/1模型中顾客策略行为与社会最优问题。Vijayashree和Pavithra [12]进一步考虑运营故障对带N策略系统的影响，假设服务台在故障期内仍以较低速率运行，并采用母函数法求得稳态概率。另一方面，考虑到即使处于休假状态，系统也常需以较低速率提供必要服务，Servi和Finn [13]提出工作休假概念，并将其引入排队模型；Tian等[14]对M/M/1工作休假模型进行求解并给出随机分解结果。Li等[15]研究了含多工作假期与假期中断的单服务台M/M/1队列，在系统状态信息给定下分析顾客基于线性奖励–成本的加入或犹豫行为。Ma等[16]研究了带伯努利休假中断的单服务台马尔可夫工作休假排队系统，分析了到达顾客在四种不同信息水平下的策略性加入行为。杨喜娟等[17]研究了带启动时间、工作休假与工作故障的M/M/1/N可修排队系统。Bouchentouf等[18]聚焦于有限源多服务台系统，结合故障，修复与止步行为，并引入多同步工作休假机制。Yang等[19]研究了含准备时间、工作休假与工作故障的M/M/1/N排队系统，构建二维连续时间马尔可夫链与有限状态QBD过程，分析了参数对吞吐量与队列长度等指标的影响，并建立了成本优化模型。然而，上述研究多将故障模式，N策略与服务台的工作模式分别或两两结合进行探讨。在高度复杂的现实服务场景中，这些因素往往同时存在且相互耦合，而对此进行综合建模与分析的文献尚不多见。为此，本文旨在构建一个集成上述多重机制的综合模型进行研究。

基于上述研究背景，本文面向一类具有显著启动成本且难以避免随机故障的单服务台服务系统，针对其运行中常见的需求波动与可靠性退化问题构建排队模型。高负荷时段顾客集中到达容易引发拥堵与弃办，而低负荷时段持续满负荷运行又会带来能耗与资源浪费，因此有必要引入按需启停与降速服务等运行控制机制。为此，本文在主服务台的服务过程中引入N策略与工作休假机制，以刻画服务设施在不同负荷下的启停切换与低速服务行为。同时，考虑服务设施可能出现两类异质故障，一类为不完全故障导致服务能力下降，另一类为完全故障导致服务中断，并假设两类故障具有差异化的发生与修复过程。此外，为提升系统在故障期间的可用性，进一步引入备用服务台以实现辅助或替代服务，刻画备用服务台在特定运行阶段可能发生的次级故障。在此基础上，本文构建二维连续时间马尔科夫链并分析其稳态与性能表现，进而为此类可调节、高可靠服务系统的运行控制与资源配置提供定量参考与决策依据。

2. 模型描述

本文的模型描述如下：

(1) 顾客到达过程服从参数为$\lambda$ 的泊松过程。

(2) 系统配置一个主服务台，采用先到先服务(FCFS)规则逐个服务顾客。主服务台可进入工作休假状态且是不可靠的。当系统为空且主服务台完好时，系统进入工作休假状态，此时主服务台以低服务速率${\mu }_0$ 提供服务$\left({\mu }_0<\lambda \right)$ ；当系统中顾客数达到阈值N时，主服务台切换至正常工作状态，服务速率提升至${\mu }_1\left({\mu }_1>{\mu }_0\right)$ 。在正常工作状态下，主服务台可能发生两类故障：一类为不完全故障，使主服务台服务速率降为${\mu }_2\left({\mu }_2<{\mu }_1\right)$ ；另一类为完全故障，此时主服务台停止提供服务。主服务台两类故障的发生间隔以及相应修理时间分别服从参数为${\epsilon }_1、{\epsilon }_2、{\xi }_1、{\xi }_2$ 的指数分布。

(3) 当主服务台发生故障时，系统启用备用服务台。当主服务台发生不完全故障时，备用服务台以服务速率${\mu }_3$ 进行辅助服务$\left({\mu }_3<{\mu }_2\right)$ ，为了模拟备用服务台的性能局限性，假设其无法完全弥补主服务台的能力损失，即两者的合计服务速率满足${\mu }_2+{\mu }_3<{\mu }_1$ 。当主服务台发生完全故障时，备用服务台以服务速率${\mu }_4$ 单独提供服务。此外，在主服务台完全故障且由备用服务台独立服务期间，备用服务台可能发生不完全故障，其故障发生间隔与修理时间分别服从参数为${\epsilon }_3、{\xi }_3$ 的指数分布；发生不完全故障后，备用服务台服务速率降为${\mu }_5\left({\mu }_5>\lambda \right)$ 。而当主服务台不完全故障，备用服务台辅助服务时由于系统强度较低，此时备用服务台认为是完全可靠的，不会发生故障。

(4) 系统配置1名修理工，服务台的维修仅在系统为空时进行。当主服务台与备用服务台同时需要维修时，修理工优先修复主服务台。若主服务台已修复而备用服务台仍处于故障状态，则系统暂停对外服务，直至备用服务台修复完成。

(5) 假设主服务台的两类故障不会同时发生。当系统中仅有1名顾客且两台服务台均可工作时，优先由主服务台提供服务。此外，顾客到达间隔、服务时间，以及服务台故障发生间隔与维修时间相互独立。

令$\left\{\left(Q\left(t\right),J\left(t\right)\right)|t>0\right\}$ 为系统处于时刻t的二维状态向量，其中$Q\left(t\right)$ 为系统中的顾客数，$J\left(t\right)$ 为系统的运行状态，$J\left(t\right)$ 的具体含义如下：

$J\left(t\right)=\{\begin{array}{ll}0,系统处于工作休假期；\hfill & \hfill \\ 1,主服务台正常工作，备用服务台未启用；\hfill & \hfill \\ 2,主服务台不完全故障，备用服务台辅助工作；\hfill & \hfill \\ 3,主服务台完全故障，备用服务台辅助工作；\hfill & \hfill \\ 4,主服务台完全故障，备用服务台不完全故障；\hfill & \hfill \\ 5,主服务台已修理，备用服务台未修理好。\hfill & \hfill \end{array}$

显然$\left\{\left(Q\left(t\right),J\left(t\right)\right)|t>0\right\}$ 是一个连续时间的马尔科夫过程，其状态空间为

$\Omega =\left\{\left(0,0\right)\right\}\cup \left\{\left(i,0\right)|1\le i\le N-1\right\}\cup \left\{\left(i,1\right)|i\ge 1\right\}\cup \left\{\left(i,j\right)|i\ge 0,j=2,3,4\right\}\cup \left\{\left(0,5\right)\right\}.$

3. 均衡止步策略及稳态概率分布

3.1. 均衡止步策略

设顾客单位时间等待成本为$C_w$ ，服务完成可获得收益$R$ 。完全可视意味着顾客可以知晓系统中的顾客数与服务台状态的全部信息，为保证系统被激活，当系统处于休假状态其收益一定为正。系统具有进队阈值$n_e\left(1\right)$ ，当顾客到达时，若系统中的顾客数小于等于$n_e\left(1\right)$ 则进入系统，否则离开。

令$T\left(i,j\right)$ 为顾客平均逗留时间，此时该顾客处在系统的第$i$ 个位置且系统处于$j$ 状态。在完全可视下当所有顾客都遵循进队阈值$n_e\left(1\right)$ 时，其原始状态空间$\Omega$ 转变为有限状态空间

${\Omega }_f=\left\{\left(0,0\right)\right\}\cup \left\{\left(i,0\right)|1\le i\le N-1\right\}\cup \left\{\left(i,1\right)|1\le i\le n_e\left(1\right)+1\right\}\cup \left\{\left(i,j\right)|0\le i\le n_e\left(1\right)+1,j=2,3,4\right\}\cup \left\{\left(0,5\right)\right\}.$

定理1 对于上述带有$N$ 策略、工作休假、备用服务台和两类故障的可修排队系统，在完全可视下有且仅有止步阈值$n_e\left(1\right)=\lfloor x_e\rfloor$ ，$x_e$ 为下面方程的唯一解：

$\rho \left[1-{\left(\frac{{\mu }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\right)}^{x_e}\right]+\frac{x_e{\rho }_1}{{\epsilon }_1+{\epsilon }_2}+\frac{{\rho }_0}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}=\frac{R}{C_w},$

其中$\lfloor x\rfloor$ 为不超过$x$ 的最大整数。

证明：根据模型假设和全期望公式，平均逗留时间$T\left(i,1\right)$ 的表达式如下所示

$T\left(1,1\right)=\frac{1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}+\frac{{\epsilon }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\cdot \frac{1}{{\mu }_2}+\frac{{\epsilon }_2}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\left(\frac{1}{{\mu }_4+{\epsilon }_3}+\frac{{\epsilon }_3}{{\mu }_4+{\epsilon }_3}\cdot \frac{1}{{\mu }_5}\right);$ (1)

$T\left(i,1\right)=\frac{1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}+\frac{{\mu }_{1}T\left(i-1,1\right)}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}+\frac{{\epsilon }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\cdot \frac{i}{{\mu }_2+{\mu }_3}+\frac{{\epsilon }_2}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\left(\frac{i}{{\mu }_4+{\epsilon }_3}+\frac{{\epsilon }_3}{{\mu }_4+{\epsilon }_3}\cdot \frac{i}{{\mu }_5}\right),i\ge 2.$ (2)

令

${\rho }_0=1+\frac{{\epsilon }_1}{{\mu }_2}+\frac{{\epsilon }_2\left({\epsilon }_3+{\mu }_5\right)}{{\mu }_5\left({\epsilon }_3+{\mu }_4\right)},{\rho }_1=\frac{{\epsilon }_1}{{\mu }_2+{\mu }_3}+\frac{{\epsilon }_2\left({\epsilon }_3+{\mu }_5\right)}{{\mu }_5\left({\epsilon }_3+{\mu }_4\right)}.$

则$T\left(1,1\right)=\frac{{\rho }_0}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}$ ，迭代式(2)得

$T\left(i,1\right)=\rho \left[1-{\left(\frac{{\mu }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\right)}^{i-1}\right]+\frac{\left(i-1\right){\rho }_1}{{\epsilon }_1+{\epsilon }_2}+\frac{{\rho }_0}{{\mu }_1+{\epsilon }_1+{\epsilon }_2},i\ge 1.$ (3)

其中

$\rho =\frac{{\mu }_1{\rho }_0}{\left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right)\left({\epsilon }_1+{\epsilon }_2\right)}-\frac{{\mu }_1{\rho }_1}{{\left({\epsilon }_1+{\epsilon }_2\right)}^2}-\frac{{\epsilon }_1{\mu }_3}{\left({\epsilon }_1+{\epsilon }_2\right)\left({\mu }_2+{\mu }_3\right){\mu }_2}.$

在完全可视情形下，顾客进入系统时系统处于状态$\left(i,1\right)$ ，则该顾客接受服务后获得的收益为$S\left(i,1\right)=R-C_{w}T\left(i,1\right)$ ，将式(3)代入得

$S\left(i,1\right)=R-C_w\left[\rho \left(1-{\left(\frac{{\mu }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\right)}^i\right)+\frac{i{\rho }_1}{{\epsilon }_1+{\epsilon }_2}+\frac{{\rho }_0}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\right],i\ge 1.$

若要证$S\left(i,1\right)=0$ 有唯一解，只需证$T\left(i,1\right)$ 是单调递增的：

$T\left(i,1\right)-T\left(i-1,1\right)=\rho {\left(\frac{{\mu }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\right)}^{i-2}\left(1-\frac{{\mu }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}\right)+\frac{{\rho }_1}{{\epsilon }_1+{\epsilon }_2}\text{, }i\ge 2.$ (4)

由${\rho }_1=\frac{{\epsilon }_1}{{\mu }_2+{\mu }_3}+\frac{{\epsilon }_2\left({\epsilon }_3+{\mu }_5\right)}{{\mu }_5\left({\epsilon }_3+{\mu }_4\right)}>0$ 。

若$\rho \ge 0$ ，显然$T\left(i,1\right)$ 单调递增；若$\rho <0$ ，则$T\left(i,1\right)-T\left(i-1,1\right)$ 单调递增，只需证$T\left(2,1\right)-T\left(1,1\right)>0$ 。

将$i=2$ 代入式(4)整理得

$T\left(2,1\right)-T\left(1,1\right)=\frac{{\mu }_1{\rho }_0}{{\left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right)}^2}+\frac{{\rho }_1}{{\mu }_1+{\epsilon }_1+{\epsilon }_2}-\frac{{\epsilon }_1{\mu }_3}{\left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right)\left({\mu }_2+{\mu }_3\right){\mu }_2},$ (5)

由${\mu }_2>{\mu }_3$ ，可知

${\rho }_1>\frac{{\epsilon }_1{\mu }_3}{\left({\mu }_2+{\mu }_3\right){\mu }_2}.$

由此，便证明了$T\left(2,1\right)-T\left(1,1\right)>0$ ，即$T\left(i,1\right)$ 是单调递增的，所以方程$S\left(i,1\right)=0$ 一定存在唯一解，将其用$x_e$ 表示，方程两边同除$C_w$ ，移项得定理1，定理1证毕。

状态$J\left(t\right)=1$ 为主服务台正常工作阶段，服务速率最高，顾客等待时间最短，因此该状态下的进队阈值$n_e\left(1\right)$ 反映了系统在最佳服务条件下的顾客容纳上限。无论观察到系统处于何种状态，只要队列长度$i>n_e\left(1\right)$ ，则进入系统的预期收益，顾客的最优选择是止步。为构建一个可分析的纯阈值均衡策略，我们假设当队列长度$i\le n_e\left(1\right)$ 时，顾客选择进入系统。这一假设在系统可靠性高，主服务台正常工作状态占主导的运营场景下是合理的。

3.2. 稳态概率分析

在完全可视情形下，当所有顾客都遵循进队阈值$n_e\left(1\right)$ 时，系统的状态空间由$\Omega$ 收缩为有限集合${\Omega }_f$ ，并且该连续时间的马尔可夫链是不可约的，因而稳态分布存在且唯一。记

令$\pi =\left\{{\pi }_{i,j}:\left(i,j\right)\in {\Omega }_f\right\}$ 表示系统的稳态分布向量。

定理2 对于上述带有$N$ 策略，工作休假且含两类故障的M/M/1 + 1复合排队系统，当所有顾客都遵循进队阈值$n_e\left(1\right)$ 时，其稳态分布$\pi$ 满足如下显式表达式：

${\pi }_{n,1}=\{\begin{array}{ll}\left(x_3^n-x_4^n\right)B_1,\hfill & 1\le n\le N;\hfill \\ b\left(a{x}_3^n+x_4^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$

${\pi }_{n,3}=\{\begin{array}{ll}\left(E_1{x}_5^n+E_2{x}_6^n+k_1{x}_3^n+k_2{x}_4^n\right)B_1,\hfill & 0\le n\le N;\hfill \\ \left(E_3{x}_5^n+E_4{x}_6^n+ab{k}_1{x}_3^n-b{k}_2{x}_4^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$

${\pi }_{0,5}=\frac{{\xi }_2}{{\xi }_3}\psi B_1.$

证明：在完全可视下当所有顾客都遵循进队阈值$n_e\left(1\right)$ 时，系统的状态转移图为(图1)：

Figure 1. State transition diagram of the system

图1. 系统状态转移图

由状态转移图可得平衡方程为：

$\lambda {\pi }_{0,0}={\mu }_0{\pi }_{1,0}+{\mu }_1{\pi }_{1,1}+{\xi }_1{\pi }_{0,2}+{\xi }_2{\pi }_{0,3}+{\xi }_3{\pi }_{0,5}.$ (6)

$\left(\lambda +{\mu }_0\right){\pi }_{n,0}=\lambda {\pi }_{n-1,0}+{\mu }_0{\pi }_{n+1,0},1\le n\le N-2.$ (7)

$\left(\lambda +{\mu }_0\right){\pi }_{N-1,0}=\lambda {\pi }_{N-2,0}.$ (8)

$\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right){\pi }_{1,1}={\mu }_1{\pi }_{2,1}.$ (9)

$\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right){\pi }_{n,1}=\lambda {\pi }_{n-1,1}+{\mu }_1{\pi }_{n+1,1},2\le n\le N-1.$ (10)

$\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right){\pi }_{N,1}=\lambda {\pi }_{N-1,1}+{\mu }_1{\pi }_{N+1,1}+\lambda {\pi }_{N-1,0}.$ (11)

$\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right){\pi }_{n,1}=\lambda {\pi }_{n-1,1}+{\mu }_1{\pi }_{n+1,1},N+1\le n\le n_e\left(1\right).$ (12)

$\left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right){\pi }_{n_e\left(1\right)+1,1}=\lambda {\pi }_{n_e\left(1\right),1}.$ (13)

$\left(\lambda +{\xi }_1\right){\pi }_{0,2}={\mu }_2{\pi }_{1,2}.$ (14)

$\left(\lambda +{\mu }_2\right){\pi }_{1,2}=\lambda {\pi }_{0,2}+\left({\mu }_2+{\mu }_3\right){\pi }_{2,2}+{\epsilon }_1{\pi }_{1,1}.$ (15)

$\left(\lambda +{\mu }_2+{\mu }_3\right){\pi }_{n,2}=\lambda {\pi }_{n-1,2}+\left({\mu }_2+{\mu }_3\right){\pi }_{n+1,2}+{\epsilon }_1{\pi }_{n,1},2\le n\le n_e\left(1\right).$ (16)

$\left({\mu }_2+{\mu }_3\right){\pi }_{n_e\left(1\right)+1,2}=\lambda {\pi }_{n_e\left(1\right),2}+{\epsilon }_1{\pi }_{n_e\left(1\right)+1,1}.$ (17)

$\left(\lambda +{\xi }_2\right){\pi }_{0,3}={\mu }_4{\pi }_{1,3}.$ (18)

$\left(\lambda +{\epsilon }_3+{\mu }_4\right){\pi }_{n,3}=\lambda {\pi }_{n-1,3}+{\epsilon }_2{\pi }_{n,1}+{\mu }_4{\pi }_{n+1,3},1\le n\le n_e\left(1\right).$ (19)

$\left({\epsilon }_3+{\mu }_4\right){\pi }_{n_e\left(1\right)+1,3}=\lambda {\pi }_{n_e\left(1\right),3}+{\epsilon }_2{\pi }_{n_e\left(1\right)+1,1}.$ (20)

$\left(\lambda +{\xi }_2\right){\pi }_{0,4}={\mu }_5{\pi }_{1,4}.$ (21)

$\left(\lambda +{\mu }_5\right){\pi }_{n,4}=\lambda {\pi }_{n-1,4}+{\epsilon }_3{\pi }_{n,3}+{\mu }_5{\pi }_{n+1,4},1\le n\le n_e\left(1\right).$ (22)

${\mu }_5{\pi }_{n_e\left(1\right)+1,4}=\lambda {\pi }_{n_e\left(1\right),4}+{\epsilon }_3{\pi }_{n_e\left(1\right)+1,3}.$ (23)

${\xi }_3{\pi }_{0,5}={\xi }_2{\pi }_{0,4}.$ (24)

由式(7)为常系数线性齐次差分方程，特征方程为${\mu }_0{x}^2-\left(\lambda +{\mu }_0\right)x+\lambda =0$ ，特征值为$x_1=1$ 、$x_2=\frac{\lambda }{{\mu }_0}$ ，所以${\pi }_{n,0}=A_1+A_2{\left(\frac{\lambda }{{\mu }_0}\right)}^n，0\le n\le N-1$ ，其中$A_1$ 和$A_2$ 为待定系数。

将$n=N-1、N-2$ 代入式(8)，得方程组

$\{\begin{array}{l}{A}_1+A_2={\pi }_{0,0},\hfill \\ A_1+A_2{\left(\frac{\lambda }{{\mu }_0}\right)}^{N-1}=\frac{\lambda }{\lambda +{\mu }_0}\left(A_1+A_2{\left(\frac{\lambda }{{\mu }_0}\right)}^{N-2}\right).\hfill \end{array}$

解得

$A_1=\frac{{\pi }_{0,0}}{1-{\left(\frac{{\mu }_0}{\lambda }\right)}^N},A_2=\frac{{\left(\frac{{\mu }_0}{\lambda }\right)}^N{\pi }_{0,0}}{{\left(\frac{{\mu }_0}{\lambda }\right)}^N-1}.$

所以

${\pi }_{n,0}=\frac{1-{\left(\frac{{\mu }_0}{\lambda }\right)}^{N-n}}{1-{\left(\frac{{\mu }_0}{\lambda }\right)}^N}{\pi }_{0,0},0\le n\le N-1.$ (25)

同理，式(10)特征方程为${\mu }_1{x}^2-\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right)x+\lambda =0$ 特征值为

$x_3=\frac{\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2+\sqrt{{\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right)}^2-4\lambda {\mu }_1}}{2{\mu }_1}$ ，$x_4=\frac{\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2-\sqrt{{\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right)}^2-4\lambda {\mu }_1}}{2{\mu }_1}.$

所以

${\pi }_{n,1}=B_1{x}_3^n+B_2{x}_4^n,1\le n\le N.\left(B_1、B_2为待定系数\right)$ (26)

将上式(26)代入式(9)，得方程组

$\{\begin{array}{l}{B}_1{x}_3+B_2{x}_4={\pi }_{1,1},\hfill \\ B_1{x}_3^2+B_2{x}_4^2=\frac{\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2}{{\mu }_1}{\pi }_{1,1}.\hfill \end{array}$

解得

$B_1=\frac{{\mu }_1}{\sqrt{{\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right)}^2-4\lambda {\mu }_1}}{\pi }_{1,1},B_2=-\frac{{\mu }_1}{\sqrt{{\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right)}^2-4\lambda {\mu }_1}}{\pi }_{1,1}.$

所以

${\pi }_{n,1}=\left(x_3^n-x_4^n\right)B_1,1\le n\le N.$

当$N+1\le n\le n_e\left(1\right)$ 时，同理得到

${\pi }_{n,1}=D_1{x}_3^n+D_2{x}_4^n,N\le n\le n_e\left(1\right)+1.\left(D_1、D_2为待定系数\right)$

由临界值$n=N$ 与式(13)得方程组

$\{\begin{array}{l}{D}_1{x}_3^N+D_2{x}_4^N=\left(x_3^N-x_4^N\right)B_1,\hfill \\ \left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right)\left(D_1{x}_3^{n_e\left(1\right)+1}+D_2{x}_4^{n_e\left(1\right)+1}\right)=\lambda \left(D_1{x}_3^{n_e\left(1\right)}+D_2{x}_4^{n_e\left(1\right)}\right).\hfill \end{array}$

令

$a=\frac{\lambda -\left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right)x_4}{\left({\mu }_1+{\epsilon }_1+{\epsilon }_2\right)x_3-\lambda }{\left(\frac{x_4}{x_3}\right)}^{n_e\left(1\right)},b=\frac{x_3^N-x_4^N}{a{x}_3^N+x_4^N}.$

解得$D_2=b{B}_1$ 、$D_1=ab{B}_1$ ，所以

${\pi }_{n,1}=\{\begin{array}{ll}\left(x_3^n-x_4^n\right)B_1,\hfill & 1\le n\le N,\hfill \\ b\left(a{x}_3^n+x_4^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$ (27)

将式(25)和(27)代入式(11)得到

${\pi }_{0,0}=\frac{1-{\left(\frac{{\mu }_0}{\lambda }\right)}^N}{\lambda -{\mu }_0}c{B}_1,$ (28)

其中

$c=\left(\lambda +{\mu }_1+{\epsilon }_1+{\epsilon }_2\right)\left(x_3^N-x_4^N\right)-\lambda \left(x_3^{N-1}-x_4^{N-1}\right)-{\mu }_{1}b\left(a{x}_3^{N+1}+x_4^{N+1}\right).$

将式(28)代入(25)得

${\pi }_{n,0}=\frac{1-{\left(\frac{{\mu }_0}{\lambda }\right)}^{N-n}}{\lambda -{\mu }_0}c{B}_1,0\le n\le N-1.$ (29)

令${\mu }_2^{*}={\mu }_2+{\mu }_3$ ，式(16)变为

${\pi }_{n,2}-{\pi }_{2,2}=\frac{\lambda }{{\mu }_2^{*}}\left({\pi }_{n-1,2}-{\pi }_{1,2}\right)-\frac{{\epsilon }_1}{{\mu }_2^{*}}\left({\pi }_{2,1}+\cdots +{\pi }_{n-1,1}\right),3\le n\le n_e\left(1\right)+1.$ (30)

由式(14)和(15)得

${\pi }_{1,2}=\frac{\lambda +{\xi }_1}{{\mu }_2}{\pi }_{0,2}{\pi }_{2,2}=\frac{1}{{\mu }_2^{*}}\left[\frac{\left(\lambda +{\mu }_2\right)\left(\lambda +{\xi }_1\right)}{{\mu }_2}-\lambda \right]{\pi }_{0,2}-\frac{{\epsilon }_1}{{\mu }_2^{*}}{\pi }_{1,1},$

将其代入式(30)得

${\pi }_{n,2}=\frac{\lambda }{{\mu }_2^{*}}{\pi }_{n-1,2}+\frac{{\xi }_1}{{\mu }_2^{*}}{\pi }_{0,2}-\frac{{\epsilon }_1}{{\mu }_2^{*}}\left({\pi }_{1,1}+\cdots +{\pi }_{n-1,1}\right),3\le n\le n_e\left(1\right).$ (31)

将式(17)与(27)代入上式(31)整理得

${\pi }_{0,2}=d{B}_1,$

其中

$d=\frac{{\epsilon }_1}{{\xi }_1}\left(z_1+z_2\right),z_1=\frac{x_3\left(1-x_3^N\right)}{1-x_3}-\frac{x_4\left(1-x_4^N\right)}{1-x_4}$

令$M_n=\frac{{\displaystyle \sum _{i=1}^n{\pi }_{i,1}}}{B_1}$ ，迭代上式(31)得，当$2\le n\le n_e\left(1\right)+1$ 时

${\pi }_{n,2}=\left[\frac{{\xi }_1}{{\mu }_2^{*}-\lambda }+\left(\frac{\lambda +{\xi }_1}{{\mu }_2}-\frac{{\xi }_1}{{\mu }_2^{*}-\lambda }\right){\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n-1}\right]d{B}_1-\frac{{\epsilon }_1}{{\mu }_2^{*}}{\displaystyle \sum _{j=1}^{n-1}{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n-1-j}}{M}_j{B}_1.$ (32)

由状态1的表达式(27)得

$M_n=\{\begin{array}{ll}\frac{x_3\left(1-x_3^n\right)}{1-x_3}-\frac{x_4\left(1-x_4^n\right)}{1-x_4},\hfill & 1\le n\le N,\hfill \\ M_N+b\left(\frac{a{x}_3^{N+1}\left(1-x_3^{n-N}\right)}{1-x_3}+\frac{x_4^{N+1}\left(1-x_4^{n-N}\right)}{1-x_4}\right),\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$

令

$\begin{array}{l}{f}_1\left(n\right)=\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n-1}\right)\frac{x_3-x_4}{\left({\mu }_2^{*}-\lambda \right)\left(1-x_3\right)\left(1-x_4\right)}-\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}{x}_3}\right)}^{n-1}\right)\frac{x_3^{n+1}}{\left(1-x_3\right)\left({\mu }_2^{*}{x}_3-\lambda \right)}\\ +\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}{x}_4}\right)}^{n-1}\right)\frac{x_4^{n+1}}{\left(1-x_4\right)\left({\mu }_2^{*}{x}_4-\lambda \right)},1\le n\le N+1,\end{array}$

$\begin{array}{l}{f}_2\left(n\right)=\frac{1}{{\mu }_2^{*}-\lambda }\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n-N-1}\right)\left(\frac{a{x}_3^{N+1}}{1-x_3}+\frac{x_4^{N+1}}{1-x_4}\right)-\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}{x}_4}\right)}^{n-N-1}\right)\frac{x_4^{n+1}}{\left(1-x_4\right)\left({\mu }_2^{*}{x}_4-\lambda \right)}\\ -\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}{x}_3}\right)}^{n-N-1}\right)\frac{a{x}_3^{n+1}}{\left(1-x_3\right)\left({\mu }_2^{*}{x}_3-\lambda \right)},N+2\le n\le n_e\left(1\right)+1.\end{array}$

将$M_n$ 的表达式代入式(32)，整理得

(33)

式(19)的特征方程为$x^2-\frac{\lambda +{\epsilon }_3+{\mu }_4}{{\mu }_4}x+\frac{\lambda }{{\mu }_4}=0$ ，令

$e_1=\frac{\lambda +{\epsilon }_3+{\mu }_4}{{\mu }_4},e_2=\frac{\lambda }{{\mu }_4}.$

特征值为

$x_5=\frac{e_1+\sqrt{e_1^2-4e_2}}{2},x_6=\frac{e_1-\sqrt{e_1^2-4e_2}}{2}.$

设齐次方程对应的通解为

${\pi }_{n,3}^{hom}=\{\begin{array}{ll}\left(E_1{x}_5^n+E_2{x}_6^n\right)B_1,\hfill & 0\le n\le N,\hfill \\ \left(E_3{x}_5^n+E_4{x}_6^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$ (34)

由于

${\pi }_{n,1}=\{\begin{array}{ll}\left(x_3^n-x_4^n\right)B_1,\hfill & 1\le n\le N,\hfill \\ b\left(a{x}_3^n+x_4^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$

所以设特解形式为

${\pi }_{n,3}^{spe}=\{\begin{array}{ll}\left(k_1{x}_3^n+k_2{x}_4^n\right)B_1,\hfill & 0\le n\le N,\hfill \\ \left(k_3{x}_3^n+k_4{x}_4^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$ (35)

令$e_3=\frac{{\epsilon }_2}{{\mu }_4}$ ，将上式(35)与${\pi }_{n,1}$ 通项(27)代入(19)，得当$1\le n\le N$ 时

$\{\begin{array}{l}{k}_1\left(x_3^2-e_1{x}_3+e_2\right)=-e_3{x}_3,\hfill \\ k_2\left(x_4^2-e_1{x}_4+e_2\right)=e_3{x}_4.\hfill \end{array}\Rightarrow \{\begin{array}{l}{k}_1=\frac{-e_3{x}_3}{x_3^2-e_1{x}_3+e_2},\hfill \\ k_2=\frac{e_3{x}_4}{x_4^2-e_1{x}_4+e_2}.\hfill \end{array}$

当$N+1\le n\le n_e\left(1\right)+1$ 时

$\{\begin{array}{l}{k}_3\left(x_3^2-e_1{x}_3+e_2\right)=-ab{e}_3{x}_3,\hfill \\ k_4\left(x_4^2-e_1{x}_4+e_2\right)=-b{e}_3{x}_4.\hfill \end{array}\Rightarrow \{\begin{array}{l}{k}_3=\frac{-ab{e}_3{x}_3}{x_3^2-e_1{x}_3+e_2}=ab{k}_1,\hfill \\ k_4=\frac{-b{e}_3{x}_4}{x_4^2-e_1{x}_4+e_2}=-b{k}_2.\hfill \end{array}$

所以状态3的通解${\pi }_{n,3}={\pi }_{n,3}^{hom}+{\pi }_{n,3}^{spe}$ ，其表达式为

${\pi }_{n,3}=\{\begin{array}{ll}\left(E_1{x}_5^n+E_2{x}_6^n+k_1{x}_3^n+k_2{x}_4^n\right)B_1,\hfill & 0\le n\le N,\hfill \\ \left(E_3{x}_5^n+E_4{x}_6^n+ab{k}_1{x}_3^n-b{k}_2{x}_4^n\right)B_1,\hfill & N+1\le n\le n_e\left(1\right)+1.\hfill \end{array}$ (36)

将式(27)、(36)代入方程(18)、(19)和(20)，整理得线性方程组

$\left[\begin{array}{cccc}{A}_{11}& A_{12}& 0& 0\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ 0& 0& A_{43}& A_{44}\end{array}\right]\left[\begin{array}{c}{E}_1\\ E_2\\ E_3\\ E_4\end{array}\right]=\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\end{array}\right]$ (37)

其中

${\alpha }_1=\lambda +{\xi }_2,{\alpha }_2={\mu }_4,{\alpha }_3=\lambda +{\epsilon }_3+{\mu }_4,{\alpha }_4={\epsilon }_3+{\mu }_4,$

$A_{11}={\alpha }_1-{\alpha }_2{x}_5,A_{12}={\alpha }_1-{\alpha }_2{x}_6,A_{21}={\alpha }_3{x}_5^N-\lambda x_5^{N-1},A_{22}={\alpha }_3{x}_6^N-\lambda x_6^{N-1},A_{23}=-{\alpha }_2{x}_5^{N+1},$

$A_{24}=-{\alpha }_2{x}_6^{N+1},A_{31}=-\lambda x_5^N,A_{32}=-\lambda x_6^N,A_{33}={\alpha }_3{x}_5^{N+1}-{\alpha }_2{x}_5^{N+2},A_{34}={\alpha }_3{x}_6^{N+1}-{\alpha }_2{x}_6^{N+2},$

$A_{43}={\alpha }_4{x}_5^{n_e(1)+1}-\lambda x_5^{n_e(1)},A_{44}={\alpha }_4{x}_6^{n_e(1)+1}-\lambda x_6^{n_e(1)},$

$F_1={\alpha }_2\left(k_1{x}_3+k_2{x}_4\right)-{\alpha }_1\left(k_1+k_2\right),$

$F_2=\lambda \left(k_1{x}_3^{N-1}+k_2{x}_4^{N-1}\right)+{\epsilon }_2\left(x_3^N-x_4^N\right)+{\alpha }_2\left(ab{k}_1{x}_3^{N+1}-b{k}_2{x}_4^{N+1}\right)-{\alpha }_3\left(k_1{x}_3^N+k_2{x}_4^N\right),$ $F_3=\lambda \left(k_1{x}_3^N+k_2{x}_4^N\right)+{\epsilon }_{2}b\left(a{x}_3^{N+1}+x_4^{N+1}\right)+{\alpha }_2\left(ab{k}_1{x}_3^{N+2}-b{k}_2{x}_4^{N+2}\right)-{\alpha }_3\left(ab{k}_1{x}_3^{N+1}-b{k}_2{x}_4^{N+1}\right),$

$F_4=\lambda \left(ab{k}_1{x}_3^{n_e(1)}-b{k}_2{x}_4^{n_e(1)}\right)+{\epsilon }_{2}b\left(a{x}_3^{n_e(1)+1}+x_4^{n_e(1)+1}\right)-{\alpha }_4\left(ab{k}_1{x}_3^{n_e(1)+1}-b{k}_2{x}_4^{n_e(1)+1}\right).$

若系统稳定存在则$E_1、E_2、E_3、E_4$ 必存在且唯一，即系数矩阵非奇异：

$\left[\begin{array}{c}{E}_1\\ E_2\\ E_3\\ E_4\end{array}\right]={\left[\begin{array}{cccc}{A}_{11}& A_{12}& 0& 0\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ 0& 0& A_{43}& A_{44}\end{array}\right]}^{-1}\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\end{array}\right]$

将式(27)、(29)、(33)、(36)和式(24)代入(6)得

${\pi }_{0,4}=\psi B_1,$

其中

$\psi =\frac{1}{{\xi }_2}\left(c-{\mu }_1\left(x_3-x_4\right)-{\xi }_{1}d-{\xi }_2\left(E_1+E_2+k_1+k_2\right)\right).$

类似于状态2，迭代式(22)得

${\pi }_{n,4}-{\pi }_{1,4}=\frac{\lambda }{{\mu }_5}\left({\pi }_{n-1,4}-{\pi }_{0,4}\right)-\frac{{\epsilon }_3}{{\mu }_5}\left({\pi }_{1,3}+\cdots +{\pi }_{n-1,3}\right),2\le n\le n_e\left(1\right)+1.$ (38)

令$G_n=\frac{{\displaystyle \sum _{i=1}^n{\pi }_{i,3}}}{B_1}$ ，并将${\pi }_{0,4}$ 和式(21)代入上式(38)得

${\pi }_{n,4}-\frac{\lambda }{{\mu }_5}{\pi }_{n-1,4}=\frac{{\xi }_2}{{\mu }_5}\psi B_1-\frac{{\epsilon }_3}{{\mu }_5}{G}_{n-1}{B}_1,2\le n\le n_e\left(1\right)+1.$ (39)

迭代上式(39)得

${\pi }_{n,4}=\left[\frac{{\xi }_2}{{\mu }_5-\lambda }+\left(1-\frac{{\xi }_2}{{\mu }_5-\lambda }\right){\left(\frac{\lambda }{{\mu }_5}\right)}^n\right]\psi B_1-\frac{{\epsilon }_3}{{\mu }_5}{\displaystyle \sum _{i=1}^{n-1}{\left(\frac{\lambda }{{\mu }_5}\right)}^{n-i-1}}{G}_i{B}_1,1\le n\le n_e\left(1\right)+1.$ (40)

由状态3的通项式(36)得

$G_i=\{\begin{array}{l}\frac{D_1{x}_5\left(1-x_5^i\right)}{1-x_5}+\frac{D_2{x}_6\left(1-x_6^i\right)}{1-x_6}+\frac{k_1{x}_3\left(1-x_3^i\right)}{1-x_3}+\frac{k_2{x}_4\left(1-x_4^i\right)}{1-x_4},1\le i\le N,\\ G_N+\frac{D_3{x}_5^{N+1}\left(1-x_5^{i-N}\right)}{1-x_5}+\frac{D_4{x}_6^{N+1}\left(1-x_6^{i-N}\right)}{1-x_6}+\frac{ab{k}_1{x}_3^{N+1}\left(1-x_3^{i-N}\right)}{1-x_3}-\frac{b{k}_2{x}_4^{N+1}\left(1-x_4^{i-N}\right)}{1-x_4},N+1\le i\le n_e\left(1\right)+1.\end{array}$

令

$q_2\left(n\right)=\left[\frac{{\xi }_2}{{\mu }_5-\lambda }+\left(1-\frac{{\xi }_2}{{\mu }_5-\lambda }\right){\left(\frac{\lambda }{{\mu }_5}\right)}^n\right],$

将$G_i$ 的表达式代入式(40)，整理得

最后，利用正规化条件求解$B_1$ ：

${\displaystyle \sum _{i=0}^{N-1}{\pi }_{i,0}}+{\displaystyle \sum _{i=1}^{n_e\left(1\right)+1}{\pi }_{i,1}}+{\displaystyle \sum _{i=0}^{n_e\left(1\right)+1}{\displaystyle \sum _{j=2}^4{\pi }_{i,j}}}+{\pi }_{0,5}=1.$

令

${\Phi }_3\left(x,\mu ,n\right)=\frac{x^2}{\left(\lambda -\mu x\right)\left(1-x\right)}\left[\frac{\lambda }{\mu -\lambda }\left(1-{\left(\frac{\lambda }{\mu }\right)}^n\right)-\frac{x\left(1-x^n\right)}{1-x}\right],$ ${\delta }_1=\left[1+\frac{n_e\left(1\right)+1}{{\mu }_2^{*}-\lambda }{\xi }_1+\left(\frac{\lambda +{\xi }_1}{{\mu }_2\left({\mu }_2^{*}-\lambda \right)}-\frac{{\mu }_2^{*}{\xi }_1}{{\left({\mu }_2^{*}-\lambda \right)}^2}\right)\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n_e\left(1\right)+1}\right)\right]d,$

${\delta }_2=\left[\frac{n_e\left(1\right)+2}{{\mu }_5-\lambda }{\xi }_2+\left(\frac{{\mu }_5}{{\mu }_5-\lambda }-\frac{{\mu }_5{\xi }_2}{{\left({\mu }_5-\lambda \right)}^2}\right)\left(1-{\left(\frac{\lambda }{{\mu }_5}\right)}^{n_e\left(1\right)+2}\right)\right]\psi ,$

$\begin{array}{c}\gamma =\frac{E_1\left(1-x_5^{N+1}\right)}{1-x_5}+\frac{E_2\left(1-x_6^{N+1}\right)}{1-x_6}+\frac{k_1\left(1-x_3^{N+1}\right)}{1-x_3}+\frac{k_2\left(1-x_4^{N+1}\right)}{1-x_4}+\frac{E_3{x}_5^{N+1}\left(1-x_5^{n_e\left(1\right)-N+1}\right)}{1-x_5}\\ +\frac{E_4{x}_6^{N+1}\left(1-x_6^{n_e\left(1\right)-N+1}\right)}{1-x_6}+\frac{ab{k}_1{x}_3^{N+1}\left(1-x_3^{n_e\left(1\right)-N+1}\right)}{1-x_3}-\frac{b{k}_2{x}_4^{N+1}\left(1-x_4^{n_e\left(1\right)-N+1}\right)}{1-x_4},\end{array}$

$S_1=\frac{x_3-x_4}{\left({\mu }_2^{*}-\lambda \right)\left(1-x_3\right)\left(1-x_4\right)}\left[N-\frac{\lambda }{{\mu }_2^{*}-\lambda }\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^N\right)\right]-{\Phi }_3\left(x_3,{\mu }_2^{*},N\right)+{\Phi }_3\left(x_4,{\mu }_2^{*},N\right),$

${\tau }_1=\frac{\lambda }{{\mu }_2^{*}-\lambda }\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n_e\left(1\right)-N}\right)f_1\left(N+1\right),{\tau }_2=\frac{M_N}{{\mu }_2^{*}-\lambda }\left(n_e\left(1\right)-N-\frac{\lambda }{{\mu }_2^{*}-\lambda }\left(1-{\left(\frac{\lambda }{{\mu }_2^{*}}\right)}^{n_e\left(1\right)-N}\right)\right),$

${\tau }_3=\frac{\lambda }{{\mu }_5-\lambda }\left(1-{\left(\frac{\lambda }{{\mu }_5}\right)}^{n_e\left(1\right)-N}\right)t_1\left(N+1\right),{\tau }_4=\frac{G_N}{{\mu }_5-\lambda }\left(n_e\left(1\right)-N-\frac{\lambda }{{\mu }_5-\lambda }\left(1-{\left(\frac{\lambda }{{\mu }_5}\right)}^{n_e\left(1\right)-N}\right)\right),$

$\begin{array}{l}{\beta }_1=\left(\frac{E_1{x}_5}{1-x_5}+\frac{E_2{x}_6}{1-x_6}+\frac{k_1{x}_3}{1-x_3}+\frac{k_2{x}_4}{1-x_4}\right)\left[\frac{N}{{\mu }_5-\lambda }-\frac{\lambda }{{\left({\mu }_5-\lambda \right)}^2}\left(1-{\left(\frac{\lambda }{{\mu }_5}\right)}^N\right)\right]\\ -E_1{\Phi }_3\left(x_5,{\mu }_5,N\right)-E_2{\Phi }_3\left(x_6,{\mu }_5,N\right)-k_1{\Phi }_3\left(x_3,{\mu }_5,N\right)-k_2{\Phi }_3\left(x_4,{\mu }_5,N\right),\end{array}$

将各状态稳态概率依次求和，整理得

${\displaystyle \sum _{i=0}^{N-1}{\pi }_{i,0}}=\left[\frac{N}{\lambda -{\mu }_0}-\frac{{\mu }_0}{{\left(\lambda -{\mu }_0\right)}^2}\left(1-{\left(\frac{{\mu }_0}{\lambda }\right)}^N\right)\right]c{B}_1,{\displaystyle \sum _{i=0}^{n_e\left(1\right)+1}{\pi }_{i,1}}=\left(z_1+z_2\right)B_1,$ ${\displaystyle \sum _{i=0}^{n_e\left(1\right)+1}{\pi }_{i,2}}=\left[{\delta }_1-{\epsilon }_1\left(S_1+b{S}_2+{\tau }_1+{\tau }_2\right)\right]B_1\text{, }{\displaystyle \sum _{i=0}^{n_e\left(1\right)+1}{\pi }_{i,3}}=\gamma B_1,$

${\displaystyle \sum _{i=0}^{n_e\left(1\right)+1}{\pi }_{i,4}}=\left[{\delta }_2-{\epsilon }_3\left({\beta }_1+{\beta }_2+{\tau }_3+{\tau }_4\right)\right]B_1\text{, }{\pi }_{0,5}=\frac{{\xi }_2}{{\xi }_3}\psi B_1.$

$B_1={\left[\left(\frac{N}{\lambda -{\mu }_0}-\frac{{\mu }_0}{{\left(\lambda -{\mu }_0\right)}^2}\left(1-{\left(\frac{{\mu }_0}{\lambda }\right)}^N\right)\right)c+z_1+z_2+{\delta }_1+{\delta }_2+\gamma -{\epsilon }_1\left(S_1+b{S}_2+{\tau }_1+{\tau }_2\right)-{\epsilon }_3\left({\beta }_1+{\beta }_2+{\tau }_3+{\tau }_4\right)+\frac{{\xi }_2}{{\xi }_3}\psi \right]}^{-1}.$

3.3. 特例分析

为了验证本文所得稳态概率分布的正确性，并揭示模型与经典排队系统的内在联系，本节考虑退化情形：通过设置部分参数，可验证本文的稳态概率公式退化为经典的M/M/1/K模型。

不考虑设备故障以及N策略，$n_e\left(1\right)+1$ 取K，状态空间收缩为${\Omega }_f=\left\{i|1\le i\le K\right\}$ 。为保证系统稳定假

设$\rho =\frac{\lambda }{{\mu }_1}<1$ ，并将$\left(N,n_e\left(1\right),{\mu }_0,{\epsilon }_1,{\epsilon }_2,{\epsilon }_3\right)=\left(1,K-1,0,0,0,0\right)$ 代入定理2得

$x_3=1,x_4=\rho =\frac{\lambda }{{\mu }_1}.$

将$x_3、x_4$ 代入$a、b$ 的表达式得

$a=0,b=\frac{1-\rho }{\rho }.$

所以

${\pi }_n=\{\begin{array}{ll}\left(1-\rho \right)B_1,\hfill & n=1,\hfill \\ \left(1-\rho \right){\rho }^{n-1}{B}_1,\hfill & 2\le n\le K.\hfill \end{array}$

同理

$c={\mu }_1-\lambda .$

代入到定理2${\pi }_{n,0}$ 通项中得到

${\pi }_0=\frac{1-\rho }{\rho }{B}_1$

最后利用正规化条件求解$B_1$ ：

${\pi }_0+{\displaystyle \sum _{i=1}^K{\pi }_i}=1.$

解得

$B_1=\frac{\rho }{1-{\rho }^{K+1}}.$

整理得

${\pi }_n=\frac{\left(1-\rho \right){\rho }^n}{1-{\rho }^{K+1}},\text{ 0}\le n\le K.$

由此可见，本文的解析结果在$\left(N,n_e\left(1\right),{\mu }_0,{\epsilon }_1,{\epsilon }_2,{\epsilon }_3\right)=\left(1,K-1,0,0,0,0\right)$ 的退化条件下，能够严格导出经典M/M/1/K队列的稳态概率分布，从而验证了理论推导的正确性与一般性。

4. 稳性能指标与数值分析

4.1. 性能指标

(1) 系统平均队长为

$E\left(L\right)={\displaystyle \sum _{n=1}^{N-1}n}{\pi }_{n,0}+{\displaystyle \sum _{n=1}^{n_e\left(1\right)+1}{\displaystyle \sum _{j=1}^{4}n}}{\pi }_{n,j}.$

(2) 系统平均等待队长为

$E\left(L_q\right)={\displaystyle \sum _{n=2}^{N-1}\left(n-1\right){\pi }_{n,0}}+{\displaystyle \sum _{n=3}^{n_e\left(1\right)+1}\left(n-2\right){\pi }_{n,2}}+{\displaystyle \sum _{n=2}^{n_e\left(1\right)+1}\left(n-1\right)}\left({\pi }_{n,1}+{\pi }_{n,3}+{\pi }_{n,4}\right).$

(3) 顾客有效排队率为

${\lambda }_e=\lambda \left({\displaystyle \sum _{n=0}^{N-1}{\pi }_{n,0}}+{\displaystyle \sum _{n=0}^{n_e\left(1\right)}{\displaystyle \sum _{j=1}^4{\pi }_{n,j}}}\right).$

(4) 顾客平均逗留时间为

$E\left(W\right)=\frac{E\left(L\right)}{{\lambda }_e}.$

(5) 顾客平均等待时间为

$E\left(W_q\right)=\frac{E\left(L_q\right)}{{\lambda }_e}.$

(6) 令$C_{m_p}、C_{m_c}、C_{b_p}$ 分别为单位时间内维修主服务台不完全故障、完全故障和备用服务台的成本，$C_{u_m}、C_{u_b}$ 为主服务台和备用服务台单位时间使用成本，系统单位时间平均社会收益为

$\Psi ={\lambda }_{e}R-C_{w}E\left(L\right)-C_{u_m}{\displaystyle \sum _{i=1}^{n_e\left(1\right)+1}{\displaystyle \sum _{j=1}^2{\pi }_{i,j}}}-C_{u_b}\left({\displaystyle \sum _{i=1}^{n_e\left(1\right)+1}{\displaystyle \sum _{j=2}^4{\pi }_{i,j}}}-{\pi }_{1,2}\right)-C_{m_p}{\pi }_{0,2}-C_{m_c}\left({\pi }_{0,3}+{\pi }_{0,4}\right)-C_{b_p}{\pi }_{0,5}.$

4.2. 数值分析

本节通过数值实验分析各参数对进队阈值，平均逗留时间和社会收益的影响。设置基准参数为$R=10$ 、$C_w=1$ 、$N=4$ 、${\mu }_0=0.5$ 、${\mu }_1=3$ 、${\mu }_2=1.2$ 、${\mu }_3=0.8$ 、${\mu }_4=2.5$ 、${\mu }_5=1.2$ 、$C_{u_m}=4$ 、$C_{u_b}=3$ 、$C_{m_p}=20$ 、$C_{m_c}=25$ 、$C_{b_p}=20$ 、${\epsilon }_1=0.2$ 、${\epsilon }_2=0.05$ 、${\epsilon }_3=0.2$ 、${\xi }_1=0.2$ 、${\xi }_2=0.1$ 、${\xi }_3=0.2$ 。在后续分析中，未提及变化的参数均视为定值。

4.2.1. 进队阈值

首先关注两类服务台工作速率和不完全故障率对进队阈值的影响，分别设置${\mu }_1\in \left\{3,4,5\right\}$ 、${\epsilon }_1\in \left[0.05,0.5\right]$ ；${\mu }_4\in \left\{2.1,2.5,2.9\right\}$ 、${\epsilon }_3\in \left[0.05,0.5\right]$ ，得到进队阈值的变化趋势图2和图3。如图2所示，随着${\epsilon }_1$ 的逐渐增大，$n_e\left(1\right)$ 逐渐减小。这与实际相符，${\epsilon }_1$ 越大，系统处于主服务台高速率正常服务的概率越小，从而$n_e\left(1\right)$ 逐渐减小。同时对于相同的${\epsilon }_1$ ，主服务台正常服务速率${\mu }_1$ 越大，$n_e\left(1\right)$ 越大。从图3可知，${\mu }_4$ 和${\epsilon }_3$ 对$n_e\left(1\right)$ 的影响非常小，因为主服务台发生完全故障的概率很小，所以备用服务台单独服务的机会较少，对进队阈值的影响小。

Figure 2. Variation of $n_e\left(1\right)$ with ${\mu }_1$ and ${\epsilon }_1$

图2. $n_e\left(1\right)$ 随${\mu }_1$ 和${\epsilon }_1$ 的变化

Figure 3. Variation of $n_e\left(1\right)$ with ${\mu }_4$ and ${\epsilon }_3$

图3. $n_e\left(1\right)$ 随${\mu }_4$ 和${\epsilon }_3$ 的变化

4.2.2. 平均逗留时间

平均逗留时间是评价系统的一个重要指标，参数设置与上两幅图的完全相同。从下图4图5可以看出随着故障率${\epsilon }_2$ 和${\epsilon }_3$ 的增大，$E\left(W\right)$ 不断上升；同时对于同一故障率，服务台的速率越高，$E\left(W\right)$ 越小。该现象符合实际情况，故障率越小，系统处于高服务率的时间越长，服务率越高，系统队长越短，逗留时间越小。

Figure 4. Variation of $E\left(W\right)$ with ${\mu }_1$ and ${\epsilon }_2$

图4. $E\left(W\right)$ 随${\mu }_1$ 和${\epsilon }_2$ 的变化

Figure 5. Variation of $E\left(W\right)$ with ${\mu }_4$ and ${\epsilon }_3$

图5. $E\left(W\right)$ 随${\mu }_4$ 和${\epsilon }_3$ 的变化

4.2.3. 社会收益

最后讨论各参数对收益$\Psi$ 的影响。首先考虑$N$ 值的影响，设置$N\in \left[1,10\right]$ 、${\mu }_1\in \left\{3,4,5\right\}$ ，分别设置$\left({\epsilon }_1,{\epsilon }_2,{\epsilon }_3\right)=\left(0.2,0.05,0.2\right)$ ；$\left({\epsilon }_1,{\epsilon }_2,{\epsilon }_3\right)=\left(0.4,0.1,0.4\right)$ ，得到图6和图7。从下图可知，随着$N$ 的增大，收益$\Psi$ 都是先上升后下降，并且在两倍故障率下有右偏的倾向。这是因为当$N$ 值较小时，启动频繁成本高；当$N$ 较大时，系统激活困难，处于工作休假期的顾客增多，$\Psi$ 下降。当系统容易发生故障时，适当增大$N$ 可以增大$\Psi$ ，因此制定适当的$N$ 尤为重要。

Figure 6. Variation of $\Psi$ with ${\mu }_1$ and $N$

图6. $\Psi$ 随${\mu }_1$ 和$N$ 的变化

Figure 7. Variation of $\Psi$ under double failure rate

图7. 两倍故障率下$\Psi$ 的变化

分别设置${\mu }_1\in \left[3,4\right]$ ，${\epsilon }_1\in \left[0.05,0.5\right]$ ；${\mu }_3\in \left[0.5,1.1\right]$ ，${\epsilon }_1\in \left[0.05,0.5\right]$ ；${\mu }_4\in \left[2.1,2.9\right]$ ，${\epsilon }_3\in \left[0.05,0.5\right]$ 。得到社会平均收益$\Psi$ 随各服务率和故障率的变化趋势图8。随着故障率的上升，$\Psi$ 逐渐减小，因为维修成本所占比重高。同时随着${\epsilon }_1$ 的上升，$\Psi$ 随${\mu }_1$ 的增长趋势逐渐放缓，此时主服务台处于${\mu }_1$ 服务速率的概率降低，${\mu }_1$ 的边际收益下降；${\mu }_3$ 则相反，当${\epsilon }_1$ 上升时，系统更容易处于备用服务台辅助工作的状态上，因此${\mu }_3$ 的边际收益递增。$\Psi$ 随${\mu }_4$ 变化平稳，对${\epsilon }_3$ 敏感。${\epsilon }_3$ 较小时，$\Psi$ 随其近似二次下降；${\epsilon }_3$ 进一步增大后下降趋缓，原因是：${\epsilon }_3$ 上升会提高备用服务台不完全故障频率，增加系统低服务能力时长与拥塞、修理成本；后续故障状态概率的占比已较高，边际影响递减。

(a) $\Psi ({\mu }_1,{\epsilon }_1)$ (b) $\Psi ({\mu }_3,{\epsilon }_1)$ (c) $\Psi ({\mu }_4,{\epsilon }_3)$

Figure 8. Variation of $\Psi$ with service rates and failure rates

图8. $\Psi$ 随各服务率和故障率的变化

5. 仿真模拟

为验证本文推导的稳态性能指标与社会收益解析表达式的正确性，本节采用离散事件仿真对系统进行数值验证。到达率$\lambda$ 在区间[0.6, 1.0]内以步长0.02均匀选取21个观测点，其余参数与上一节数值实验的基本参数完全相同。仿真配置如下表1：

Table 1. Simulation experiment configuration

表1. 仿真实验配置

预热期	统计期	重复试验次数	置信水平
${10}^6$	$5\times {10}^6$	100	95%

首先运行${10}^6$ 个仿真时钟单位以消除初始状态影响，并在稳定运行$5\times {10}^6$ 个仿真时钟单位内统计平均等待队长$E\left(L_q\right)$ ，平均等待时间$E\left(W_q\right)$ 及单位时间平均社会收益$\Psi$ 。对每个$\lambda$ 取100次独立重复试验，计算仿真均值及其标准化残差z = (仿真值 − 解析值)/标准误的95%置信区间，得到如图9所示的两组对比图。

从图中可以看出：各性能指标的解析曲线与离散事件仿真均值高度一致，在标准化残差的21个观测点中只有$\Psi$ 的一个点位于95%置信区间外，符合比例为95.2%。高于统计检验的置信水平，说明本文稳态性能指标与社会收益解析表达式得到了数值验证。需要指出的是，部分指标在局部区间呈现非单调变化，主要是因为均衡阈值策略$n_e\left(1\right)$ 随到达率$\lambda$ 的调整引起系统在不同服务模式之间的占比变化，从而改变有效服务能力与等待结构；$E\left(W_q\right)$ 刚开始出现缓慢下降，以及高负荷区域等待成本快速上升导致$\Psi$ 出现回落，均与上述机制一致。

(a) $E\left(L_q\right)$ (b) $E\left(W_q\right)$ (c) $\Psi$

(d) $z_{E(L_q)}$ (c) $z_{E(W_q)}$ (f) $z_{\Psi }$

Figure 9. Simulation and residual analysis

图9. 仿真模拟与残差分析

6. 结论

本文研究了带N策略、工作休假、主服务台两类故障及备用服务台的M/M/1 + 1排队系统。通过建立连续时间马尔可夫链模型，推导了稳态概率分布与顾客进队阈值策略，并计算了平均队长、逗留时间及社会收益等性能指标。数值实验表明：主服务台故障率与维修策略对系统性能影响显著，而合理设置N策略可优化社会收益，研究结果为具有时段性客流与故障异质性的服务系统提供了运行调控与资源配置的理论依据。

此外，本研究还将向以下两个方向拓展：一是在系统信息不可视的情况下，顾客仅能依据历史经验或平均状态进行决策，此时系统模型可扩展为拟生灭过程；二是可引入顾客优先级机制，如区分普通顾客与优先顾客，研究其在故障，休假等状态下的均衡行为，以进一步贴近实际运营中的灵活调度与收益管理需求。