《1 前言》

1 前言

在军事装备、民用高技术产品等领域，对于具 有较高可靠性，并且流通量较低的昂贵可修备件， 管理者通常选择多级（S - 1，S）库存策略。Sherbrooke[1] 最早根据这种库存策略开发了METRIC模 型，他选择费用约束下的期望缺货量最小为优化目 标，假设各级维修能力无限。Grave[2] 针对 Sherbrooke 研究中设计了一种便于理解的近似计算方 法。文献[3，4]则针对METRIC模型无限修理能力 的假设，研究了维修台有限情况下的多级可修备件 库存模型。付兴方等[5] 研究了两级库存系统备件分 配优化方法。Lau和Song[6] 采用非平稳Poisson过程 研究了维修能力有限的两级可修库存系统。这些 研究均采用Poisson过程描述部件的失效过程，假设 维修时间等服从指数分布，这样的假设具有一定的 不合理性[7] 。Kim等[8] 为了放宽假设条件，进一步假 设维修时间服从一般分布，维修台有限；但该研究 仍是建立在需求到达是Poisson流的基础上，因此假 设条件仍略显严格。

本文针对已有研究假设条件过于严格的情况，提出采用马尔可夫到达过程（MAP）描述备件的需 求情况，考虑有限维修设施的情况，假设故障件维 修时间、备件运输时间以及采购时间均服从一般分 布，采用phase-type（PH）分布进行表示；建立了一个 描述能力更强的两级可修备件（S -1，S）库存模型。 采用MAP描述需求到达过程时，（S -1，S）策略通常 并不是最优策略，这是因为当需求过程比较平稳 时，适当增加订货延迟时间是有助于降低成本的。 但在军事和工业领域，备件缺货除了对费用的影 响，往往还会导致任务失败或增加设备的安全风 险，因此采用（S -1，S）策略更加稳妥，也更符合实际 情况。

利用MAP和PH分布开展研究具有如下好处。 1）将备件需求过程表示为MAP形式便可以利 用“MAP类在[0，+∞）的概率空间上，对平稳点过程 类似稠密[9] ”这一性质。该性质使得MAP类涵盖了 许多常用的点过程[10] ，如泊松过程、PH 更新过程 等；并且实际的备件需求往往具有突发性，而MAP 本身结构的特点[11] 使它能够很好的描述和分析突 发现象。因此，MAP非常适合描述备件的需求发生 规律。

2）采用 PH 分布表示一般分布是因为：任何分 布总可以选择一个适当的PH分布把它拟合到任意 精确的程度[12] ，因此PH分布可以作为许多假设分布 的一般表述，采用 PH 分布描述维修时间等能有效 避免指数分布等经典分布的局限性。

《2 预备知识》

2 预备知识

首先介绍PH分布和MAP的有关概念。

定义 1[7] [0,+∞) 上的概率分布函数 F (∙) 称为 PH分布，当且仅当它是一个有限状态马尔科夫过程 的吸收时间分布，有分布函数

式（1）中，T 为 m 阶方阵； α=(α1 ，α2 ，...，αm ) 为其瞬态的初始概率向量； e 为元素均为1的 m 阶列向量； (α，T ) 称为该PH分布的 m 阶表示。PH分布中的每一个瞬态称为相位。

定义2 [11] 一个有限不可约马尔可夫链具有状态空间 S ={1，2，...，m } ，设 D 为该马尔可夫链的无穷小生成元。状态i的逗留时间服从参数为 λi 的指数分布，在该状态即将结束时，有下列两个事件之一发生：

1）存在一个马尔可夫链，使得从状态 i 到 j（ ）的转移中，有一个事件到达的概率为 ；

2）存在一个马尔可夫链，使得从状态 i 到 j（； ）的转移中，没有事件到达的概率为 。

可以看到只有通过一次事件到达，才能使得状态 i 返回到自身。并且有：

定义矩阵和，其 中 ：，;,且;,, 。

和 分别是控制没有事件到达和有单个事件到达的状态转换矩阵，它们均为 m 阶子随机矩阵，可以给出无穷小生成元 。

然后，定义 为在 (0，t ] 内事件到达的次数；J(t ) 为在时刻 t 嵌入马尔可夫链所处的状态，其状态空间为 {:1≤≤m } 。则 {，} 为一个两维马尔可夫过程，其状态空间为 {(0，1≤≤m } 。而 {，} 被称为马尔可夫到达过程。

MAP是半马氏过程的特例，其状态转移概率矩阵为：

设 是无穷小生成元为 的稳态概率向量，则有：   。并且 表示在 MAP 进入稳态后，单位时间内事件到达的期望个数。

下面给出Kronecker乘积的定义，它在后面的模型中被大量使用。

定义3 [13] 如果 A 和 B 分别为 和 的矩阵 ，它们的 Kronecker 乘 积 则为 的矩阵，且：

《3 问题描述》

3 问题描述

本文假设如下。

1）一个中心仓库通常负责 N 个基地的备件供应和维修，在中心仓库和基地分别存有 和 （=1，...，N ）个某型备件，并且在中心仓库和每个基地均有一个维修站。

2）基地 （ =1，...，N ）的故障到达是一个MAP，具有 m 阶表示。

3）当基地 i 中出现一个备件需求时，若该基地中仍有库存，则立刻满足该需求。基地 备件需求处理过程示意图见图1。该更换下来的故障件有 r 的概率可在基地进行维修，否则需交给中心仓库处理。

4）故障件交还中心仓库得同时，基地向中心仓库申领一个备件，若中心仓库有库存则立即下发，否则需要等待维修完成或订货到达；

5）在中心仓库，故障件有 p 的概率可修， 1-p 的概率报废；当发生报废时，中心仓库需向生产厂家订购；

6）基地 （ =1，...，N ）的维修站维修时间服从PH分布，具有 阶表示 ；中心仓库维修站维修时间服从PH分布，具有 阶表示 ；所有维修站均采用先到先服务原则；

7）采购时间服从 PH 分布，具有 阶表示。

《图1》

图1 基地 i 备件需求处理过程示意图

Fig. 1 The disposal process for spare parts demand of base i

8）从中心仓库向各基地的运输时间服从PH分布，具有 阶表示 ；而从基地向中心仓库运送故障件的时间则可以计入中心仓库的维修时间或采购时间中，不单独列出。

9）所有节点均选择（S -1，S）补充策略。

《4 系统特性》

4 系统特性

下面首先对整个系统的备件供应流程进行分析。

《4.1 备件供应流程》

4.1 备件供应流程

根据假设条件，可以将整个备件供应流程划分为4个部分，如图2所示。

《图2》

图2 备件供应流程划分示意图

Fig. 2 The partition of spare parts supply flow

令 、 和 分别表示基地 i 的工作部件数、维修队长和缺货量；  和  表示在中心仓库的等待维修队长和采购队列队长，则 表示中心仓库不可用件数量；而中心仓库的缺货量为 BO0 ，其中基地 在中心仓库未满足订单的比例为 ； OH0 和 分别表示中心仓库和各基地的可用件数量；在运往基地 途中的备件数为 。下面研究这些参数之间的关系：

由上述关系式可知， BO0 和 可以表示为多个随机变量的卷积形式。而又由假设条件可知、 、  和 均是MAP/PH/1排队系统的队长，则通过研究MAP/PH/1排队系统可知它们所服从的分布。

《4.2 MAP/PH/1排队系统》

4.2 MAP/PH/1排队系统

令=( , , )，其中、、分别为在时刻 t 的系统顾客数，顾客到达过程相位和服务台工作相位；另假设服务时间服从 k 阶 PH分布。则为连续时间马尔可夫链，其状态 空 间 Ω 可 以 表 示 为 ： Ω=⋃ 。 其 中 ，={(0，)：1≤m }表 示 系 统 内 没 有 顾 客 ； ={(s，，j )：1s，1m，1j k } 表示系统内有 s 个顾客，其中一个顾客正在接受服务。

通常称 为状态空间的宏状态。对宏状态按字典序排列，可以给出该马尔可夫链的无穷小生成元矩阵：

式（9）中， ； ； ； ； ；。

对维修台工作情况进行分析，发现维修台的状态转移过程亦为一个MAP，其中无穷小生成元矩阵为 ；因此，存在一组稳态概率向量 使得：。

定义矩阵，则 为系统队列长的状态转移率矩阵。

定理1 矩阵 和 均是不可约随机阵，并且矩阵 的稳态概率向量。

证明 由 和 的定义已知它们均为不可约随机阵，而：

因此矩阵 的稳态概率向量。定理得证。

由MAP和PH分布性质知，顾客到达率和服务率分别为：。

Neuts [14，15] 给出了求矩阵 Q 稳态概率向量的一个重要矩阵 —— 率阵 R ， R 是矩阵方程的最小非负解。

定理2 不可约马氏链 { } 是正常返的，当且仅当：

证明 根据文献[7]定理2.4知若该马氏链正常返 ，则 必 须{ } 的 R 的谱半径 sp(R )<1 ，即，该式正是经典排队系统条件向矩阵空间的推广。定理得证。

在定理1、定理2的基础上，运用矩阵几何解理论，不难得出：

推论 1 若 ρ <1 ，则矩阵 Q 的稳态概率向量θ =(θ0 ，θ1 ，θ2 ，...) 存在，且，k >1；而和满足下面的矩阵方程：

并且由下式归一化后可得 θ ：

θ 就是系统在任意时刻的队长分布函数，因此容易给出任意时刻队长的各阶矩参数。根据推论1可以获得各维修队列、采购队列的平均队长及方差等参数。

此外，备件从中心仓库向基地 i 运送的过程也可以看作是一个MAP/PH/1排队系统，这是因为备件进入运输渠道是一个MAP，而运输时间可以视为服务台服务时间，因此根据推论1容易确定运输中的备件个数分布。

《4.3 缺货量分布》

4.3 缺货量分布

记 为中心仓库不可用备件数，则中心仓库缺货量的期望值可以表示为：

而 BO0 的两阶矩为：

BO0 的方差为：

在到达过程为MAP的情况下，由MAP性质 [15]可知：

令 为基地 的不可用备件数，则基地 的期望缺货量则为：

而 BOi 的两阶矩为：

则BOi 的方差为：

在获得各基地缺货量的期望和方差后，Grave [2]通过大量数据分析认为缺货量分布可以近似表示为参数 和 的负二项分布，因此本文亦直接采用该结论，由式（20）、式（21）可确定出 、 的值。

《5 库存优化模型》

5 库存优化模型

定义 PBO( ， ) 为在给定 和 值后，基地 发生缺货的概率，由4.3节获得的缺货量分布容易得到：

因为各基地是一个并联系统，因此该多级库存系统的备件满足率可以表示为：

以总费用 c total 最小为目标函数，确定合适的S *=(S 0 ，S 1 ，...，S N ) 使得整个系统满足备件满足率约束。令 c 1 、 c 2 、 c 3 、 c 4 和 c 5 分别为备件单价，单位备件库存费用，单位备件基地维修费用，单位备件中心仓库维修费用和单位缺货费用，因此该优化模型可以表示为如下形式：

式（23）中，和 分别表示各库存点库存量和维修队列长。

显然，备件满足率与单个库存点备件数之间存在 单 调 递 增 关 系 ；但 A (S 0 ，S 1 ，...，S N ) 与{(S 0 ，S 1 ，...，S N )} 的关系由于涉及多个相互独立的变量，就非常复杂了。对这类非线性整数规划模型而言，许多成熟的智能算法均能有效获得最优解。本文对智能算法不做深入研究，而直接采用遗传算法进行求解。

下面通过示例演示MAP和PH分布能够有效描述两级可修备件的（S -1，S）库存模型。

《6 算例》

6 算例

某中心仓库负责4个基地的备件供应和维修保障任务。基地 i 的某型部件备件需求到达过程可以表示为MAP（）；而各基地的维修时间分布、向中心仓库运输故障部件的时间分布可以分别表示为PH（）和PH（），表1、表2是各基地的需求、维修及运输数据。

 

中心仓库的维修时间分布为PH（），采购备件的时间分布为PH，其中：

而在中心仓库可修的概率则为 p =0.79 。c 1 、 c 2 、 c 3 、 c 4 和 c 5 分别为备件单价，单位备件年库存费用，单位备件基地维修费用，单位备件中心仓库维修费用和单位缺货费用， c 1 =26 800 元， c 2 =90 元，c 3 =440 元， c 4 =1 900 元， c 5 =7 000 元。

《表1》

表1 各基地备件需求到达过程

Table 1 The spare parts demand process of each base

《表2》

表2 各基地备件维修时间分布及运输时间分布

Table 2 The distributions of repair time and transit time for each base

下面利用4.3节获得的缺货量分布，研究中心仓库和各基地备件数对整个系统的备件满足率的影响。由于本优化模型包含5个自变量，因此为了便于更直观地了解备件数与备件满足率之间的关系，以中心仓库备件数 S 0 、基地1备件数 S 1 同时对系统备件满足率的影响为例进行说明。

如图3所示，令其他基地备件数均为0时，可以看到整个系统的备件满足率早期随着 S 0 、 S 1 增加而增加幅度较大；但当 S 0 8 并且S 12 时，备件满足率就保持在0.711 2附近，说明要使满足率继续增加，必须调整其他基地的备件量。当其他基地的备件数均设为1时，系统备件满足率有大幅提升，迅速达到了0.937 4附近。而当其他基地的备件数均为5时，满足率上升则没有早期那么明显了，稳定在0.999 9附近。从图3可以看到，当给定其他基地的备件数时，S 0 、S 1 对满足率的影响十分相似。同时也表明无论中心仓库还是基地的备件数增加，均会使得备件满足率增加，只是增加的幅度与系统内其他仓库的备件数量密切相关；如果只增加 S 0 ,⋯,S4中任意一个，整个系统的备件满足率增长则十分有限。而当S 0 ,...,S4 均较大时，系统的备件满足率非常接近1，但此时仅购置费就会大幅上升。

《图3》

图3 中心仓库和基地1备件数对备件满足率的影响

Fig. 3 The influence of the spare amount of central depot and base 1 on the system spare fill rate

在对备件满足率的分析过程中，人工对各库存点的库存量进行调整，发现当 (S 0 ，...，S4 )=(6，0，2，2，1)时，满足率为0.9092；而 (S 0 ，...，S4 )=(5，1，1，1，1)时，满足率达到了0.921 5。这说明为了提高系统的备件满足率，除了增加总的备件量外，还必须合理分配各单位的备件数。

《7 结语》

7 结语

本文分析了以往多级可修备件库存模型中存在的问题，用MAP描述各基地的故障到达流，将维修与采购时间、运输时间均假设服从一般分布，建立了一个适应性更好的多级可修备件（S -1，S）模型，并给出了该系统缺货量分布。该类模型涉及到的库存优化算法是库存模型研究的一个重要方向，结合该类库存问题的实际背景和特点，开发更合理高效的启发式算法对本文模型在实际工作中获得有效应用有着重要的意义。本文建立的模型在计算时涉及到矩阵的运算，而目前高性能计算机和矩阵解析方法的应用能对大型矩阵的运算提供良好的支持，所以本文建立的模型具有很好的实用
价值。