基于Kriging元模型的多级供应链稳健优化设计

朱连燕,吴 锋,欧阳林寒

(1.南京科技职业学院 基础科学部,江苏 南京 210048;2.安徽工程大学 经济管理学院,安徽 芜湖 241005;3.南京航空航天大学 经济管理学院,江苏 南京 210016)

0 引言

随着工业4.0的推进及经济全球化发展,产品供应链涉及的企业数量愈加庞大,动态性强、复杂化等特点使得多级供应链研究受到越来越多学者的关注[1-2]。相较于二级供应链,多级供应链可以比较直观地反映商品的运作流程,更具实用性和代表性,对供应链管理更有现实意义[3]。从流程的角度,多级供应链可看作是一个动态的多阶段复杂过程系统,仿真建模技术与优化方法相结合已逐渐成为分析、优化这样复杂系统的有效手段之一[4]。然而,由于多级供应链系统自身的复杂性和动态性,仿真模型的运行往往十分耗时[5]。仿真模型只能给出某种输入的输出值,输入与输出之间的解析函数关系往往未知,仿真模型的这种黑箱性使得常见的确定性优化方法失效[6];而且仿真模型的昂贵运行成本降低了优化方法的效率,增加了寻优的复杂性[7]。为了提高仿真试验效率,降低其运行成本,基于仿真试验和元模型的供应链优化方法应运而生[8]。

目前多项式元模型以其结构简单而成为众多研究者研究供应链的首选[9-11]。SHI等[12]针对汽车零部件供应链的多响应优化问题,提出一种基于离散事件仿真建模、序贯分支法和响应面方法的优化方法,以最大限度地提高物流绩效。SHABAN等[13]考虑到供应链中订单和库存的变化,提出一种将仿真建模和多项式响应曲面相结合的混合方法,试图寻找当需求和库存方差之和达到最小时指数平滑策略参数的最优值。TANG等[14]将arena仿真和多项式元模型相结合,基于稳健参数设计思想,有效解决具有随机需求和随机提前期的供应链库存策略优化问题。YANG等[15]在文献[14]研究的基础上,将均值条件风险值与响应曲面法相结合,解决了具有风险偏好的多响应供应链库存策略优化问题。然而,多项式元模型受预设结构的限制对于强非线性问题难以处理,并且对多峰函数问题近似精度较差[16]。因此,对于多级供应链这类强非线性的仿真优化问题并不适合。

随着研究的不断深入,Kriging元模型逐渐引起研究者的关注[17-19]。该模型对于强非线性输入输出关系有很好的拟合和预测能力,同时能够提供这种近似不确定性的度量,被逐渐应用于解决供应链相关问题的研究[20-21]。针对具有多源能力、异步订货、不确定需求和随机提前期的一般供应链库存控制问题,YE等[22]通过Kriging元模型对于供应链中的每个库存系统性能进行区域估计,在此基础上,提出一种集成的Kriging元模型的优化模型,该方法在满足服务水平约束的条件下,能够有效降低运营成本,大大减少了采样和优化过程中的计算时间。考虑到决策制定者的风险态度,朱连燕等[23]结合Kriging元模型和条件风险值准则,提出了基于Kriging模型的均值—条件风险值优化策略,有效解决了具有风险规避特性的库存管理问题。针对供应链优化问题,朱连燕等[24]提出了基于Kriging元模型的供应链优化方法,有效解决了三级供应链的整体优化问题,但未考虑到不确定性因素对供应链整体绩效的影响。随着不确定因素的增加,MEDINA-GONZALEZ等[25]结合多参数规划法,提出一种基于Kriging元模型的数据驱动决策框架理论,解决了多目标生物能源供应链网络的优化管理问题,有效降低了求解最优值的复杂度。

通过上述文献分析发现,基于Kriging元模型的优化方法研究供应链相关问题,涉及多个绩效响应指标的多级供应链整体优化的相关研究并不多见。事实上,多级供应链涉及的相关方较多,需求复杂,不确定性因素对整体绩效的作用不可忽视,研究多个绩效响应指标的供应链整体稳健优化更加符合现实,在优化过程中如何平衡多个绩效响应指标则尤为重要。因此,本文以四级供应链为例,首先基于仿真模型的试验数据分别构建供应链绩效响应——整体利润的均值响应和标准差响应及整体库存的均值响应和标准差这四个Kriging元模型;其次,考虑到多个绩效响应之间的平衡,在所构建的Kriging元模型的基础上,利用满意度函数法构建综合绩效的稳健优化策略;最后利用遗传优化算法进行求解,确定在综合绩效满足一定约束条件时,综合绩效波动最小对应的最优参数值。

1 多级供应链arena仿真模型

现实生活中的产品供应链多数为多级供应链[26]。啤酒游戏是一种典型的多级供应链系统,本文以啤酒游戏的供应链系统为例,假设一个啤酒游戏供应链系统由1个制造商、2个分销商、4个批发商、8个零售商组成,具体结构如图1所示。仿真假设零售商反馈时间为0,批发商(或分销商)反馈时间为4天,制造商的反馈时间为4天,制造商生产周期为1天,供应链中各成员均采用定期订货策略。以供应链的总体利润和总体库存为优化目标,总体利润由销售收入、订货成本、库存成本、销售损失及惩罚4个部分构成,总体库存则为各级企业的总计现有库存。为了利用arena仿真软件对所描述的供应链系统进行建模,首先要明确其输入输出情况,从系统的角度给出该供应链系统的示意图如图2所示。

图1 多级供应链的结构示意图

图2 多级供应链的系统模型

在供应链的实际运营中,顾客需求的波动导致生产水平和库存水平的波动,同时整个供应链的订货时间、运输时间和订货提前期具有一定的随机性。仿真模型主要分为“需求产生”,“下订单”,“收到订单”,“延迟订单”,“订货”,“发货”,“收货”,“完成订单”等8大功能模块,总体仿真流程逻辑模型如图3所示。假设在仿真运行规定的时间内需求水平和价格结构固定,将订货点和订货数量看作可控因子(或决策变量),顾客的需求量和供应的提前期为噪声因子(或环境因子),arena仿真模型记录供应链中各节点企业的库存水平和订单数量,并动态展示了供应链中的订单传递和完成情况及库存变动情况。仿真模型为终态(Terminating)仿真模型,考察供应链两年的实际运营情况,仿真模型运行720个仿真日(每日工作24小时),系统重复运行30次,另有10个仿真日的系统Warm-up时间,系统运行时制造商、分销商、批发商和零售商的初始库存水平分别设置为300、500、300和150,它们的订货点分别设置为200、400、200和100,订货数量分别设置为700、300、150和100。系统的每日需求量由固定和波动两部分组成,其中固定值为25,波动部分服从参数为25的泊松分布,仿真系统的具体参数设置、提前期的参数设置、供应链的各级价格参数设置及具体参数运行结果参考文献[27]。

图3 仿真流程逻辑总图

2 基于Kriging元模型的多绩效响应供应链系统的稳健优化策略

2.1 Kriging元模型

Kriging元模型是一种基于随机过程的广义回归模型,它不仅对于非线性输入输出关系具有很好的拟合效果,还能够给出这种拟合误差的度量[28]。Kriging元模型将黑箱函数y(x)看作包含回归项的随机过程的Y(x)的一次实现。对于给定的样本点x,对应的函数响应值为y(x),则y(x)为随机函数Y(x)的可能取值之一,其一般表达式可记为:

(1)

Cov(Z(xi),Z(xj))=σ2R(xi,xj,θ)。

(2)

其中:σ2为随机过程Z(x)的方差;R(xi,xj,θ)表示给定的相关函数,一般选取高斯函数、指数函数等作为相关函数;θ表示相关函数的参数,通过优化θ,能够自适应调节样本点xi,xj之间的空间相关性,本文选取高斯函数作为相关函数,其表达式为:

(3)

可以通过最大似然估计法得到参数σ2、θ、βj的估计值。游海龙等[29]指出,Kriging元模型是一种基于黑箱理论的插值模型,不具有明确的显函数形式,经常作为复杂系统仿真模型的近似模型,利用传统梯度类优化方法优化Kriging元模型存在困难,而遗传算法和Kriging元模型相结合则能够很好地解决这类优化问题。

2.2 基于Kriging元模型的稳健优化策略

事实上,根据第1章构建的供应链仿真模型,供应链这类优化问题一般可表示为:

s.t.

L1≤x1≤L2,

L3≤x2≤L4,

L5≤x3≤L6,

L7≤x4≤L8,

L9≤x5≤L10,

L11≤x6≤L12,

L13≤x7≤L14,

L15≤x8≤L16。

(4)

其中:式(4)中的目标函数为综合绩效fl(x),此时fl(x)=Ee[ωl(x,e)],l=1,2(l=1时,计算供应链的整体库存(IN),l=2时,计算供应链的整体利润(TP))。ω1(x,e)表示给定可控因子x、噪声因子e时对应的供应链整体库存的仿真输出值,f1(x)表示噪声因子e服从某种分布时对应供应链整体库存的仿真输出均值;ω2(x,e)表示给定可控因子x、噪声因子e时对应的供应链整体利润的仿真输出值,f2(x)表示噪声因子e服从某种分布时对应供应链整体利润的仿真输出均值。x=(x1,x2,…,x8)表示可控因子(或决策因子),Θ∈R8为可控因子的可行域,满足约束条件中所规定的数值区域。本文中的可控因子为制造商、分销商、批发商和零售商的订货点(reorder point)和订货数量(reorder quantity),其中:x1表示零售商的订货点;x2表示批发商的订货点;x3表示分销商的订货点;x4表示制造商的订货点;x5表示零售商的订货数量;x6表示批发商的订货数量;x7表示分销商的订货数量;x8表示制造商的订货数量。可控因子之间满足如下约束条件:x3≥x2≥x1;x7≥x6≥x5。Li(i=1,2,…,16)为可控因子的取值范围。e表示噪声因子(或环境因子),以供应链模型中市场需求的波动和订货提前期的波动为噪声因子。YANG等[15]指出,批发商的订货点、制造商的订货点、零售商的订货量及制造商的订货量对于供应链的绩效响应库存和利润影响最大,因此可将模型(4)中的8个可控因子减少至以下4个关键可控因子,模型(4)可转化为:

s.t.

L3≤d1≤L4,

L5≤d2≤L6,

L9≤d3≤L10,

L15≤d4≤L16。

(5)

其中:d=(d1,d2,d3,d4)T表示上述4个关键因子所构成的向量,E1(d)=f1(d)表示整体库存(IN)的仿真输出均值,E2(d)=f2(d)表示整体利润(TP)的仿真输出均值,D∈R4表示可行域,是约束条件中所规定的数值区域。L3,L4表示可控因子d1的取值范围;L5,L6表示可控因子d2的取值范围;L9,L10表示可控因子d3的取值范围;L15,L16表示可控因子d4的取值范围。

由于噪声因子对供应链系统产生的影响不可忽视,在供应链优化过程中不仅需要关注供应链绩效的最优解,还需要度量最优解的稳健性。鉴于此,考虑到随机需求和提前期的不确定性,基于稳健参数设计的思想,构建多级供应链的稳健优化策略来优化供应链系统中的可控因子水平组合以便降低噪声因子对供应链系统综合绩效的影响则是本文需要解决的问题。此时,优化问题(5)可转化为:

s.t.

L3≤d1≤L4,

L5≤d2≤L6,

L9≤d3≤L10,

L15≤d4≤L16,

E1(d)≤H1,

E2(d)≥H2。

(6)

其中:S1(d)和S2(d)分别表示供应链系统绩效库存和利润的标准差响应;目标函数是S1(d)和S2(d)复合绩效的最大化;E1(d),E2(d)分别表示供应链的绩效库存(IN)和利润(TP)均值响应模型;H1,H2分别表示库存和利润的均值取值范围。考虑到供应链仿真模型的绩效响应具有多重性,即库存(IN)和利润(TP),由于库存是望小的质量特性,利润是望大的质量特性,为了在稳健优化过程中平衡这两个响应之间的冲突性,采用满意度函数法[30]构建多响应稳健的优化目标函数。即先分别建立S1(d)和S2(d)各自的满意度函数g(y1)和g(y2):

(7)

其中:U1为库存这一绩效响应的规格上限,T1为其目标值。当响应值小于目标值T1时,则其满意度g(y1)=1,当响应值大于规格上限U1时。其满意度g(y1)=0;L2为利润的规格下限,T2为其目标值,当响应值超出其目标值T2时,其满意度函数为g(y2)=1;当响应值小于规格下限L2时,其满意度函数为g(y2)=0。然后采用几何平均的思想定义综合满意度函数,从而模型(6)的优化目标函数可表示为:

(8)

其中:y1=S1(d),y2=S2(d)分别表示供应链整体库存和整体利润的标准差模型,为了提高仿真效率和优化成本,首先利用少量的仿真输入输出数据,分别构建S1(d),S2(d)的近似模型——Kriging元模型,类似地,基于仿真试验数据分别构建库存和利润的均值E1(d),E2(d)的Kriging元模型,然后在此基础上,构建稳健优化策略模型。

3 仿真结果

3.1 库存和利润的Kriging元模型建立

借鉴稳健参数设计的思想和arena仿真的试验数据,分别拟合两个绩效响应(库存IN和利润TP)的均值响应和标准差响应的Kriging元模型,共构建4个Kriging元模型。其具体实施步骤如下:

(9)

(10)

步骤2元模型的构建。利用步骤1获得的试验数据,分别构建均值响应模型E(d)和标准差响应模型S(d)的Kriging元模型。

步骤3元模型的验证。利用“留一”交叉验证法[21]验证元模型是否有效,若有效,则转步骤4,否则转步骤1,直至所构建的元模型有效。

步骤4稳健优化策略的构建。基于满意度函数法,利用步骤2所构建的元模型,构建基于Kriging元模型的均值-标准差稳健优化策略模型,此时模型(6)可转化为模型(11):

s.t.

L3≤d1≤L4,

L5≤d2≤L6,

L9≤d3≤L10,

L15≤d4≤L16,

E1(d)≤H1,

E2(d)≥H2。

(11)

步骤5采用bootstrap抽样方法给出稳健优化解的置信区间以度量噪声因子的波动对于优化解的影响。

上述步骤如图4所示。

图4 基于Kriging元模型的多级供应链的稳健优化流程图

3.2 仿真试验及分析

表1 Kriging元模型的部分参数

图5 基于Kriging元模型的库存和利润绩效均值响应和标准差响应散点图

针对本文的仿真案例,结合优化算法优化模型(11),可将模型(11)转化为:

s.t.

150≤d1≤300,

300≤d2≤600,

75≤d3≤150,

600≤d4≤900,

E1(d)≤2000,

E2(d)≥200000。

(12)

g(S1(d))=

(13)

g(S2(d))=

(14)

其中:123.97和645.89为库存(IN)的300个试验点标准差的最小值和最大值;27 073.15和40 000为利润(TP)的300个试验点标准差的最小值和最大值。Kriging元模型是一种基于黑箱理论的插值模型,不具有明确的显函数形式,利用遗传算法优化Kriging元模型能够得到理想的优化结果[29]。因此本文采用遗传优化算法求解模型(12),得到最优可控因子水平设置为(0.8,-1,-0.35,0.68),综合满意度为0.821,最优可控因子水平及最优响应值(IN和TP)的结果如表2和表3所示。

表2 基于Kriging元模型的可控因子最优水平设置

表3 基于Kriging元模型最优水平设置对应的响应值

与基于Kriging元模型的优化方法类似,给出利用基于多项式元模型的优化方法求解模型(12),利用仿真试验数据构建的4个多项式元模型,其预测精度分别为2.2%、1.38%、2.4%、2.6%,对比表1的最后一列,表明Kriging元模型具有更高的预测精度。表4和表5给出了基于多项式元模型得到的最优可控因子水平及最优响应值。此时,达到最优解对应的综合满意度为0.469,显然低于基于Kriging元模型的综合满意度0.821。通过对比表2和表4可知,基于Kriging的最优因子组合绩效预测值略优于多项式最优绩效预测值,通过表3和表5数值试验结果对比,显示了Kriging元模型的稳健解优于多项式元模型。

表4 基于多项式模型的可控因子最优水平设置

表5 基于多项式模型的最优水平设置对应的响应值

一般情形下,最优结果不一定最稳定,故需要对其进行稳健性分析。根据需求量和提前期的分布,用拉丁超立方抽样抽取1 000个样本点,对表2和表4中可控因子分别取步长±1、±0.1、±0.01,在各自的可控因子水平区间内等距离取1 000个点,采用上述优化策略,可以得到3种步长下的1 000组最优解的均值和标准差,结果如表6所示,由表6可以看出,基于Kriging元模型的优化解的均值和仿真模型得到的满意度均值最接近,而且标准差最低,表明该优化策略稳健性能最优。

表6 不同优化策略的稳健性能比较

4 基于bootstrap的稳健优化

针对利润这一绩效响应,如图6所示为基于多项式元模型和基于Kriging元模型的1 000个bootstrap稳健最优解的样本估计值和置信区间,类似地,如图7所示为针对库存这一绩效响应基于不同元模型(多项式和Kriging)的1 000个bootstrap稳健最优解的结果。符号“+”和“o”分别代表基于多项式元模型和Kriging元模型的最优解的1 000个bootstrap样本估计值,颜色相同的虚线间区域表示该最优解的bootstrap置信区间,其中:红色虚线表示基于Kriging元模型得到的最优解的bootstrap置信区间,绿色虚线表示基于多项式回归模型得到的最优解的bootstrap置信区间。如在图6中,基于多项式元模型和Kriging元模型的利润(TP)最优解的置信区间分别为[198398.3421,209259.8911],[197619.8651,208744.2148],相应地,在图7中,基于多项式元模型和Kriging元模型的库存(IN)最优解的置信区间分别为[1207.976,1281.382],[1024.397,1100.542]。

图6 基于不同元模型的利润bootstrap估计点、置信区间及预测值

图7 基于不同元模型的库存bootstrap估计点、置信区间及预测值

直线表示最优解所对应的元模型的预测值,其中:蓝色直线表示基于Kriging元模型的预测值,黑色直线表示多项式元模型的预测值。从图6中可以看出,元模型的预测值(直线)介于置信区间(虚线)内,说明所构建的元模型具有较高的可信度,所拟合的两类元模型能够较好地反映啤酒供应链系统的绩效情况。虽然图6中基于两类元模型得到的各自bootstrap估计值的分布并没有显著差异,但纵观总体,图6中基于多项式元模型(“+”型符号)的利润(TP)bootstrap估计值绝大多数位于基于Kriging元模型的bootstrap估计值(“o”型符号)上方,水平方向上的分布并没有Kriging元模型的bootstrap估计值集中,而垂直方向上却略高于Kriging元模型的bootstrap估计值,表明多项式虽然受到了稳健优化式约束条件的限制,但并不满足目标函数标准差最小化的要求,导致基于多项式元模型最优解的稳健性劣于基于Kriging元模型最优解的稳健性,说明了基于Kriging元模型得到的最优解能使绩效响应的波动幅度处于较小范围,更加稳健。图7中基于多项式元模型的预测值(黑色直线)落在置信区间外,而基于Kriging元模型的预测值(蓝色直线)却落在置信区间内,说明对于库存(IN)这一绩效响应,相较于多项式元模型,Kriging元模型具有较高的预测精度。这说明Kriging元模型非常适合拟合具有非线性关系的四级供应链系统,而多项式元模型拟合这样的系统并不太合适。

5 结束语

本文考虑多级供应链的非线性、复杂性及不确定性因素所引起的系统波动,结合Kriging元模型的优势,提出了一种基于Kriging元模型的综合满意度稳健优化设计方法。首先针对多级供应链中各节点企业需求及提前期的波动对供应链整体绩效的影响,借鉴稳健参数设计的思想,将arena仿真和Kriging元模型相结合,构建基于Kriging元模型的均值—标准差模型;其次,针对供应链的各绩效响应指标之间的冲突性,利用满意度函数法,构建基于Kriging元模型的综合满意度稳健优化策略,进而确定系统综合绩效最优时稳健最优参数值,并和基于多项式元模型的优化方法进行比较;最后,采用非参数bootstrap方法度量不确定性因素对优化结果的影响,并和基于多项式元模型的稳健优化结果进行对比分析。仿真结果表明:①和多项式元模型相比较,基于Kriging元模型的综合满意度稳健优化方法无论在预测性能和优化结果上都更胜一筹;②基于Kriging元模型所得的优化结果更加稳健;③将满意度函数法和Kriging元建模思想相融合,能够有效平衡多级供应链中多个绩效指标的冲突性,在保证供应链系统整体绩效达到最优的同时,尽可能地降低不确定性因素对供应链绩效的影响,为多级供应链的稳健运营提供了一定的理论依据和决策参考。然而,由于本文所考虑的供应链系统仅是单一产品的多级供应链,而多产品的供应链运营可能面临更为复杂的环境,如何基于Kriging元模型的优化方法解决多产品的供应链稳健优化将是下一步研究的方向。