基于模型参考预测控制的时滞闭环供应链仿真分析

郭海峰，李 白

(沈阳理工大学 自动化与电气工程学院，沈阳 110159)

电子商务环境下，正向销售物流存在电子商务直销与传统线下销售两种渠道，随着第三方回收商的加入，第三方逆向物流与正向销售物流形成的双渠道闭环供应链使得供应链结构更加复杂，管理难度增加，而供应链中各环节的随机时滞增加了整个供应链系统的不确定性，这种系统网络的复杂性和不确定性使得牛鞭效应进一步加大，牛鞭效应的抑制问题更加困难。第三方逆向物流能够使企业集中精力于核心业务上，为企业节约物流成本，并且能够加快退货处理和废品回收，提升企业对顾客的水平服务，提高顾客满意度。陈朝芬等提出了基于循环经济的第三方回收物流模型[1]。Lu等研究了两级回收再制造供应链的选址问题[2]。费威对由制造商、零售商、第三方回收以及零供回收一体化四种回收模式下的利润进行了分析[3]。李凌乐等分析了第三方物流在以第三方物流为主导的物流服务供应链中的战略地位，并指出必须要有效地共享链条信息，从库存控制入手，降低仓储成本[4]。张成堂等针对供应商、第三方物流服务商与零售商组成的3层供应链系统，建立独立与联合决策模型，结果表明联合决策下系统利润达到最优，并采用Nash协商收益分配机制激励第三方物流服务商提高物流服务水平[5]。孙小清等在市场需求不确定情况下，针对一个上游物流服务需求商和下游3PLE构成的两级供应链模型，运用博弈论Bertrand模型对第三方物流企业进行定价策略分析，指出信息共享和企业联盟策略对第三方物流产业具有积极的影响[6]。梅中亚针对双渠道闭环供应链的问题，考虑第三方物流企业同时提供正向物流和逆向回收物流服务，对比和分析了集中决策和分散决策中供应链的利润情况，得出两种模式下利润的差异[7]。Savaskan等研究了闭环供应链三种回收模式下的定价与渠道选择问题，运用博弈理论得出零售商回收模式时的最优定价[8]。Kaya等研究了制造商对回收商的激励和最优生产决策，设计了基于转移支付的线性合同协调机制[9]。杨天剑等研究了存在价格竞争下的第三方回收平台定价策略，使供应链在“成本共担+二步定价”契约下实现完美协调[10]。肖冰将作业成本法应用于第三方物流企业中，给出了第三方物流企业在应用作业成本法进行成本核算的过程、方法等[11]。郭海峰等提出了一种可以有效处理复杂非线性闭环供应链网络动态系统的模型预测控制算法，使系统得到良好的动态性能，抑制其中的牛鞭效应[12]。周浩然等建立了引入第三方物流后的药品供应链系统动力学模型，并通过敏感性分析对模型进行优化改进，该方法可以降低库存水平，有效提高供应链运作效率[13]。庞燕等基于第三方物流配送模式，构建了一种开放式车辆路径问题数学模型，其优化目标为总行驶距离最短和车辆数最少，并设计了一个改进的两阶段禁忌搜索算法进行求解[14]。代应等在新零售背景下，针对第三方物流及零售商对市场需求预测不同的情形，建立了物流服务水平影响市场需求的Stackelberg模型，分析了物流服务水平影响市场占有率及市场总容量的契约模型[15]。关博等针对C2C电子商务第三方退换货逆向物流建立0-1MILP数学模型并进行数据模拟得出成本最小化时的设施选址情况和运输量，绘制物流分配规划图，该方法可以有效降低成本并提高效率[16]。综上所述，现有研究是从物流模型、选址优化、回收利润、定价策略和激励机制等角度来考虑第三方逆向物流，对第三方逆向物流的双渠道时滞闭环供应链的库存控制策略研究相对较少。

本文在相关研究的基础上，研究了第三方回收模式下的双渠道时滞闭环供应链，主要工作有：(1)考虑到各个节点时滞，针对第三方回收模式下的双渠道闭环供应链，建立闭环供应链库存平衡模型；在此基础上通过状态重组，建立增广状态空间模型，解决状态之间的复杂联系；(2)建立供应链系统的预测模型，推导状态估计和预测的一般公式，给出基于模型预测控制方法的生产/订购策略求解方法和模型预测控制KWIK算法的一般步骤；(3)通过建立SIMULINK仿真，讨论该方法对第三方回收模式下的双渠道时滞闭环供应链中生产、订购和库存波动的平滑作用，以及对牛鞭效应的抑制作用。

1 第三方逆向物流的双渠道时滞闭环供应链动态模型

按文献[17]建立第三方逆向物流的双渠道时滞闭环供应链如图1所示。

图1 第三方逆向物流的双渠道时滞闭环供应链

正向渠道包括生产商直销渠道和零售商分销渠道；逆向渠道是消费者通过第三方回收商退货给生产商，这里假设所有售出产品质量相同。其中α为电子商务渠道占有率，β为退货率；将消费者持有的商品当作虚拟库存，考虑生产和运输提前期，建立生产商、零售商、回收商和消费者库存平衡方程如下：

x3(k+1)=x3(k)+u2(k-τ1(k))-u1(k)-αd(k)+x1(k-τ2(k))

(1)

x2(k+1)=x2(k)+u1(k-τ3(k))-(1-α)d(k)

(2)

x1(k+1)=βx0(k-τ4(k))

(3)

x0(k+1)=x0(k)+(1-α)d(k-τ5(k))+αd(k-τ6(k))-βx0(k)

(4)

式中：xi(k)，i=0、1、2、3分别表示消费者、回收商、零售商和生产商在时刻k的库存水平；uj(k)，j=1、2分别代表零售商向生产商订货数量和生产商生产产品的数量；d(k)为在k时刻消费者的需求量；τm(k)，m=1、2、…、6分别代表生产商生产的、由回收商到生产商的、由生产商到零售商的、由消费者退货到回收商的、由零售商运输到消费者、生产商运输到消费者的生产/运输时滞。系统动态可写成如下矩阵形式。

(5)

式中：x(k)=(x2(k)x1(k)x0(k))T，

u(k)=(u2(k)u1(k)d(k))T，

τ=max(τ1，τ2，…，τ6)，

考虑提前期，k时刻的状态与k-1时刻的状态有关，通过状态重组建立增广状态空间模型，解决状态之间的复杂联系，这更具有鲁棒性，令

X(k)=(x(k)Tx(k-1)Tx(k-2)T…x(k-τ)Tu(k-1)Tu(k-2)T…u(k-τ)T)T

(6)

则增广状态空间方程描述为

X(k+1)=AX(k)+Bu(k)

(7)

Y(k)=CX(k)

(8)

式中：

Y(k)表示时刻k各个节点的库存量。

2 SIMULINK环境下模型预测控制方法

模型预测控制(MPC)是一种基于模型预测、滚动优化、矫正的过程控制算法，具有良好的跟踪性能和较强的鲁棒性。建立如图2所示的SIMULINK环境下第三方逆向物流的双渠道时滞闭环供应链模型预测控制结构图，其中包含参考输入、需求模型、MPC控制器和被控对象模型4个部分。

图2 第三方逆向物流的双渠道时滞闭环供应链模型预测控制结构图

2.1 模型预测控制方法

2.1.1 预测模型

在X(k)可测时，模型(7)可直接作为预测模型使用。设从k时刻起，输入发生Nc步变化，而后保持不变。由模型(7)可以预测出在u(k)，u(k+1)，…，u(k+Nc-1)作用下未来预测时域Np(Np>Nc)个时刻的系统状态。由模型(8)可以推导预测输出，则预测模型可以表示为以下矩阵形式。

(9)

式中：

X(k+i|k)表示在k时刻对k+i时刻库存状态的预测值；X(k|k)=X(k)为当前时刻各个节点的库存状态；Nc为控制时域；Np为预测时域。

2.1.2 状态估计

当闭环供应链系统的状态不可直接观测时，可使用预测模型(9)来进行估计，采用卡尔曼滤波器作为最优状态估计器，如式(10)和式(11)所示。

(10)

(11)

2.1.3 目标函数

在k时刻的优化目标是：使供应链系统各个节点的库存状态尽可能接近目标库存状态，希望产能和订货量不要剧烈变化，因此k时刻可以采用式(12)作为优化性能指标。

(12)

(13)

式中：G=2DpTQDp+2R，

gT=[2DpTQT(CpX(k)-R(k))]T

2.2 SIMULINK环境下KWIK二次优化算法

2.2.1 KWIK二次优化算法

二次规划问题(13)可以按照文献[18]提出的KWIK算法求解，二次规划问题(13)的K-T条件如下。

(14)

当Hessian矩阵G正定且M列满秩，式(14)的解可以表示为

(15)

式中：

(16)

由于G正定，可以由Cholesky因子表示，分解为G=LL-1，于是

(17)

将式(17)继续分解，可由矩阵T表示为以下形式

(18)

由式(17)和式(18)，式(16)可化简为

(19)

如果正确地确定了有效集，可以将式(19)代入式(15)来求得二次规划问题的最优解u*和μ*。

设uL为二次规划问题在当前有效集SL下的最优解，则

(20)

(21)

如果约束条件的集合为Ω，uL+1为二次优化问题在有效集SL+1=SL∪ΩKNEXT下的最优解，KNEXT表示添加到SL上的约束索引。uL+1也是有效集SL上的可行点，因此可以得到

(22)

g=-LLTuL+1+MLμL+1+aKNEXTt

(23)

式中为了确保可行性，取t=μKNEXT。

将式(22)和式(23)代入式(20)可以得到

(24)

将式(22)和式(23)代入式(21)可以得到

(25)

因此，随着约束条件的加入，有如下递归公式。

(26)

(27)

2.2.2 KWIK二次优化算法步骤

算法中NACT表示有效约束的个数，KNEXT表示加入SL的约束索引，KDROP表示从SL中删除的约束索引。

步骤0.求无约束时的最小值。

u=-G-1g=-L-TL-1g，令T1=L-T，NACT=0。

步骤1.检查是否违反了无效的不等式约束。如果所有的约束都满足，当前的解u是可行解并且是最优解，结束。否则，将KNEXT设置为最违反约束的索引，μ+=[μ0]T。

步骤3.计算步长t

(1)计算在对偶空间中不违背对偶可行性前提下的最大步长t1。如果NACT=0或者r≤0，那么t1=；否则，t1=min{μj/rj：rj>0}，KDROP=j。

(2)计算在原始空间中使得约束KNEXT可行的最小步长t2。如果|z|=0，那么t2=，否则

(3)步长t=min(t1，t2)。

步骤4.迭代计算

(1)初始空间和对偶空间均不可行：t=，二次规划问题没有可行解，结束。

(2)在对偶空间中搜索：如果t2=，t1为有限值，那么t=t1，则μ+=μ++t[-r1]T，删除约束KDROP，更新TL和WL，令NACT=NACT-1，转步骤2。

(3)在初始空间和对偶空间中搜索：u=u+tz，μ+=μ++t[-r1]T，如果t=t2，则在有效集SL上加入约束KNEXT，更新TL和WL，令NACT=NACT+1，转步骤1；如果t=t1，删除约束KDROP，更新TL和WL，令NACT=NACT-1，转步骤2。

2.3 牛鞭效应

牛鞭效应是评价供应链性能的最重要的绩效指标，由文献[17]可知，牛鞭效应等于供应链上游节点接收到需求量的方差与下游节点需求量的方差之比

(28)

3 仿真分析

设生产商库存标称值为300，零售商库存标称值为200，回收商库存标称值为0，消费者持有量标称值为50，生产商初始库存为200，零售商初始库存为120，回收商初始库存为0，消费者持有量为60，生产商生产量的约束为：0≤u1(k)≤100，零售商订购量的约束为：0≤u2(k)≤40，电子商务占有率α=0.3，退货率β=0.2，u1(k)、u2(k)的权值变化率为0.2，输出权值为1、0.8、0、0.1，控制时域为3，预测时域为15。这里考虑两种消费者需求模式：第一种消费者需求为随机扰动，为一随机数，取值范围为0≤d(k)≤25；第二种消费者需求为正弦函数，d(k)=sink+25，采样时间为1。

消费者需求波动分别如图3和图4所示。

图3 消费者随机扰动需求量

图4 消费者正弦函数需求量

通过模型预测控制器不断的滚动优化，在10～15个生产周期后生产商、零售商、回收商和消费者的库存水平趋近于设定的标称值，生产商、零售商、回收商和消费者在随机需求和正弦需求下的库存量波动曲线如图5和图6所示。

图6 正弦需求下的库存水平

图5 随机需求下的库存水平

在MPC控制下生产商生产量和零售商订货量经过初期扰动，大约10个周期后逐渐趋于平滑，接近标称值。生产商生产量和零售商订货量在随机需求和正弦需求下的波动曲线如图7和图8所示。

图7 随机需求下的生产/订购量

图8 正弦需求下的生产/订购量

牛鞭效应始终存在于供应链系统中，无法完全消除，但是随着模型预测控制器的不断优化，牛鞭效应逐渐被抑制。系统的牛鞭效应在随机需求和正弦需求下的波动曲线如图9和图10所示。

图9 随机需求下的牛鞭效应

图10 正弦需求下的牛鞭效应

文献[17]研究了Internet下双渠道闭环供应链，在忽略各环节生产/运输时滞和需求波动对输出影响的前提下，可以在保证系统稳定的前提下得到良好的动态性能，使系统中的牛鞭效应得到抑制。本文在其基础上引入第三方回收逆向物流，考虑各环节生产/运输时滞，在随机需求和正弦需求输入作用下，保证良好的动态性能，降低库存波动降低系统中的牛鞭效应。

4 结论

第三方回收商的出现，虽然可以降低生产商的回收物流成本，提升对顾客的服务水平，提高顾客满意度，但对整个闭环供应链库存水平也有一定影响，如果不制订合适的生产/订购策略，可能会导致供应链系统牛鞭效应过大。本文考虑第三方逆向物流的双渠道时滞闭环供应链的现实场景，根据物料平衡方程，建立了供应链增广状态空间模型，推导了系统的预测模型，设计了模型预测控制器，建立SIMULINK仿真模型，通过仿真分析证明了该方法能够在随机需求、正弦信号需求下保证系统稳定，有效抑制生产/运输时滞和客户需求带来的不确定性扰动，得到良好的动态性能，使牛鞭效应得到抑制。