自动化学报  2018, Vol. 44 Issue (2): 270-279   PDF    
双源采购跳跃-扩散库存控制模型
娄山佐1, 田新诚1, 吴颖颖1     
山东大学控制科学与工程学院 济南 250061
摘要: 供应中断和退货会引发库存短缺和剧烈波动,所以,如何缓解它们的影响,成为当前企业管理者亟待解决的难题.在采用双源采购策略防御库存短缺和跳跃-扩散过程描述库存水平变化条件下,利用连续时间Markov链、水平穿越和鞅理论,分别确定了库存水平分布及循环的期望费用和时间函数,在此基础上,构建了系统长程平均费用率模型.最后,仿真结果表明,供应商的可靠性和中断类型,对最优控制策略和系统费用产生较大影响.另外,双源采购策略能够有效缓解供应中断对库存的影响,尤其是,当供应商的可靠性较低或中断类型均为频率低持续时间长时.
关键词: 库存控制     供应中断     双源采购     跳跃-扩散     水平穿越    
A Jump-diffusion Inventory Control Model with Dual-sourcing
LOU Shan-Zuo1, TIAN Xin-Cheng1, WU Ying-Ying1     
School of Control Science and Engineering, Shandong University, Jinan 250061
Manuscript received : July 7, 2016, accepted: November 28, 2016.
Foundation Item: Supported by National Natural Science Foundation of China (61703241), Major Science and Technology Innovation Project of Shandong Province (2017CXGC0601), and Key Research and Development Project of Shandong Province (2016ZDJS02B03)
Author brief: TIAN Xin-Cheng  Professor at the School of Control Science and Engineering, Shandong University. His research interest covers motion control, mechatronics, intelligent control and optimization;
WU Ying-Ying  Lecturer at the School of Control Science and Engineering, Shandong University. Her research interest covers scheduling optimization of automatic picking system, logistics system planning, design and simulation
Corresponding author. WU Ying-Ying   Lecturer at the School of Control Science and Engineering, Shandong University. Her research interest covers scheduling optimization of automatic picking system, logistics system planning, design and simulation
Recommended by Associate Editor ZHAO Qian-Chuan
Abstract: Supply disruptions and returns may lead to stock shortage and fluctuation. So it is a very difficult problem for today's managers to mitigate their influence. Under the condition that a dual-sourcing strategy is utilized to tackle the stock shortage and a jump-diffusion process is adopted to express the inventory level process, the stationary distribution of the inventory level, as well as the expected cycle cost and time functions are derived by applying continuous-time Markov chain, level-crossing and martingale theorems. Subsequently, the functions are employed to develop a long-run average cost rate model. Finally, numerical results show that the suppliers' reliability and the nature of the disruptions have a big impact on the optimal policy and cost. Moreover, the dual-sourcing strategy can mitigate the influence of supply disruptions on inventory effectively, in particular, when the suppliers' reliability are lower or their disruptions are rare and long in nature.
Key words: Inventory control     supply disruption     dual sourcing     jump-diffusion     level-crossing    

随着市场经济发展, 建立在供应稳定可靠和库存单调变化假设之上的传统库存管理方法, 遇到现实严峻的挑战.如受雾霾和交通堵塞等事件影响, 我国零售企业缺货率约为$10\%$, 每年损失830亿元[1].与此同时, 用户退货又引起库存变化失去单调性, 呈现随机波动状态.如$2012$年唯品会上零售企业平均退货率高达$20\%$[2].在此环境下, 如何有效控制库存, 成为近年来供应链管理研究的热点.

与本文相关的有供应中断防御问题和库存波动控制问题两方面成果, 下面分别给予综述.

针对供应中断防御问题, 有许多成果[3].这里, 仅给出与双源采购有关成果.自Gerchak等[4]率先利用双源采购策略应对随机供应问题以来, 不少学者基于经典报童模型, 将该策略应用到防御供应中断问题.在假设不可靠(即可能中断)和可靠两供应商补货情况下, 文献[5]分析了库存防御、双源采购和应急订货策略, 对缓解供应中断的有效性; 文献[6]又考虑使用快运方式策略; 文献[7]还结合供应能力的期权契约.另外, 在假设供应商均不可靠情况下, 文献[8]通过比较单源和双源采购系统性能, 指出双源采购策略更能有效防御中断造成的库存短缺.最近, 文献[9-10]在零售商可利用中断信息和学习功能环境下得到同样结论.以上成果均针对周期盘点系统.而对应的连续盘点系统, 虽然实际应用广泛, 但因包含多个随机因素, 目前, 仅检索到两篇文献.在假设供应和中断持续时间均服从指数分布条件下, 文献[11]构建了单源和双源采购库存模型.随后, 文献[12]又将该模型推广到供应时间为Erlang分布环境.然而, 上述所有成果均假设库存为单调非增变化.

近年来, 为应对退货等引发的库存波动, 人们常用扩散过程近似库存水平变化, 简称扩散库存.如在用布朗运动描述库存水平过程条件下, 文献[13]基于脉冲控制法, 给出线性库存和短缺费的最优订货策略; 文献[14]将其拓展到凸库存费环境; 文献[15]更是包含了凹订货费用.当库存水平起伏很大时, 有学者又用跳跃-扩散(即Lévy)过程描述库存水平变化, 给出最优的库存策略[16]及处理速率和缓冲器大小[17].但这些成果均没涉及供应中断对库存的影响.

很遗憾, 因建模困难, 与供应中断和库存波动变化均有关的成果很少.文献[18]研究了库存水平含扩散项不可靠生产系统的控制问题.针对随机提前期造成的供应中断, 文献[19]给出最优脉冲库存控制策略.另外, 文献[20]构建了单源采购扩散库存控制模型.随后, 文献[21]又分析了跳跃-扩散库存的应急控制问题.然而, 当短缺费和应急采购费较高时, 这些成果很难取得理想效果.

为解决供应中断和退货环境下连续盘点库存问题, 本文利用水平穿越理论和鞅方法, 构建了双源采购跳跃-扩散库存控制模型, 并通过仿真实验, 证实此环境下双源采购策略的有效性, 同时, 分析了关键因素对系统的影响.最后, 给出有关库存管理启示.

1 问题的描述

一零售商销售某种商品, 并承诺用户不满意可退货.据统计, 用户平均需求率为$\mu$.退货到达为复合Poisson过程$Y(t)=\sum_{n=0}^{N(t)}V_n$, 这里, $N(t)$为参数是$\lambda$的Poisson过程, 表示到时间$t$退货到达批次; $V_1, V_2, \cdots$为独立的且参数是$\upsilon$的指数分布随机变量序列, 表示每批退货数量.为使系统稳定运行, 满足$\mu{-}{\lambda}/{\upsilon}{>}0$.设退回产品无缺陷, 经检查和重新包装等处理后(处理时间可忽略不计), 直接进入库存, 与新品一样满足未来用户需求.

为防御供应中断, 该零售商实施双源采购策略.设供应商$i$ (用$S_i$表示)的供应状态(用ON表示)和中断状态(用OFF表示)持续时间, 分别服从独立的参数为$\gamma_i$$\theta_i$的指数分布, $i=1, 2$, 且两供应商间的状态也相互独立.那么, 系统供应商运行情况, 可用$4$个状态的连续时间Markov链$I=(I(t))_{t\geq0}$描述, 这里,

$ \begin{equation*} I(t)= \begin{cases} 0, & S_1 为 {\rm{ON}} 且 S_2 为 {\rm{ON}} \\ 1, &\ S_1 为 {\rm{ON}} 且 S_2 为 {\rm{OFF}}\\ 2, & S_1 为 {\rm{OFF}} 且 S_2 为 {\rm{ON}}\\ 3, & S_1 为 {\rm{OFF}} 且 S_2 为 {\rm{OFF}} \end{cases} \end{equation*} $

零售商采用扩展的$(s, S)$订货策略.在提前期为$0$假设下, 当库存降到订货水平$s$时, $S_1$和(或) $S_2$若没发生中断, 则补货量分别为$q_{1}$和(或) $q_{2}$; 若均发生中断, 则补货量据中断结束时的库存水平确定.

具体订货过程, 按订货时刻$t$系统所处状态, 可分下列4种情况: 1)若$I(t)=0$, 则$S_1$$S_2$将库存补充到水平$\bar{q}_{0}=q_{1}{+}q_{2}{+}s$; 2)若$I(t)=1$, 则$S_1$将库存补充到水平$\bar{q}_{1}=q_{1}{+}s$; 3)若$I(t)=2$, 则$S_2$将库存补充到水平$\bar{q}_{2}=q_{2}{+}s$; 4)若$I(t)=3$, 则库存从水平$s$继续运行, 期间如果库存水平为$0$, 则需求丢失, 造成缺货费用.因$S_1$$S_2$中断同时结束的概率为$0$, 故当$S_1(S_2)$中断先结束时, 若库存水平小于$s$, 则$S_1(S_2)$将库存补充到水平$\bar{q}_{1}(\bar{q}_{2})$; 若库存水平大于$s$, 则$S_1$$S_2$不再补货.当库存又降到水平$s$时, 零售商再次向$S_1$$S_2$提出订货, 如此重复循环运行.假设$S_i$的固定和可变补货费用分别为$K_i$$k_i$, $i=1, 2$; 单位产品退货费用(包括产品的成本, 以及检查、换包装和重新入库等费用)为$\delta$; 单位产品的缺货费用为$\pi$; 单位时间内单位产品的库存费用为$h$.确定$q_1$$q_2$$s$, 使系统单位时间的长程平均费用最小.

2 模型的构建

若库存从水平$x_{0}$开始运行, 因它为$0$时需求丢失, 故补货发生前库存水平变化, 可表为在$0$点反射的跳跃-扩散过程$W=(W(t))_{t\geq0}$, 这里, $W(t)=x_{0}+X(t){+}L(t)$, 其中, $X(t)=Y(t){-}\mu t$是一个跳跃-扩散过程; $L(t)=-{\rm min}_{{0{\leq}\varsigma{\leq} t}}[(x_{0}{+}X(\varsigma))\wedge0]$是局部时间过程, 表示到时间$t$的短缺量.

据跳跃-扩散过程强Markov性知, 二维过程$( I, W )$是Markov过程, 且$W$还是一个更新过程.更新点可选$S_1$$S_2$补货完成, 即过程$(I, W)$的状态点$(0, \bar{q}_{_0})$.循环周期等于相邻更新点间的时间.图 1给出一次循环库存水平典型的样本路径(设$q_2 < q_1$).

图 1 一次循环库存水平典型的样本路径 Figure 1 A typical sample path of the inventory level in one cycle

定义$C_{{i0}}$$T_{{i0}}$分别表示状态为$i$时库存降到水平$s$ (若$i{\neq}3$, 因提前期为$0$, 故库存水平实际为$\bar{q}_{i}$), 自此直到下次循环开始的期望费用和时间.则据更新报酬定理知, 系统长程平均费用率为

$ TC(q_{1}, q_{2}, s)=\frac{{{\rm{E}}\left[ {{\rm{一次循环的费用}}} \right]}}{{{\rm{E}}\left[ {{\rm{一次循环的时间}}} \right]}}\;\;\;\;\;\;\;\;\;=\frac{C_{{00}}}{T_{{00}}}. $

显然, 要求最优$q_1$$q_2$$s$, 需确定循环的期望费用$C_{{00}}$和时间$T_{{00}}$.下面, 先求与它们有关的函数.

3 有关函数的求解

为确定系统短缺、库存、退货和补货费用, 需求解它们的期望量.又因在状态$3$供应商中断结束时库存水平是随机的, 故计算期望补货量, 需确定此时库存水平的分布函数.另外, 还要确定系统状态间的转移函数.下面, 分三部分给出它们的具体形式.

3.1 确定有关期望量和时间函数

设在状态$3$过程$W$从水平$s$开始到一供应商中断结束为止的运行时间为$\tau$.因$S_1$$S_2$中断持续时间, 分别服从独立的参数为$\theta_1$$\theta_2$的指数分布, 故据随机过程理论知, $\tau$服从参数为$\theta=\theta_1{+}\theta_2$的指数分布, 且$S_1$ ($S_2$)中断先结束的概率为$\vartheta_1=\theta_1/\theta$ ($\vartheta_2=\theta_2/\theta$).因此, 期望时间${\rm E}[\tau]=1/\theta$.

另设$W$从任一水平$z (z{>}s)$出发首次到达水平$s$的时间为$\tau(z)$, 即$\tau(z)={\rm inf}\{t{>}0{:} W(t)=s|W(0)=z\}.$

下面, 先给出对应时间$\tau$的有关期望量.在此基础上, 给出对应时间$\tau(z)$的有关量.

引理1.$\Pi={\rm E}[L(\tau)]$$R={\rm E}\int_0^\tau {\rm d}Y(t)$$H={\rm E}\int_0^\tau W(t){\rm d}t$.则时间$\tau$内的期望短缺量、退货量和库存量分别为

$ \begin{equation} \Pi=\frac{{\rm e}^{-\alpha_1(\theta)s}}{\alpha_1(\theta)} \end{equation} $ (1)
$ \begin{equation} R=\frac{\lambda}{\upsilon\theta} \end{equation} $ (2)
$ \begin{equation} H=\frac{s+\Pi}{\theta}+\frac{\lambda-\mu\upsilon}{\upsilon\theta^2} \end{equation} $ (3)

这里, $\alpha_1(\theta)$为常数, 其值由式(5)给出.

证明.确定上述函数采用的工具为下列Kella-Whitt鞅$M(t)$, 详细了解请参阅文献[22].

$\begin{equation*} M(t)=(\varphi(\alpha)-\beta)\int_0^t{\rm e}^{-\alpha W(\omega)- \beta\omega}{\rm d}\omega+{\rm e}^{-\alpha W(0)}- \end{equation*} $
$\begin{equation} {\rm e}^{-\alpha W(t)-\beta t}-\alpha\int_0^t {\rm e}^{-\beta\omega}{\rm d}L(\omega) \end{equation} $ (4)

这里, 跳跃-扩散指数$\varphi(\alpha)=\mu\alpha-\lambda\alpha/(\upsilon+\alpha)$.

$\varphi(\alpha){-}\theta=0$, 可得两根如下:

$ \begin{equation} \alpha_{1}(\theta)=\frac{\lambda{+}\theta{-}\mu\upsilon{+}\sqrt{ (\lambda{+}\theta{-}\mu\upsilon)^2{+} 4\theta\mu\upsilon}}{2\mu} \end{equation} $ (5)
$ \begin{equation} \alpha_{2}(\theta)=\frac{\lambda{+}\theta{-}\mu\upsilon{-}\sqrt{ (\lambda{+}\theta{-}\mu\upsilon)^2{+} 4\theta\mu\upsilon}}{2\mu} \end{equation} $ (6)

令式(4)中$\beta$$0$, 并对停时$\tau$应用最优抽样定理.因${\rm E}[M(0)]={\rm E}[M(\tau)]$, 故有:

$ \begin{equation} \varphi(\alpha){\rm E}\hspace{-1.0mm}\int_0^\tau\hspace{-1.0mm}{\rm e}^{-\alpha W(t)}{\rm d}t={\rm E}[{\rm e}^{-\alpha W(\tau)}]- {\rm e}^{-\alpha s}+\alpha{\rm E}[L(\tau)] \end{equation} $ (7)

又因

$ \begin{align} {\rm E}[{\rm e}^{-\alpha W(\tau)}]=&\int_0^\infty{\rm E}[{\rm e}^{-\alpha W(t)}]\theta{\rm e}^{-\theta t}{\rm d}t=\notag\\ &\theta{\rm E}\int_0^\tau{\rm e}^{-\alpha W(t)}{\rm d}t \end{align} $ (8)

将式(8)代入式(7), 化简得:

$ \begin{equation} (\varphi(\alpha)-\theta){\rm E}\int_0^\tau{\rm e}^{-\alpha W(t)}{\rm d}t=-{\rm e}^{-\alpha s}+\alpha\Pi \end{equation} $ (9)

首先, 令式(9)中的$\alpha$$\alpha_1(\theta)$, 左边等于$0$, 故期望短缺量

$ \begin{equation} \Pi=\frac{{\rm e}^{-\alpha_1(\theta)s}}{\alpha_1(\theta)} \end{equation} $ (10)

其次, 据$Y(t)$为复合Poisson过程得:

$ \begin{equation}{\rm E}[{\rm d}Y(t)]={\rm d}{\rm E}[Y(t)]=\left(\frac{\lambda}{\upsilon}\right){\rm d}t \end{equation} $ (11)

故期望退货量

$ \begin{equation} R={\rm E}\int_0^\tau\frac{\lambda}{\upsilon}{\rm d}t=\frac{\lambda}{\upsilon\theta} \end{equation} $ (12)

最后, 将式(9)两边同除$\varphi(\alpha){-}\theta$得:

$ \begin{equation} {\rm E}\int_0^\tau{\rm e}^{-\alpha W(t)}{\rm d}t=\frac{-{\rm e}^{-\alpha s}+ \alpha\Pi}{\varphi(\alpha)-\theta} \end{equation} $ (13)

对式(13)两边关于$\alpha$求导, 并令$\alpha=0$, 得期望库存量

$ \begin{equation} H=-\frac{{\rm d}}{{\rm d}\alpha}\Big({\rm E}\hspace{-1.7mm}\int_0^\tau\hspace{-2.0mm}{\rm e}^{-\alpha W(t)}{\rm d}t\Big)\Big|_{\alpha=0} =\frac{s{+}\Pi}{\theta}{+}\frac{\lambda{-}\mu\upsilon}{\upsilon\theta^2} \end{equation} $ (14)

下面, 给出对应$\tau(z)$的有关期望量.因该段时间内$W$总大于$0$, 故无反射发生, 即式(4)最后一项为0.采用上面类似方法, 可得:

$ \begin{equation} {\rm E}\int_0^{\tau(z)}{\rm e}^{-\alpha W(t)}{\rm d}t=\frac{{\rm e}^{-\alpha s}-{\rm e}^{-\alpha z}}{\varphi(\alpha)} \end{equation} $ (15)

在此基础上, 令式(15)中$\alpha{\rightarrow}0$, 并对右边应用两次L.Hospital法则, 可得期望时间

$ \begin{equation} {\rm E}[\tau(z)]=\frac{\upsilon(z-s)}{\mu\upsilon-\lambda} \end{equation} $ (16)

据式(11)和式(16), 求得期望退货量

$ \begin{equation} R(z)={\rm E}\int_0^{\tau(z)} {\rm d}Y(t)=\frac{\lambda(z-s)}{\mu\upsilon-\lambda} \end{equation} $ (17)

对式(15)两边关于$\alpha$求导, 并令$\alpha=0$, 得期望库存量

$ \begin{equation} H(z) = {\rm E} \int_0^{\tau(z)}\hspace{-3.5mm} W(t){\rm d}(t) = \frac{(z^2{-}s^2)(\mu\upsilon{-}\lambda)\upsilon{+}2\lambda(z{-}s)}{2(\mu\upsilon-\lambda)^2} \end{equation} $ (18)
3.2 确定库存水平分布函数

易知, 时间$\tau$结束时的库存水平$W(\tau-)$ (因补货引发库存跳跃, 用左极限表示补货前水平)为随机变量.故求期望补货量, 需确定$W(\tau-)$的分布函数.为此, 构建一个从水平$s$开始, 周期为$\tau$的辅助更新过程$W_1=(W_1(t))_{t{\geq}0}$.当$0\leq t < \tau$时, 有$W_{1}(\hspace{-0.0mm}t\hspace{-0.0mm})\hspace{-1.0mm}=\hspace{-1.0mm}W(\hspace{-0.0mm}t\hspace{-0.0mm})$.下面, 先求$W_{1}$的平稳分布函数.

引理2.过程$W_1$的平稳密度函数为

$ \begin{equation} P_{0}=\frac{\mu(A_1{+}A_2)}{\lambda+\theta} \end{equation} $ (19)
$ \begin{equation} f_1(x)=A_1{\rm e}^{\alpha_1(\theta)x}{+}A_2{\rm e}^{\alpha_2(\theta)x}, \quad 0 < x\leq s \end{equation} $ (20)
$\begin{equation} f_2(x)=A_3{\rm e}^{\alpha_2(\theta)x}, \quad s <x \end{equation} $ (21)

这里, $\alpha_1(\theta)$$\alpha_2(\theta)$分别由式(5)和(6)给出. $A_1$$A_2$$A_3$为常数, 其值在证明中给出.

证明.根据水平穿越理论[23]知, 库存水平的稳态密度等于其上穿(或下穿)水平$x$的长程平均数.故利用PASTA原则, 通过使库存水平的系统点(简称系统点)轨迹, 上穿和下穿某一水平速率相等, 构建$x$在区间$(0, s]$内的速率平衡方程

$ \begin{equation*} \mu f_1(x)=\lambda\int_0^x{\rm e}^{-\upsilon(x-y)}f_1(y){\rm d}y+ \theta\int_0^xf_1(y){\rm d}y+ \end{equation*} $
$ \begin{equation} \lambda{\rm e}^{-\upsilon x}P_{0}+\theta P_{0} \end{equation} $ (22)

上式左边为需求引起系统点下穿水平$x$的速率.右边第$1$项是量大于$x{-}y$的退货到达, 引起位于$y\in(0, x)$系统点上穿水平$x$的速率; 第$2$项是供应商中断结束, 引起位于$y\in(0, x)$系统点上穿水平$x$的速率; 第$3$项和第$4$项分别是量大于$x$的退货到达和供应商中断结束, 引发位于水平$0$上的系统点上穿水平$x$的速率.

为求$f_1(x)$, 需将式(22)化为微分方程.定义微分算子$\langle D\rangle={\rm d}/{\rm d}x$.对式(22)运用$\langle D\rangle\langle D{+}\upsilon\rangle$得:

$ \begin{equation} \mu \ddot{f}_1(x){-}(\lambda{+}\theta{-}\upsilon\mu)\dot{f}_1(x){-}\upsilon\theta f_1(x)=0 \end{equation} $ (23)

式(23)对应的特征方程为

$ \begin{equation} \mu\omega^2{-}(\lambda{+}\theta{-}\upsilon\mu)\omega{-}\theta\upsilon=0 \end{equation} $ (24)

求解式(24)得: $\omega_1=\alpha_1(\theta)$$\omega_2=\alpha_2(\theta)$.

所以, 式(23)的通解为

$ \begin{equation} f_1(x)=A_1{\rm e}^{\alpha_1(\theta)x}+A_2{\rm e}^{\alpha_2(\theta)x} \end{equation} $ (25)

这里, $A_1$$A_2$为待定常数.

同理, 构建$x$在区间$(s, \infty)$的速率平衡方程

$\begin{equation*} \mu f_2(x)+\theta\int_x^\infty f_2(y){\rm d}y=\lambda\int_0^s{\rm e}^{-\upsilon(x-y)}f_1(y){\rm d}y+ \end{equation*} $
$ \begin{equation} \lambda\int_s^x{\rm e}^{-\upsilon(x-y)}f_2(y){\rm d}y+\lambda P_{0}{\rm e}^{-\upsilon x} \end{equation} $ (26)

式(26)各项说明参考式(22), 故不再赘述.

下面确定$f_2(x)$.因退货到达时间和批量均服从指数分布, 故据指数分布的无记忆性知, 系统点在水平$x$之上逗留时间$\tau_x$ (即系统点从上穿到下穿水平$x$间的时间)的分布独立于$x$, 且${\rm E}[\tau_x]$为常数, 设它等于$\bar{\Gamma}_x$.另外, 若$D_t(x)$表示时间$(0, t)$内系统点下穿水平$x$的次数, 那么, 它是一个更新计数过程.令$d_x$表示系统点相继两次下穿水平$x$的时间.则过程$\{D_t(x), t{\geq}0\}$的更新速率为

$ \begin{equation} \lim\limits_{t\rightarrow\infty}\frac{D_t(x)}{t} =\lim\limits_{t\rightarrow\infty}\frac{{\rm E}[D_t(x)]}{t}=\mu f_2(x)=\frac{1}{{\rm E}[d_x]} \end{equation} $ (27)

据更新过程理论, 系统点位于水平$x$之上的时间比例有下列关系

$ \begin{equation} \frac{{\rm E}[\tau_x]}{{\rm E}[d_x]}=\frac{\bar{\Gamma}_x}{{\rm E}[d_x]}=1-F_2(x) \end{equation} $ (28)

这里, $F_2(x)=\int_0^xf_2(y){\rm d}y$.

由式(27)和(28)得:

$ \begin{equation} \mu\bar{\Gamma}_xf_2(x)=1-F_2(x) \end{equation} $ (29)

故据式(29)知, $f_2(x)$可表示为

$ \begin{equation} f_2(x)=A_3{\rm e}^{-\eta x} \end{equation} $ (30)

这里, $A_3$$\eta$均为待定正常数.

在此基础上, 确定$P_0$及常数$A_1$$A_2$$A_3$$\eta$.

首先, 利用式(22)和(25), 令$x$从右边趋于$0$得:

$ \begin{equation} P_{0}=\frac{\mu(A_1+A_2)}{\lambda+\theta} \end{equation} $ (31)

将式(25)、式(30)和式(31)代入式(26), 并比较等式两边$x$函数的系数得:

$ \begin{equation} \mu+\frac{\theta}{\eta}=\frac{\lambda}{\upsilon-\eta} \end{equation} $ (32)

$ \begin{align} &\Big[\frac{\mu}{\lambda+\theta}{-} \frac{1{-}{\rm e}^{(\upsilon{+}\alpha_1(\theta))s}}{\upsilon{+} \alpha_1(\theta)}\Big]A_1+\notag\\ &\qquad \Big[\frac{\mu}{\lambda+\theta}-\frac{1{-}{\rm e}^{(\upsilon{+}\alpha_2(\theta))s}}{\upsilon{+}\alpha_2(\theta)}\Big]A_2- \notag\\ &\qquad \frac{{\rm e}^{(\upsilon-\eta)s}}{\upsilon-\eta}A_3=0 \end{align} $ (33)

据式(32)得: $\eta=-\alpha_2(\theta)$.

其次, 分布函数要满足下列标准条件

$ \begin{equation} P_{0}+\int_0^sf_1(y){\rm d}y+\int_s^\infty f_2(y){\rm d}y=1 \end{equation} $ (34)

将式(25)、式(30)和式(31)代入式(34)得:

$ \begin{align} &\left[\frac{\mu}{\lambda\hspace{-1.0mm}+\hspace{-1.0mm}\theta}-\frac{1{-}{\rm e}^{\alpha_1(\theta) s}}{\alpha_1(\theta)}\right]A_{1}+ \left[\frac{\mu}{\lambda\hspace{-1.0mm}+\hspace{-1.0mm}\theta}{-}\frac{1{-}{\rm e}^{\alpha_2(\theta)s}}{\alpha_2(\theta)}\right]A_2-\notag\\ &\qquad \frac{{\rm e}^{\alpha_2(\theta)s}}{\alpha_2(\theta)}A_3=1 \end{align} $ (35)

另外, 在水平$s$上系统点还要满足速率平衡方程

$ \begin{equation} \mu f_2(s){+}\theta\int_s^\infty f_2(y){\rm d}y{+}\theta\int_0^s f_1(y){\rm d}y{+}P_0 \theta=\mu f_1(s) \end{equation} $ (36)

结合式$(34)$, 式$(36)$可化简为

$ \begin{equation} \mu f_1(s)-\mu f_2(s)=\theta \end{equation} $ (37)

将式$(25)$和式$(30)$代入式$(37)$得:

$ \begin{equation} \mu {\rm e}^{\alpha_1(\theta)s}A_1+\mu {\rm e}^{\alpha_2(\theta)s}A_2-\mu {\rm e}^{\alpha_2(\theta)s}A_3=\theta \end{equation} $ (38)

最后, 据式(33)、式(35)和式(38)构成的方程组得:

$ \begin{equation} [A_1, A_2, A_3]^{\rm T}={\rm \Lambda}^{-1}{\pmb \ell} \end{equation} $ (39)

这里, ${\pmb \ell}=[\theta, 1, 0]^{\rm T}, $

$ \begin{array}{l} \Lambda = \\ \left[ {\begin{array}{*{20}{c}} {\mu {{\rm{e}}^{{\alpha _1}(\theta )s}}}&{\mu {{\rm{e}}^{{\alpha _2}(\theta )s}}}&{ - \mu {{\rm{e}}^{{\alpha _2}(\theta )s}}}\\ {\frac{\mu }{{\lambda + \theta }} - \frac{{1 - {{\rm{e}}^{{\alpha _1}(\theta )s}}}}{{{\alpha _1}(\theta )}}}&{\frac{\mu }{{\lambda + \theta }} - \frac{{1 - {{\rm{e}}^{{\alpha _2}(\theta )s}}}}{{{\alpha _2}(\theta )}}}&{ - \frac{{{{\rm{e}}^{{\alpha _2}(\theta )s}}}}{{{\alpha _2}(\theta )}}}\\ {}&{}&{}\\ {\frac{\mu }{{\lambda + \theta }} - \frac{{1 - {{\rm{e}}^{(\upsilon + {\alpha _1}(\theta ))s}}}}{{\upsilon + {\alpha _1}(\theta )}}}&{\frac{\mu }{{\lambda + \theta }} - \frac{{1 - {{\rm{e}}^{(\upsilon + {\alpha _2}(\theta ))s}}}}{{\upsilon + {\alpha _2}(\theta )}}}&{ - \frac{{{{\rm{e}}^{(\upsilon + {\alpha _2}(\theta ))s}}}}{{\upsilon + {\alpha _2}(\theta )}}}\\ {}&{}&{} \end{array}} \right] \end{array} $

$\tau$服从指数分布, 故据PASTA原则得, $W_1$的平稳分布函数即为待求$W(\tau-)$的分布函数.

3.3 确定状态转移函数

对固定时间$t$, 定义转移函数矩阵$P(t)=[P_{ij}(t)]$, 其中, $P_{ij}(t)=P\{I(t)=j|I(0)=i\}$, $i, j=0, {\cdots}, 3$.考虑到$W$从任一水平$z(z>s)$出发首次到达水平$s$所需时间$\tau(z)$为随机的, 下面, 给出对应$\tau(z)$的状态转移函数期望矩阵$[\phi_{_{ij}}(z)]$, 即$[\phi_{_{ij}}(z)]={\rm E}[P(\tau(z))]$.

引理3.$\varpi_k=\gamma_k{+}\theta_k$, $k=1, 2$, 和$\varpi_3=\varpi_1{+}\varpi_2$, 则

$ \begin{equation} [\phi_{_{ij}}(z)]=U\bar{D}(z)U^{-1}, i, j=0, \cdots, 3 \end{equation} $ (40)

这里

$ \begin{align*} &U=\left[\begin{array}{cccc} 1&1&1&1\\ 1&1&-\frac{\theta_2}{\gamma_2}&-\frac{\theta_2}{\gamma_2}\\ 1&-\frac{\theta_1}{\gamma_1}&1&-\frac{\theta_1}{\gamma_1}\\ 1&-\frac{\theta_1}{\gamma_1}&-\frac{\theta_2}{\gamma_2}&\frac{\theta_1\theta_2}{\gamma_1\gamma_2}\\ \end{array}\right]\end{align*} $
$ \begin{align*} &\bar{D}(z)=\left[\begin{array}{cccc} 1&0&0&0\\ 0&{\rm e}^{-\bar{\alpha}_1(z{-}s)}&0&0\\ 0&0&{\rm e}^{-\bar{\alpha}_2(z{-}s)}&0\\ 0&0&0&{\rm e}^{-\bar{\alpha}_3(z{-}s)}\\ \end{array}\right]\end{align*} $
$ \begin{align*} &U^{-1}=\left[\begin{array}{cccc} \frac{\theta_1\theta_2}{\varpi_1\varpi_2}&\frac{\gamma_2\theta_1}{\varpi_1\varpi_2}&\frac{\gamma_1\theta_2}{\varpi_1\varpi_2}&\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}\\ \frac{\gamma_1\theta_2}{\varpi_1\varpi_2}&\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}&-\frac{\gamma_1\theta_2}{\varpi_1\varpi_2}&-\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}\\ \frac{\gamma_2\theta_1}{\varpi_1\varpi_2}&-\frac{\gamma_2\theta_1}{\varpi_1\varpi_2}&\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}&-\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}\\ \frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}&-\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}&-\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}&\frac{\gamma_1\gamma_2}{\varpi_1\varpi_2}\\ \end{array}\right] \end{align*} $

其中, 常数$\bar{\alpha}_i$据式$(5)$$\theta$等于$\varpi_i$得到, $i=1, 2, 3$.

证明. 对固定时间$t$, 文献[11]给出$4$个状态连续时间Markov链的转移函数矩阵

$ \begin{equation} P(t)=UD(t)U^{-1} \end{equation} $ (41)

这里

$ \begin{matrix} &&D(t)=\left[\begin{array}{cccc} 1&0&0&0\\ 0&{\rm e}^{-\varpi_1t}&0&0\\ 0&0&{\rm e}^{-\varpi_2t}&0\\ 0&0&0&{\rm e}^{-\varpi_3t}\\ \end{array}\right] \end{matrix} $

$P_j=\lim_{t\rightarrow\infty}P_{ij}(t)$, 可得$4$个状态的长程概率

$ \begin{equation} [P_{0}, P_{1}, P_{2}, P_{3}]=\frac{1}{\varpi_1\varpi_2}[\theta_1\theta_2, \theta_1\gamma_2, \gamma_1\theta_2, \gamma_1\gamma_2] \end{equation} $ (42)

另据式(41)得

$ \begin{equation} [\phi_{{ij}}(z)]=U{\rm E}[D(\tau(z))]U^{-1} \end{equation} $ (43)

下面, 基于Kella-Whitt鞅, 确定${\rm E}[D(\tau(z))]$.

首先, 利用式(4), 对$\tau(z)$应用最优抽样定理得:

$ \begin{align} &(\varphi(\alpha) -\beta){\rm E}\int_0^{\tau(z)}{\rm e}^{-\alpha W(t)-\beta t}{\rm d}t=\notag\\ &\qquad {\rm E}[{\rm e}^{-\alpha W(\tau(z))-\beta\tau(z)}]-{\rm e}^{-\alpha z}=\notag \\ &\qquad {\rm e}^{-\alpha s}{\rm E}[{\rm e}^{-\beta\tau(z)}]-{\rm e}^{-\alpha z}. \end{align} $ (44)

其次, 取$\beta$等于$\varpi_i$, 得$\varphi(\alpha)- \varpi_i=0$正根$\bar{\alpha}_i=\alpha_{\hspace{-0.2mm}1}( \varpi_i)$.

将它代入式$(44)$, 左边为$0$, 故有:

$ \begin{equation} {\rm E}[{\rm e}^{-\varpi_i\tau(z)}]={\rm e}^{-\bar{\alpha}_i(z-s)}, \hspace{4mm} i=1, 2, 3 \end{equation} $ (45)

据式$(45)$得:

$ \begin{equation} {\rm E}[D(\tau(z))]=\bar{D}(z) \end{equation} $ (46)

最后, 将式$(46)$代入式$(43)$, 即得式$(40)$.

4 确定循环的期望费用和时间

为简化表述, 令随机变量$Z=W(\hspace{-0.2mm}\tau{-})$, $\rho_{1}=P(Z < s)$$\rho_{2}={\rm E}[Z1_{\{Z{ < }s\}}]$$\rho_{3}={\rm E}[R(Z)1_{\{Z{>}s\}}]$$\rho_{4}={\rm E}[H(Z)1_{\{Z{>}s\}}]$$\psi_{ij}={\rm E}[\phi_{ij}(Z)1_{\{Z{>}s\}}]$, $i=1, 2$, $j=0, \cdots, 3$.

另外, 令$\bar{K}_i=K_i{+}k_iq_{i}$, $i=1, 2$, $\bar{K}_{0}=\bar{K}_{1}{+}\bar{K}_{2}$$\bar{\vartheta}_j=\vartheta_1\psi_{1j}{+}\vartheta_2\psi_{2j}$, $j=1, 2, 3, $以及

$ \begin{align*} M=\left[\begin{array}{ccc} 1-\phi_{{11}}(\bar{q}_{1})&-\phi_{{12}}(\bar{q}_{1})&-\phi_{{13}}(\bar{q}_{1})\\ -\phi_{{21}}(\bar{q}_{2})&1-\phi_{{22}}(\bar{q}_{2})&-\phi_{{23}}(\bar{q}_{2})\\ -\bar{\vartheta}_1-\rho_1\vartheta_1&-\bar{\vartheta}_2-\rho_1 \vartheta_2&1-\bar{\vartheta}_3\\ \end{array}\right] \end{align*} $

这里, $\vartheta_i$值已在第3.1节给出. $\rho_{1}$$\rho_{2}$$\rho_{3}$$\rho_{4}$$\psi_{{ij}}$值, 利用引理$2$的分布函数及式$(17)$$(18)$$(40)$易得.最后, 求得矩阵$M$.考虑空间限制, 省略其结果.

定理1.利用引理1、2和3给出的结果, 可求得循环期望费用

$ \begin{equation} C_{{00}}=\sigma_{0}+[\phi_{{01}}(\bar{q}_{0}), \phi_{{02}}(\bar{q}_{0}), \phi_{{03}}(\bar{q}_{0})]M^{-1}{\pmb \epsilon} \end{equation} $ (47)

这里, ${\pmb \epsilon}=[\sigma_1, \sigma_2, \bar{\sigma}]^{\rm T}$.其中, $\bar{\sigma}$$\sigma_i$ ($i=0, 1, 2$)为常数, 其值在证明中给出. $\phi_{0i}(\bar{q}_0)$据式(40)取$z$等于$\bar{q}_0$得到.

证明.由问题描述知, 过程$W$从更新点(即状态$i$$0$水平为$\bar{q}_{0}$)出发, 到达水平$s$时, 系统供应商的$4$种状态均可能发生, 其概率分别为$\phi_{{00}}(\bar{q}_{0})$$\phi_{{01}}(\bar{q}_{0})$$\phi_{{02}}(\bar{q}_{0})$$\phi_{{03}}(\bar{q}_{0})$.故据更新过程原理得:

$ \begin{align} C_{{00}}=\, &\phi_{{00}}(\bar{q}_{0}) [hH(\bar{q}_{0})+\delta R(\bar{q}_{0})+\bar{K}_{0}]+\notag\\ &\phi_{{01}}(\bar{q}_{0})[hH(\bar{q}_{0})+\delta R(\bar{q}_{0})+\bar{K}_{1}+C_{{10}}]+\notag\\ &\phi_{{02}}(\bar{q}_{0})[hH(\bar{q}_{0})+ \delta R(\bar{q}_{0})+\bar{K}_{2}+C_{{20}}]+\notag\\ &\phi_{{03}}(\bar{q}_{0})[hH(\bar{q}_{0}) +\delta R(\bar{q}_{0})+C_{{30}}]=\notag\\ &hH(\bar{q}_{0})+\delta R(\bar{q}_{0})+\sum\limits_{j=0}^2\phi_{{0j}} (\bar{q}_{0})\bar{K}_j+\notag\\ & \sum\limits_{j=1}^3\phi_{{0j}}(\bar{q}_{0})C_{{j0}} =\sigma_{0}+\sum\limits_{j=1}^3\phi_{{0j}}(\bar{q}_{0})C_{{j0}} \end{align} $ (48)

这里

$ \begin{equation*} \sigma_0=hH(\bar{q}_0)+\delta R(\bar{q}_0)+\sum\limits_{j=0}^2\phi_{0j}(\bar{q}_0)\bar{K}_j \end{equation*} $

据式$(17)$、式$(18)$和式$(40)$, 取$z$等于$\bar{q}_0$, 求得$R(\bar{q}_0)$$H(\bar{q}_0)$$\phi_{{0j}}(\bar{q}_0)$, 将它们代入上式, 即得$\sigma_0$.

下面, 分别确定式(48)中的$C_{i0}$, $i=1, 2, 3$.

首先, 针对$i=1$$2$.采用上面类似方法得

$ \begin{align} C_{{i0}}=&hH(\bar{q}_i)+ \delta R(\bar{q}_i)+\sum\limits_{j=0}^2 \phi_{{ij}}(\bar{q}_i)\bar{K}_j+\notag\\ &\sum\limits_{j=1}^3\phi_{{ij}}(\bar{q}_i)C_{{j0}} =\sigma_i+\sum\limits_{j=1}^3\phi_{{ij}}(\bar{q}_i)C_{{j0}} \end{align} $ (49)

这里

$ \begin{equation*} \sigma_i=hH(\bar{q}_i)+\delta R(\bar{q}_i)+\sum\limits_{j=0}^2\phi_{ij}(\bar{q}_i)\bar{K}_j \end{equation*} $

$\sigma_i$求解方法与$\sigma_0$相同, 故不再赘述.

其次, 对$i=3$.易知, $\tau$结束时库存水平$W ( \tau- )$, 即$Z$, 有两种情况: 1) $Z{<}s$, 此时, 若$S_1$中断先结束(概率为$\vartheta_1$), 补货量为$\bar{q}_{1}{-}Z$; 若$S_2$中断先结束(概率为$\vartheta_2$), 补货量为$\bar{q}_{2}{-}Z$. 2) $Z{>}s$, 此时, 供应商不再补货.当$W$从水平$Z$又降到水平$s$时, 系统供应商的$4$种状态均可能发生.两种情况下除包含库存、补货和退货费外, 还有短缺费.采用上面类似方法得:

$ \begin{align} &C_{{30}} = \pi\Pi + hH + \delta R + \vartheta_1\big\{{\rm E}[1_{\{Z < s\}}(K_{1} + k_1(\bar{q}_{1}-Z)+\notag\\ & C_{{10}})] + {\rm E}[1_{\{Z>s\}}(hH(Z) + \delta R(Z) + \sum\limits_{j=0}^2\phi_{{1j}}(Z)\bar{K}_j{+}\notag\\ &\sum\limits_{j=1}^3\phi_{{1j}}(Z)C_{{j0}})]\big\} + \vartheta_2\big\{{\rm E}[1_{\{Z < s\}}( K_{2} + k_2(\bar{q}_{2}-Z) +\notag\\ & C_{{20}})]{+}{\rm E}[1_{\{Z>s\}}(hH(Z) + \delta R(Z) + \sum\limits_{j=0}^2\phi_{{2j}}(Z) \bar{K}_j + \notag\\ & \sum\limits_{j=1}^3\phi_{{2j}} (Z)C_{{j0}})]\big\} = \bar{\sigma} + (\vartheta_{1}\rho_{1} + \vartheta_{1}\psi_{{11}} + \vartheta_{2}\psi_{{21}})C_{{10}}+\notag\\ & (\vartheta_{2}\rho_{1}+\vartheta_{1}\psi_{{12}}+ \vartheta_{2}\psi_{{22}})C_{{20}}+(\vartheta_{1}\psi_{{13}}+\vartheta_{2}\psi_{{23}})C_{{30}} \end{align} $ (50)

这里

$ \begin{equation*} \bar{\sigma}=\pi\Pi+\delta(R+\rho_{3})+h(H+\rho_{4})+\rho_{1}[\vartheta_1(\bar{K}_1 +k_1s)+ \end{equation*} $
$ \begin{equation*} \vartheta_2(\bar{K}_2+k_2s)]- \rho_{2}(\vartheta_{1}k_1+\vartheta_{2}k_2)+\sum\limits_{j=0}^2(\vartheta_{1}\psi_{{1j}}+ \vartheta_{_2}\psi_{{2j}})\bar{K}_j \end{equation*} $

将式(1)、式(2)和式(3)及上面求得的变量代入, 可得$\bar{\sigma}$.

据式$(49)$和式$(50)$构成的方程组, 求得

$ \begin{equation} [C_{{10}}, C_{{20}}, C_{{30}}]^{\rm T}=M^{-1}{\pmb \epsilon} \end{equation} $ (51)

最后, 将式$(51)$代入式$(48)$, 即得式$(47)$.

定理2.利用引理1 $\sim$ 3给出的结果, 可求得循环期望时间

$ \begin{equation} T_{00}={\rm E}[\tau(\bar{q}_{0})]+[\phi_{{01}}(\bar{q}_{0}), \phi_{{02}}(\bar{q}_{0}), \phi_{{03}}(\bar{q}_{0})] M^{-1}{\pmb \nu} \end{equation} $ (52)

这里, ${\pmb \nu}=\big[{\rm E}[\tau(\bar{q}_{1})], {\rm E}[\tau(\bar{q}_{2})], 1/\theta{+}{\rm E}[\tau(Z)1_{\{Z>s\}}]\big]^{\rm T}$.其中, ${\rm E}[\tau(\bar{q}_{i})]$可据式$(16)$$z$等于$\bar{q}_{i}$得到, $i=0, 1, 2$; ${\rm E}[\tau(Z)1_{\{Z>s\}}]$据式$(16)$和引理$2$的分布函数求得.考虑空间限制, 省略其结果.

证明.与确定期望费用$C_{{00}}$类似, 可求得:

$ T_{00}=\phi_{{00}}(\bar{q}_{0}){\rm E}[\tau(\bar{q}_{0})]+\phi_{{01}}(\bar{q}_{0})\big({\rm E}[\tau(\bar{q}_{0})]+T_{10}\big)+\\ \phi_{{02}}(\bar{q}_{0})\big({\rm E}[\tau(\bar{q}_{0})]+T_{20}\big)+\phi_{{03}}(\bar{q}_{0})\big({\rm E}[\tau(\bar{q}_{0})]+T_{30}\big)=\\ {\rm E}[\tau(\bar{q}_{0})]+\sum\limits_{j=1}^3\phi_{{0j}}(\bar{q}_{0})T_{{j0}} $ (53)

同理, 可求得:

$ \begin{equation} T_{{i0}}={\rm E}[\tau(\bar{q}_{i})]+\sum\limits_{j=1}^3\phi_{{ij}}(\bar{q}_{i})T_{{j0}}, i=1, 2 \end{equation} $ (54)

至于$T_{30}$, 参考$C_{30}$的确定方法, 可得:

$ \begin{align} T_{30}= &{\rm E}[\tau]+ \vartheta_{1}\big\{{\rm E}\hspace{-0.2mm}[T_{10} 1_{\{Z{ < }s\}}]+ {\rm E}[(\tau(Z)+\notag\\ &\sum\limits_{j=1}^3\phi_{{1j}}(Z)T_{j0}) 1_{\{Z{>}s\}}]\big\}+\vartheta_{2}\big\{{\rm E}[T_{20}1_{\{Z{ < }s\}}]+\notag\\ &{\rm E}[(\tau(Z)+\sum\limits_{j=1}^3\phi_{{2j}}(Z) T_{j0})1_{\{Z{>}s\}}] \big\} =\notag\\ &\dfrac{1}{\theta}+{\rm E}[\tau(Z)1_{\{Z>s\}}]+[\vartheta_{2}\psi_{{21}}+\notag\\ &\vartheta_{1}(\rho_{1}+\psi_{{11}})]T_{10} +[\vartheta_{1}\psi_{{12}}+\vartheta_{2}(\rho_{1}+\notag\\ & \psi_{{22}})]T_{20}+(\vartheta_{1}\psi_{{13}}+\vartheta_{2} \psi_{{23}})T_{30}. \end{align} $ (55)

据式$(54)$$(55)$得:

$ \begin{equation} [T_{10}, T_{20}, T_{30}]^{\rm T}=M^{-1}{\pmb \nu} \end{equation} $ (56)

最后, 将式$(56)$代入式$(53)$, 即得式$(52)$.

5 仿真试验与分析

利用定理1和2, 构建长程平均费用率优化模型

$ \begin{equation*} {\rm min} TC(q_1, q_2, s)=\frac{C_{00}}{T_{00}} \end{equation*} $

因无法判断该函数的单调性或凹凸性, 这里, 利用MATLAB非线性优化函数, 求得它的近似最优解(为便于叙述, 下面简称为最优解).

试验目的是分析退货存在时, 双源采购策略缓解中断对库存影响的效果, 以及供应商差异性和系统参数对最优订货策略$q_1^\ast$, $q_2^\ast$, $s^\ast$和费用$TC^\ast$的影响.因给定$S_i$的可靠性$\zeta_i=\frac{1/\gamma_i}{1/\gamma_i{+}1/\theta_i}$, 中断又有频率低(高)持续时间长(短)两种类型.为全面分析它的影响, 据$S_1$$S_2$可靠性和中断类型不同的组合形式, 表 1给出$8$组试验数据, 其中, 前$4$组为不同可靠性组合; 后$4$组为可靠性固定条件下不同中断类型组合.

表 1 供应商的状态参数及可靠性 Table 1 Status parameters and reliability of the suppliers

试验采用其他参数的基本值如下: $K_1=10$, $k_1=1$, $K_2=20$, $k_2=2$, $\mu=120$, $\lambda=15$, $\upsilon=0.5$, $h=0.3$, $\pi=15$, $\delta=5$.

5.1 双源采购策略的效果

针对表 1每组实验数据, 在其他参数取基本值条件下, 表 2给出零售商单源和双源采购对应的最优结果.这里, 单源采购的结果, 是假设另一供应商的中断时间分布参数为$0$条件下得到的.

表 2 单源和双源采购对应的最优控制策略和费用 Table 2 Optimal control policy and cost for single and dual sourcing

表 2知, 供应商的可靠性和中断类型, 均对系统产生较大影响.其表现为, 尽管$S_2$的补货费用比$S_1$的大, 但当它的可靠性比$S_1$高(如第$3$组), 或它的中断类型为频率高持续时间短, 而$S_1$的中断类型为频率低持续时间长(如第$6$组)时, 零售商从它比从$S_1$单源采购更经济.

另外, 针对每组试验数据, 双源采购的费用$TC^\ast$均比单源采购的费用$\overrightarrow{TC}^\ast$$\widehat{TC}^\ast$低, 尤其是, 在$S_1$$S_2$的可靠性低(如第$4$组)或中断类型均为频率低持续时间长(如第$5$组)情况下, 双源采购比单源采购的优势更明显.究其原因是, 在两个供应商可靠性低或中断类型均为频率低持续时间长情况下, 若采用单源采购策略, 为避免中断引发的库存短缺, 零售商不但要订购大量产品, 而且还要保持较高订货水平, 由此造成过多的库存费用.若采用双源采购策略, 就能起到较好地防御中断效果.这可根据双源采购下两供应商同时处入中断状态的长程概率$\frac{\gamma_1\gamma_2}{(\gamma_1{+}\theta_1)(\gamma_2{+}\theta_2)}$ (即式$(42)$$P_3$), 总小于单源采购下供应商处入中断状态的长程概率${\gamma_i}/{(\gamma_i{+}\theta_i)}, i=1, 2$, 得到证实.

5.2 供应商差异性对系统的影响

与供应商有关除状态参数外, 还有补货费用.在假设其他参数取基本值条件下, 针对每组实验数据, 表 3给出$S_2$相对$S_1$的固定和可变补货费用变化, 对应的最优结果.

表 3 $K_2$$k_2$变化对应的最优控制策略和费用 Table 3 Optimal control policy and cost for varying $K_2$ and $k_2$

首先, 据表 3知, 与人们预期相同, 当$k_2$$1$ (即$S_2$$S_1$的可变补货费相同)时, 增加$K_2$, 导致每次从$S_2$的固定订货费增大.为缓解$TC^\ast$大幅提高, 零售商采取增加订货量$q_1^\ast$$q_2^\ast$及减小订货水平$s^\ast$的策略, 从而降低从$S_2$订货的频率.当$K_2$$10$ (即$S_2$$S_1$的固定补货费相同)时, 增大$k_2$, 引起$S_2$订货费用显著增加, 为此, 零售商降低从$S_2$的订货量$q_2^\ast$和订货水平$s^\ast$, 相应增加从$S_1$的订货量$q_1^\ast$, 结果$TC^\ast$增大了.

其次, 表 3给出, 在$K_2$$10$$k_2$$1$, 即$S_2$$S_1$补货费用相同情况下, 可靠性高的供应商, 其补货量大(如第$2$$3$组); 若可靠性也相同, 则中断类型为频率高持续时间短的供应商, 具有较高的补货量(如第$6$$7$组); 若补货费用、可靠性和中断类型都相同, 即$S_2$$S_1$无差异, 此时, 它们的补货量才相同(如第$1$$4$$5$$8$组).

最后, 对比表 2表 3知, 无论$S_1$$S_2$的可靠性、中断类型及补货费用如何变化, $TC^\ast$均小于$\overrightarrow{TC}^\ast$$\widehat{TC}^\ast$.说明不管供应商差异性多大, 双源采购策略均优于单源采购策略.因此, 本文将文献[8-10]的结论, 推广到包含退货的更一般环境.

5.3 系统关键参数的灵敏度分析

其他参数取基本值条件下, 针对每组实验数据, 表 4给出库存费$h$和短缺费$\pi$变化, 对应的最优结果.

表 4 $h$$\pi$变化对应的最优控制策略和费用 Table 4 Optimal control policy and cost for varying $h$ and $\pi$

对比表 2表 4知, 与人们期望相同, 增大$h$, 引起$q_1^\ast$$q_2^\ast$$s^\ast$减小; 而增大$\pi$, 又造成$q_1^\ast$$q_2^\ast$$s^\ast$增加.但作用的主要对象不同.因$S_1$的补货费用较低, 零售商从其订货量较多, 故增大$h$, 造成$q_1^\ast$较大幅度减少; 而增大$\pi$, 主要影响短缺费用, 故$s^\ast$有较大幅度的增加.但它们增大均导致$TC^\ast$增加.

最后, 因退货为不可控的, 故退货费$\delta$变化, 仅与$TC^\ast$有关, 对$q_1^\ast$$q_2^\ast$$s^\ast$无影响.需求率$\mu$、退货批次$\lambda$和批量$1/\upsilon$对系统的影响, 符合人们的预期, 即当$\mu$增大时, $q_1^\ast$$q_2^\ast$$s^\ast$$TC^\ast$均增加; 当$\lambda$$1/\upsilon$增大时, $q_1^\ast$$q_2^\ast$$s^\ast$均减小, 而$TC^\ast$增加(因退货费$\delta$大于库存费$h$).考虑空间限制, 省略其计算结果.

6 结论

供应中断和退货环境下库存控制问题, 因包含多个随机因素(如供应和中断的持续时间及退货的批次和批量), 致使利用传统建模方法很难解决, 从而造成该问题在理论上缺乏深入研究.

在将库存水平变化表示为跳跃-扩散过程条件下, 利用水平穿越理论和鞅方法, 研究了双源采购库存控制模型的构建, 这也是本文的创新之处.在此基础上, 得到新环境下有关库存控制结论如下: 1)系统最优控制策略和费用, 不但与供应商可靠性有关, 而且还与中断类型有关. 2)当单位产品短缺费用较高时, 双源采购优于单源采购, 尤其是, 当两个供应商的可靠性低或中断类型均为频率低持续时间长时. 3)除供应中断外, 供应商的差异性及关键参数, 均会对最优控制策略和系统费用产生较大影响.

还可拓展研究如下: 1)本文假设中断是独立的, 但有时并非如此, 如暴风雪灾会导致两供应商同时中断.因此, 研究中断相关情况下双源采购问题更有实际意义; 2)本文假设需求和退货参数是固定的, 实际中它们很可能变化, 如顾客得知零售商供应中断时, 会增加需求和减少退货.因此, 研究这些参数与中断状态相关环境下库存问题更有使用价值.通过深入研究, 以期更好解决该问题.

参考文献
1
How to reduce the retail shortage rate[Online], available: http://www.chinawuliu.com.cn, May 6, 2016.
(如何降低零售业的缺货率[Online], 下载地址: http://www.chinawuliu.com.cn, May 6, 2016.)
2
To cope with high return rates vipshop contrarian price logistics[Online], available: http://www.chinawuliu.com.cn, May 6, 2016.
(应对高退货率唯品会逆势提价物流[Online], 下载地址: http://www.chinawuliu.com.cn, May 6, 2016.)
3
Snyder L V, Atan Z, Peng P, Rong Y, Schmitt A J, Sinsoysal B. OR/MS models for supply chain disruptions:a review. ⅡE Transactions, 2016, 48(2): 89-109.
4
Gerchak Y, Parlar M. Yield randomness, cost tradeoffs, and diversification in the EOQ model. Naval Research Logistics, 1990, 37(3): 341-354. DOI:10.1002/nav.v37.3
5
Tomlin B. On the value of mitigation and contingency strategies for managing supply chain disruption risks. Management Science, 2006, 52(5): 639-657. DOI:10.1287/mnsc.1060.0515
6
Qi L, Lee K. Supply chain risk mitigations with expedited shipping. Omega, 2015, 57(12): 98-113.
7
Tian Jun, Zhang Hai-Qing, Wang Ying-Luo. Emergency supplies purchasing model based on capacity option contract with dual purchasing sources. Systems Engineering-Theory & Practice, 2013, 33(9): 2212-2219.
( 田军, 张海青, 汪应洛. 基于能力期权契约的双源应急物资采购模型. 系统工程理论与实践, 2013, 33(9): 2212-2219. DOI:10.12011/1000-6788(2013)9-2212)
8
Silbermayr L, Minner S. A multiple sourcing inventory model under disruption risk. International Journal of Production Economics, 2014, 149(1): 37-46.
9
Zhu S X. Analysis of dual sourcing strategies under supply disruptions. International Journal of Production Economics, 2015, 170(1): 191-203.
10
Silbermayr L, Minner S. Dual sourcing under disruption risk and cost improvement through learning. European Journal of Operational Research, 2016, 250(1): 226-238. DOI:10.1016/j.ejor.2015.09.017
11
Parlar M, Perry D. Inventory models of future supply uncertainty with single and multiple suppliers. Naval Research Logistics, 1996, 43(2): 191-210. DOI:10.1002/(ISSN)1520-6750
12
Gürler Ü, Parlar M. An inventory problem with two randomly available suppliers. Operations Research, 1997, 45(6): 904-918. DOI:10.1287/opre.45.6.904
13
Harrison J M, Sellke T M, Taylor A J. Impulse control of Brownian motion. Mathematics of Operations Research, 1983, 8(3): 454-466. DOI:10.1287/moor.8.3.454
14
Dai J G, Yao D C. Brownian inventory models with convex holding cost, Part 1:average-optimal controls. Stochastic Systems, 2013, 3(2): 442-499. DOI:10.1287/11-SSY041
15
Yao D C, Chao X L, Wu J C. Optimal control policy for a Brownian inventory system with concave ordering cost. Journal of Applied Probability, 2015, 52(4): 909-925. DOI:10.1239/jap/1450802743
16
Benkherouf L, Bensoussan A. Optimality of an (s, S) policy with compound Poisson and diffusion demands:a quasi-variational inequalities approach. SIAM Journal on Control and Optimization, 2009, 48(2): 756-762. DOI:10.1137/080715883
17
Li X D, Tang D, Wang Y J, Yang X W. Optimal processing rate and buffer size of a jump-diffusion processing system. Annals of Operations Research, 2014, 217(1): 319-335. DOI:10.1007/s10479-013-1521-2
18
Song Chun-Yue, Li Ping. Numerical solution for optimal production control of unreliable production systems with diffusion terms. Control Theory & Applications, 2009, 26(7): 709-714.
( 宋春跃, 李平. 含扩散项不可靠生产系统最优生产控制的数值求解. 控制理论与应用, 2009, 26(7): 709-714.)
19
Muthuraman K, Seshadri S, Wu Q. Inventory management with stochastic lead times. Mathematics of Operations Research, 2015, 40(2): 302-327. DOI:10.1287/moor.2014.0671
20
Lou Shan-Zuo, Tian Xin-Cheng. Modeling and control for inventory with stochastic supply disruptions and returns. Acta Automatica Sinica, 2014, 40(11): 2436-2444.
( 娄山佐, 田新诚. 随机供应中断和退货环境下库存问题的建模与控制. 自动化学报, 2014, 40(11): 2436-2444.)
21
Lou Shan-Zuo, Tian Xin-Cheng. Contingent control of inventory under stochastic supply disruptions and returns. Acta Automatica Sinica, 2015, 41(1): 94-103.
( 娄山佐, 田新诚. 随机供应中断和退货环境下库存的应急控制. 自动化学报, 2015, 41(1): 94-103.)
22
Kella O, Boxma O. Useful martingales for stochastic storage processes with Lëvy-type input. Journal of Applied Probability, 2013, 50(2): 439-449. DOI:10.1239/jap/1371648952
23
Brill P H. Level Crossing Methods in Stochastic Models. New York, USA: Springer, 2008.