关于易腐产品的M/M/1 生产服务库存模型的最优控制

王 赛, 岳德权, 张媛媛, 田瑞玲

(燕山大学理学院，秦皇岛 066004)

1 引言

随着经济的迅速发展，企业面临的市场竞争越来越激烈，生产制造企业需要尽可能的降低库存管理成本才能有可观的经济效益.对于生产制造企业来说，库存不仅用来减少顾客需求，也可以用来调节生产，使得生产能够平稳的进行.因此生产库存模型及其库存的控制管理一直是库存理论的重要研究内容.

生产库存策略是生产库存系统的重要研究问题.施文武等[1]研究了一种多周期随机需求的生产库存控制系统，给出了费用函数，并设计算法求得最优控制策略.Benjaafar 等[2]研究了带有不耐烦顾客的生产库存系统，确定了最佳生产策略，分析了最优库存控制策略对于各项性能指标的影响.Krishnamoorthy 和Viswanath[3]研究了具有(s,S)策略的带有随机服务时间的生产服务库存系统，利用拟生灭过程方法得到了系统的稳态分布和费用函数.Baek 和Moon[4]研究了具有(s,S)策略的多个服务员的生产服务库存系统，给出了系统的性能指标和费用函数，并通过数值例子确定了系统最优控制策略.Axs¨ater[5]研究了具有(s,S)策略且容量有限的生产库存系统，假设系统生产时间服从伽马分布，利用M/G/1 的相关理论，得到了最优的生产策略.Nair 和Jose[6]研究了具有(s,S)策略，两种服务模式和重试需求的生产系统，求出了系统的性能指标，并建立了费用函数.

在上述关于生产库存系统的研究中，都没有考虑产品的易腐性.我国每年因为易腐性产品造成的损失之和高达千亿[7]，所以易腐库存系统的分析是很重要的.易腐性产品是指那些必须在有限时间内售出，否则将发生变质、损坏、挥发、过期且必须进行清仓处理的商品，如生鲜食品、水果、蔬菜、海鲜等.其显著特点是在储存和流通的过程中产品的数量会因腐烂变质、挥发、失效等逐渐减少.

易腐品的库存问题已经引起了学者的广泛关注.Kouki 等[8]研究了具有多种类型产品的易腐库存问题，并且给出了系统的最优库存控制策略.Sivakumar[9]研究了具有(s,S)策略，带有有限个顾客的易腐品库存系统，根据稳态方程求出了稳态概率向量，并建立了费用函数.Vijaya Laxmi 和Soujanya[10]分析了带有服务中断、重试需求和负顾客的易腐品库存模型，并利用矩阵分析法求出了模型在稳态下的联合概率分布，讨论了系统的各项性能指标，确定了费用函数并得到了库存的最优控制方案.Shophia Lawrence 等[11]研究了顾客有限的易腐产品的库存系统.假设服务时间服从位相分布，产品寿命服从指数分布，利用拉普拉斯变换方法推导出了系统的性能指标，构建了费用函数.Manuel 等[12]研究了具有(s,S)策略的两类客户的易腐产品库存系统，求得了稳态概率向量，并给出了费用函数的表达式.

在上述关于易腐产品的库存控制文献中，都是关注库存系统的性能分析和库存策略问题.目前在生产服务库存系统中，关于易腐产品的生产库存策略的研究工作较少.Sangeetha 等[13]研究了带有重试需求的易腐品生产库存系统.假设系统的顾客数量有限，基于(s,S)策略建立半马尔可夫过程，利用线性规划方法求出系统的最优生产速率，并给出数值例子进行分析.本文讨论了易腐产品的M/M/1 生产服务库存模型.假设顾客数量是无限的，规定当系统库存为零时，允许顾客进入系统.根据拟生灭过程理论，我们求得了系统的平衡条件，进一步得到了稳态概率向量矩阵几何解，性能指标及成本函数.最后给出数值例子，利用遗传算法实现最优解的有效搜索，得到了系统的最优控制策略.

2 模型描述

我们考虑具有(s,S)生产策略的易腐产品的马尔可夫生产服务库存模型，可描述为：

顾客需求的到达服从参数为λ 的泊松过程，顾客到达系统后形成一个队列，每位顾客的库存需求量为一个单位库存.系统中只有一个服务员，服务过程需要一定的时间，系统服务时间服从参数为µ的指数分布，服务规则为FCFS.

系统采用(s,S)生产策略，即当系统的库存水平下降到s 时，生产系统立即启动，并且每次只生产一个，直至库存水平达到最大值S 时，才停止生产，其中s < S.生产时间服从参数为η 的指数分布.

产品寿命服从参数为γ 的指数分布，变质后的产品将不再具有价值.系统是无损失的，即当系统库存为零时，顾客可以进入系统等待.假设需求到达过程，服务过程和生产过程是相互独立的.

3 系统平衡条件

3.1 状态过程

定义系统的状态过程为ψ = {(N(t),I(t),C(t)); t ≥0}，其中N(t)表示时刻t 系统中的顾客数，I(t)表示时刻t 系统的库存水平，C(t)表示时刻t 系统的生产状态，其中C(t) = 0 表示生产系统关闭，C(t) = 1 表示生产系统开启.根据(s,S)策略，C(t)应满足如下关系

过程ψ ={(N(t),I(t),C(t)); t ≥0}的状态空间为

将状态按字典序排列可得过程ψ ={(N(t),I(t),C(t)); t ≥0}的无穷小生成元如下

其中

其中

D1是(2S −s)×S 维矩阵，D2是(2S −s)×(S −s)维矩阵，C 是(2S −s)×(2S −s)维矩阵，I 是(2S −s)维单位阵，I1是(2S −s −1)维单位阵.

3.2 系统平衡条件

由矩阵Q 的结构可知过程ψ = {(N(t),I(t),C(t)); t ≥0}是拟生灭过程[14].令H =C+B+A，则有

H =(H1H2),

其中

L1(i)=−(η+µ+iγ), 1 ≤i ≤S −1,

L2(i)=−(µ+iγ), s+1 ≤i ≤S, G(i)=iγ+µ, 1 ≤i ≤S,

H1是(2S −s)×S 维矩阵，H2是(2S −s)×(S −s)维矩阵.

定义H 的平稳概率向量为

π =(π(0,1),π(1,1),··· ,π(S −1,1),π(s+1,0),π(s+2,0),··· ,π(S,0)),

平稳概率向量π 满足如下条件

其中e1是所有元素为1 的(2S −s)维列向量.通过计算整理，概率向量π 的各个分量可表示为

其中

由文献[14]知，过程ψ ={(N(t),I(t),C(t)); t ≥0}正常返的充分必要条件为πAe1<πCe1.经计算整理可得

通过计算，πAe1<πCe1等价于

所以，公式(1)就是系统稳态平衡的充分必要条件.

注1 根据给出的πAe1和πCe1的表达式，可以看出πAe1表示系统的到达率，πCe1表示当库存非零时系统的服务率，即系统的有效服务率.根据排队论的知识，当系统的到达率λ 小于系统的有效服务率时，系统达到平衡状态.

4 矩阵几何解

本节利用拟生灭过程理论，首先给出系统的稳态概率分布的矩阵几何解，然后讨论率阵R 和边界状态概率的计算问题.定义稳态概率为

在稳态条件ρ < 1 下，Q 的平稳概率向量P 存在.相应于Q 的分块结构，稳态概率向量P 分块如下

P =(P0,P1,··· ,Pi,···),

其中

Pi={P(i,0,1),P(i,1,1),··· ,P(i,S −1,1),P(i,s+1,0),··· ,P(i,S,0)}, i ≥0.

稳态向量P 满足如下平衡方程

其中e 是元素都为1 的适当维数的列向量.根据公式(2)，我们可以得到如下方程组

由文献[15]知，系统的稳态概率向量有矩阵几何解

其中R 是矩阵二次方程

的最小非负解，其谱半径sp(R)<1.P0是方程组的唯一正解，其中0 矩阵是元素都是0 的(2S −s)维行向量，I 是(2S −s)维单位矩阵.

为了计算系统的平稳向量，我们需要求出方程(7)的最小非负解R.在计算R 时，我们采用改进的循环简约算法.循环简约算法在文献[16,17]中有详细的介绍，循环简约算法步骤如下：

步骤1 令

步骤2 矩阵^B(j)和率阵R 由以下循环迭代方法求得

直至迭代得到R 值的差的范数满足‖R(j+1)−R(j)‖≤ε，从而求得率阵R；

步骤3 由方程(8)及Pi=P0Ri, i ≥0，即可求得平稳分布.

注2 为计算系统的平稳概率向量，我们需要求解矩阵二次方程(7)的最小非负解.在本文中，矩阵B、矩阵C 和矩阵A 均不是特殊矩阵，若直接利用方程(7)求解率阵R 则计算过程过于复杂，不易直接求出.例如，当采用(1,2)策略时，为求解R 需要求解一个九元二次方程组，求解过程相当复杂.因此我们使用上述算法，可近似求出率阵R 的数值计算结果.

例1 给定系统的参数值s = 1, S = 2, λ = 3.5, η = 5, γ = 0.05, µ = 8，根据率阵R 的计算方法，可近似求出R 如下

利用方程组(8)求解向量P0，结果如下

P0=(0.04132849666420,0.07048216383424,0.09789189421422),

再由公式(6)可以计算出系统的稳态概率向量.

5 稳态性能指标及成本分析

5.1 稳态指标

根据上述稳态概率向量的表达式，易得系统稳态性能指标，具体如下：

1) 平均队长

其中δ1是元素为1 的(2S −s)维列向量.

2) 平均等待队长

3) 平均库存

其中δ2为(2S −s)维列向量，且δ2=(0,1,2,··· ,S −1,s+1,s+2,··· ,S)T.

4) 平均易腐率

5) 平均生产率

其中δ3=(a b)T，a 是元素都为1 的S 维行向量，b 是元素都为0 的(S −s)维行向量.

6) 平均生产启动率

其中δ4为第(S+1)个元素为1，其余元素为0 的(2S −s)维列向量.

5.2 成本函数分析

假设该库存系统的成本主要由库存保管成本、生产成本、设备启动所需要的固定成本、产品腐烂的损失成本以及系统中每个顾客的等待成本组成.令单位时间单位库存的保管成本为Cinv，单位时间的生产成本为Crp，每单位产品腐烂的成本损失为Cp，系统中每个顾客的等待成本Clq，每次重启设备所需要的固定成本为K，则系统的成本函数C(s,S)为

显然，系统的成本函数是关于库存控制决策变量s, S 的非线性函数，且决策变量是离散的整数型变量.因为成本函数的表达式具有高复杂性，所以对于最小值的计算有一定难度.对于具体的数值例子，本文采用遗传算法求解最优值.

5.3 遗传算法

遗传算法(GA)是模拟自然界遗传、变异、适者生存的进化思想来求解优化问题，对目标函数和约束条件的要求较低.对于存在有限的可行解空间的纯整数规划问题，可采用适当的遗传算法得到最优解，本文采用了丰建荣等[18]提出的利用二进制编码的遗传算法，具体步骤如下：

步骤1 随机产生初始种群，个体数目一定.如果个体无效，那么重新随机生成初始个体，直至有效；

步骤2 评价群体的适应度，找到最好的染色体；

步骤3 用轮盘赌策略对每一代种群中的染色体进行选择；

步骤4 进行个体交叉操作，判断个体的有效性.如果是无效个体，则随机生成一个交叉位置进行交叉，直至有效；

步骤5 进行个体变异操作.如果产生无效个体，则重新进行变异操作，直至有效个体产生；

步骤6 由选择、交叉和变异产生新一代种群；

步骤7 对新种群适应度评价，找到最好的染色体，将它与上一次进化中最好的染色体比较，记录每一代进化中最好的适应度和平均适应度；

步骤8 如果不满足算法终止的条件，转到步骤3.否则，输出当前最优个体，算法结束.

6 数值分析

本节通过数值例子来考察系统参数的敏感性.根据上一节设置的遗传算法，设置种群大小为50，交叉率为0.6，变异率为0.01，算法终止条件为100 次(通过反复实验，程序迭代100 次结束，稳定收敛到最优解)，并用Matlab R2014a 编程进行数值实验.

表1 给出了参数λ 取不同的值时，系统的最优控制策略和一些性能指标的取值变化，其中参数设置为η = 5, γ = 0.05, µ = 8, Cinv= 10, Crp= 20, Cp= 50, Clq=200, K =2500.

表1 参数λ 的敏感性分析

由表1 可见，随着参数λ 的增大，系统的安全库存水平和最大库存水平都呈上升趋势.系统的平均库存、平均易腐率、平均生产速率、平均等待队长逐渐增加，平均生产启动率逐渐减小，最终成本函数逐渐增加.对于企业来说，在系统成本参数确定的条件下，市场需求增加时，库存控制的最优策略应相应的上调安全库存水平和最大库存水平.

表2 给出了参数η 取不同的值时，系统的最优控制策略和一些性能指标的取值变化，其中参数为λ=3.5, γ =0.05, µ=8, Cinv=10, Crp=20, Cp=50, Clq=200, K =2500.

表2 参数η 的敏感性分析

由表2 可见，随着参数η 的增加，系统的安全库存水平，最大库存水平均呈下降趋势.系统的平均库存、平均生产率、平均易腐率先减后增，平均等待队长逐渐减小，平均生产启动率逐渐增大，最终导致系统的成本逐渐增大.对于企业来说，可以通过技术改进，更新设备等方法提高生产率.同时当生产效率提高时，库存控制的最优策略需要下调安全库存水平和最大库存水平，企业在调整时还需要注意生产效率过高会造成产品堆积.

表3 给出了参数γ 取不同的值时，系统的最优控制策略和一些性能指标的取值变化，参数为λ=3.5, η =5, µ=8, Cinv=10, Crp=20, Cp=50, Clq=200, K =2500.

表3 参数γ 的敏感性分析

由表3 可见，随着γ 的增加，系统最大库存水平逐渐减少.系统的平均库存、平均生产启动率逐渐减少，平均生产率、产品的平均易腐率逐渐增加，平均等待队长先增后减，最终导致系统费用增加.易腐产品的寿命参数γ 越小，产品的腐烂速度则越慢，系统成本越小.对于企业来说，可优化易腐品的保存方法，延长产品寿命，进而降低库存成本.

表4 给出了参数µ取不同的值时，系统的最优控制策略和一些性能指标的取值变化，参数设置为λ=3.5, η =5, γ =0.05, Cinv=10, Crp=20, Cp=50, Clq=200, K =2500.

表4 参数µ的敏感性分析

由表4 可见，随着参数µ的增加，系统的安全库存和最大库存水平保持稳定.系统平均生产启动率逐渐增加，但增幅很小.平均生产率、平均库存和平均易腐率、平均等待队长逐渐减小，系统成本逐渐减小.在系统成本参数确定的条件下，服务速率越高，系统的成本则越低.所以对于企业而言，应该提高员工的服务效率，进而降低系统的库存成本.

7 结论

本文研究了有关易腐产品的M/M/1 生产服务库存模型，建立了关于系统中的顾客数、库存水平和生产状态三个随机变量的三维拟生灭过程.利用拟生灭过程理论，求出了系统的稳态平衡条件，给出了系统的稳态概率的矩阵几何解，并求出一些系统稳态下的性能指标.进一步构建了成本函数，针对系统成本函数的非线性、整数型变量的特点，运用遗传算法，实现了最优解的有效搜索.通过数值实验分析了系统参数的敏感性，这些分析结果，有助于生产型企业在外部参数发生变化时，根据实际情况及时调节库存控制策略，对库存管理实践具有指导借鉴意义.