孙军艳, 李晓朋, 陈智瑞, 陈泽飞
(陕西科技大学 机电工程学院, 陕西 西安 710021)
当今世界正经历新一轮大发展、大变革,不确定的市场竞争环境、突发事件的影响以及碳政策的出台导致供应链面临的风险加剧,影响着供应链网络的鲁棒性能.
不确定性是供应链网络的一个固有特征,Zhen L等[1]认为不确定性是影响供应链网络发展的重要因素.目前,学者大多基于两种类型的不确定性因素研究供应链网络的鲁棒优化,一是供应链网络运营过程中参数的不确定性,例如赵潇等[2]考虑了物流网络配送规划问题下需求的不确定性;Peidro D等[3]考虑了需求、成本和供应的不确定性;Attia A M[4]考虑了石油市场需求及价格的不确定性.二是特殊危险事件的不确定性,即供应链网络中的中断风险[5],例如Abbasi S等[6]考虑了物流设施中断的情况;Chen J等[7]考虑了新冠疫情背景下的一个三级供应链中供应商发生大规模中断的情况;王长琼等[8]考虑了制造商、分销商同时发生中断的情况.也有学者同时考虑了运营过程中参数的不确定性和设施中断风险,例如王振等[9]研究了再制造闭环物流网络中需求和回收品质量不确定以及供应中断的问题;Yin M等[10]针对不确定环境下第四方物流网络设计问题,采用情景法来描述不确定需求与中断情景;狄卫民等[11]针对三级供应链网络优化问题,考虑了制造商和分销商两级设施中断与零售商需求不确定的问题.电商时代下,由于消费者对商品型号的期望与商品本身差异、商家对商品的描述与商品本身的差异、物流过程中的商品破损与运输时效性等问题导致退货率居高不下.然而在当前的鲁棒优化模型中,不确定性因素较少考虑产品的退换货问题.
在模型构建方面,有学者以供应链的设施选址为目标,构建了网络选址模型,例如Yu H等[12]针对一个三级闭环供应链网络,建立以最小总成本和碳足迹的多目标设施选址模型;有学者以配送中心库存量及设施节点转运量的分配为目标,构建了流量分配模型,例如Pudasaini P[13]针对需求不确定的供应链网络优化问题,以运输成本及产品损失成本最小为目标建立不同需求情景下的流量分配模型;也有学者以供应链的路径优化为目标,构建了路径规划模型,例如Validi S等[14]针对制造商和分销商组成的两级供应链网络,建立最小化运营成本和碳排放的路径规划模型,有效降低了网络二氧化碳排放量.然而,选址、流量分配和路径优化分阶段的优化很难达到系统的整体最优,有学者将选址和流量分配进行联合优化,例如李进[15]针对闭环供应链网络需求不确定的问题,建立了总成本和碳排放量最小的战略选址及资源流量配置的鲁棒优化模型;有学者将选址和路径规划进行联合优化,例如Zhang Y等[16]针对供应链设施中断概率相同时,用若干个发生概率较大的中断情景集合代替所有中断情景,建立一个设施选址-路径优化模型.现有的研究较少对选址、流量分配和路径优化三者进行联合优化,例如Azizi V等[17]研究了多周期逆向物流网络需求及利润不确定性的问题,建立了两阶段随机规划模型,得到了不确定逆向物流网络最佳设施选址方案与最佳订货批量、运输路线.
电商背景下,针对商品退货率居高不下的问题,本文以市场需求和退换货产品数量为不确定因素,同时考虑设施中断风险,建立了供应链网络总成本和碳排放量最小的选址-流量分配-路径优化的双目标鲁棒优化模型,针对NSGAII算法在求解收敛速度以及保持种群多样性等方面存在缺陷,提出一种Prim-NSGAⅡ算法求解模型,本文的研究为“双碳”时代供应链网络全局优化提出新的求解思路与方法.
如图1所示的闭环供应链网络由制造中心、配送中心和客户需求点组成,正向物流中,产品由制造中心运输至各个配送中心,再由配送中心为客户点提供配送服务;逆向物流中,客户点的退换货产品先收集至配送中心进行检查处理,再将需要返回制造中心的部分产品运至制造中心进行维修再加工.不确定性因素是客户点需求量和逆向退货量,同时考虑制造中心的设施中断风险.在碳中和背景下,政府制定相关碳限额政策,企业通过绿色升级改造来减少碳排放量,决策者需要对不确定闭环供应链网络的设施选址及设施关系分配、产品流量分配、运输路径及减排技术投资方案做出决策.
图1 闭环供应链网络结构示意图
(1)客户需求为单一品种产品.
(2)候选制造中心、配送中心、客户位置已知.
(3)每个客户点仅由一个配送中心为其服务,每个客户点的退换货产品仍由其服务的配送路线进行回收.
(4)每个制造中心只选择一种减排技术投资等级.
(5)客户点的退换货产品中部分有质量缺陷,配送中心进行检查后返回至制造中心维修,退换货产品的返修率已知,维修后的产品再次进入销售网络与新产品质量无差异.
符号说明如表1所示.
表1 集合、参数定义及决策变量符号说明
1.4.1 经济目标
网络总成本最小函数Z1具体表示如下:
(1)设施建设固定成本
(1)
式(1)中:两项分别表示制造中心和配送中心设施建设固定成本;
(2)制造中心生产成本及退换货产品维修成本
(2)
(3)配送中心回收处理成本
(3)
(4)制造中心减排技术投入成本
(4)
为了加快实现碳中和目标,将对制造中心现有技术进行绿色升级改造,r代表对制造中心i进行等级为r的减排技术投入,1/2φr2是一个关于减排技术投入水平的递增凸函数,即节点设施的减排技术投入水平越高,减排技术投入成本越大,产品生产过程中产生的碳排放量越小.
(5)运输成本
(5)
式(5)中:第一项为从制造中心i到配送中心j正向物流运输成本,第二项为配送中心j到需求点k的正向物流运输成本,第三项表示退换货产品由需求点k到配送中心j的逆向物流运输成本,第四项表示需要返厂维修的退换货产品从配送中心j到制造中心i的逆向物流运输成本.
综上,
Z1=min(FC+MC+VC+GC+TC)
(6)
1.4.2 低碳目标
碳排放量最小函数Z2具体表示如下:
(1)设施建设固定碳排放量
(7)
式(7)中:两项分别表示制造中心和配送中心建设的固定碳排放量;
(2)生产过程中产生的碳排放量
(8)
(3)回收检查过程总产生的碳排放量
(9)
(4)运输过程中产生的碳排放
(10)
式(10)中:第一项表示产品从配送中心i至配送中心j运输过程中产生的碳排放量;第二项表示从配送中心j至需求点k运输过程中产生的碳排放量;第三项表示客户退换货产品从需求点k至配送中心j运输过程中产生的碳排放量;第四项表示需要返厂维修的退换货产品从配送中心j至制造中心i运输过程中产生的碳排放量.
综上,
Z2=min(FCE+MCE+VCE+TCE)
(11)
1.4.3 约束条件
闭环供应链网络总成本和碳排放量最小的鲁棒优化模型的约束条件如下公式所示:
(12)
Yjk≤Yj,∀j∈J
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
FCE+MCE+VCE+TCE≤Elim
(25)
(26)
(27)
式(12)~(14)中:代表模型逻辑关系约束,式(12)代表只有选择建设的制造中心才向配送中心提供服务,式(13)代表只有选择建设的配送中心才向客户点提供服务,式(14)为减排技术投入约束,表示每个制造中心只有一种等级的减排技术投入.
式(15)~(18)中:代表网络物流平衡约束,式(15)表示制造中心生产产品全部运输至配送中心,式(16)表示配送中心产品运输总量满足总客户点正向物流需求,式(17)表示回收产品总量必须满足客户点逆向物流需求,式(18)表示符合再制造要求的产品全部运往制造中心进行产品再制造生产.
式(19)~(21)中:代表设施能力约束,式(19)代表制造中心供应生产产品数量不超过其产品制造能力上限,式(20)代表流入配送中心的产品数量不超过其最大容量,式(21)代表配送中心回收产品数量不超过其最大回收能力.
式(22)~(24)中:代表车辆路径规划约束,式(22)~(23)代表每个需求点仅由一个配送中心为其服务,式(24)表示每个配送中心仅服务分配给该设施的客户.
式(25)中:表示闭环供应链总碳排放量不超过政府规定的限制额度.式(26)~(27)中:代表决策变量约束.
1.4.4 不确定需求处理
对于不确定的客户点需求量、不确定的退货量,采用Box不确定集描述不确定参数,基于Bertsimas D等[18]提出的鲁棒对等优化方法得到含有不确定参数的约束式(16)和(17)的鲁棒对应式为:
(28)
s.t.
(29)
pk≥0,zk≥0,∀k∈K
(30)
(31)
s.t.
(32)
rk≥0,tk≥0,∀k∈K
(33)
1.4.5 设施状态描述方法
王知津等[19]通过构建不同的情景集反应系统的不确定性,本文采用情景法描述制造中心的设施中断情景,制造中心生产状态采用参数αi表示.
(34)
假设每个制造中心发生设施中断的概率相互独立,制造中心i中断概率为pi,则制造中心正常工作的概率为1-pi,假设共有S种中断情景,则情景s(s∈S)发生的概率ps为:
(35)
式(35)中:m为发生中断的制造中心.
综上,本文提出的碳中和背景下的不确定闭环供应链网络鲁棒模型如下:
min(Z1,Z2)
(36)
尽管NSGAII算法在求解多目标供应链网络优化问题时应用广泛,但由于交叉、变异操作以及精英策略等自身的局限性,在求解收敛速度以及保持种群多样性等方面存在缺陷.所以本文提出一种Prim-NSGAII算法,通过Prim算法改善初始种群质量,从而达到加快求解收敛速度保持种群多样性的目的,算法流程如图2所示.
图2 Prim-NSGAII算法流程图
具体设计流程如下:
Step1染色体编码
本文所提模型为混合整数模型,在编码时采用0-1编码和实数编码两种方式.每一个染色体代表一种决策方案.每条染色体包含三个染色体组,染色体组1由三个基因段组成,基因段1为0-1编码,基因位代表设施编号,基因值代表设施是否被选择;基因段2采用实数编码,基因位代表设施编号,表示第I个制造中心采用第r种减排技术;基因段3采用实数编码,基因位代表设施编号,基因值代表设施i采用第l种运输方式.染色体组2采用实数编码,表示设施(I,J,K)之间的物流量.染色体组3为运输路径编码段,采用实数编码,对配送中心和需求点依次进行编码,生成规模为N的初始种群C0,染色体具体表示形式如图3所示.
图3 Prim-NSGAII算法编码过程示意图
Step2初始化种群改善
(1)对C0中每一个个体对应的配送中心、客户点形成的连通图,NG=(V1,e1),V1=J∪K,e1={(j,k)|j,k∈V1,j≠k},从J中随机选择一个点j作为顶点,放入生成树节点集U1,选择边权值最小的节点k加入到最小生成树节点集U1,从V1删去j,k.循环该过程,直至所有的客户点k被分配给配送中心.
(3)记由制造中心、配送中心、客户点形成的连通图NT=(V2,e2),V2=I∪UN,e2={(i,j)|i,j∈V2,i≠j},从UN随机选择一个生成树放入节点集合S1,选择与这个顶点边权值最小的节点i加入S1,从V2中删去i,UN,直至所有的制造中心i分配给UN.
Step3选择、自适应交叉、变异
对Ci中所有个体执行快速非支配排序操作,初始化个体rank值,令gen=1,按照如下所示的交叉概率pc和变异概率pm进行交叉变异操作,得到新的子代种群Cg.
(37)
(38)
Step4非支配排序
合并种群Ci与子代种群Cg生成新1代父代种群R1,记R1中每个个体i在该种群中支配的个体数量和支配个体所在的集合为nri,Sri,对于每一个i∈R1,j∈R1,若nri=0,记F1=F1∪{i};rank(ri)=1,若i支配j,则令Sri=Sri∪{j},对任意j∈Sri,令nj=nj-1,若nj-1=0,记F2=F2∪{j},rank(ri)=2,重复上述步骤,直至R1中所有个体被分层.
Step5拥挤度计算
(39)
Step6精英保留
将新种群Ri+1与父代种群R1合并组成种群大小为2N的新种群Gi,根据Pareto等级从小到大依次放入新父代种群Cn中,直到Gi的某一层不能全部放入新父代种群中,将该层个体按照拥挤度从小到大顺序依次放入到新父代种群中.
Step7终止条件
判断是否符合停止条件gen≥Gmax,若是,输出最优解集,否则,令gen=gen+1重复Step4.
某电商自营家电品牌根据品牌发展战略,欲拓展新的消费市场,该品牌企业筛选出6个候选的制造中心,12个候选配送中心,20个客户需求点,各候选设施点位置及客户需求点的位置已知,不同制造中心的建设成本、生产能力、生产成本、与配送中心的单位运输成本及中断概率是不同的,不同配送中心的固定建设成本、最大容量、退换货产品回收处理成本、与客户点之间的单位运输成本是不同的,以往期的销售数据作为不确定参数的名义值.
基于以上案例背景,考虑通过对现有技术进行绿色升级改造来减少碳排放量,考虑3种不同等级的减排技术,其成本系数φ为25 000,投入等级越高生产过程中产生的碳排放量越小,air分别为0.2、0.3、0.6.算例的基本参数设置如表2所示.
表2 算例的基本参数
此外企业自有2种类型车辆,分别为低排放车型与高排放车型,低排放车型单位运输成本较小,但装载能力与碳排放量较大;高排放车型单位运输成本较大,但装载能力和二氧化碳排放量较小,两种车型相关参数如表3所示.
表3 运输相关参数
图4 目标函Z1收敛曲线
图5 目标函数Z2收敛曲线
由图4和图5可知,随着迭代次数的增加,两种算法得到目标函数Z1、Z2均逐步趋于稳定,相比于传统NSGAⅡ算法,Prim-NSGAⅡ算法得到的初始解更优,搜索效率更高,收敛速度较快,在40代左右目标函数收敛,解的质量较传统NSGAⅡ算法解的质量分别提升0.60%、0.79%.
基于以上描述,进一步采用NSGAⅡ算法和Prim-NSGAⅡ算法分别求解本文所建的多目标鲁棒优化模型,得到两种算法下的Pareto前沿分布如图6所示.
图6 Pareto前沿
由图6可知,Prim-NSGAⅡ算法得到的Pareto最优解集总体优于NSGAⅡ算法,Prim-NSGAⅡ算法得到的Pareto前沿相较NSGAⅡ算法得到的Pareto前沿更加靠近Pareto最优前端,而且Pareto前沿分布更为均匀.
反向世代距离IGD通过计算近似最优解集与帕累托最优解集之间的最小欧氏距离,评价近似最优解集的多样性和收敛性,IGD指标计算公式如式(40)所示:
(40)
式(40)中:P和P*代表近似最优解和帕累托最优解集,dis(x,y)表示近似解集P中的点x和Pareto前沿面上的点y的最小欧式距离,IGD值越小,表示近似最优解集的多样性和收敛性越小.
采用空间评价指标SM用来衡量Pareto解集分布的均匀性,SM指标计算公式如式(41)所示:
(41)
在相同软硬件环境下将Prim-NSGAⅡ,传统NSGAⅡ程序运行40次,指标结果如表4所示.由表4可知,本文所提算法的IGD指标均值为0.049 2,最优值为0.041 8,明显优于NSGAⅡ算法;SM指标均值为0.525 3,最优值为0.453 7,略优于NSGAⅡ算法,但是Prim-NSGAⅡ算法的SM指标的标准差及方差均明显小于NSGAⅡ算法.由此可知,在Prim-NSGAⅡ算法中通过改进交叉变异算子,避免了个体由于不完全变化导致算法陷入局部最优,使得Prim-NSGAⅡ算法更为稳定,在收敛性及多样性的改进上具有更好的效果.
表4 IGD及SM指标统计表
通过构建满意度模糊函数在多组可行解中找到最优网络规划方案,目标函数Z1和Z2的隶属度函数如式(42)和式(43)所示:
(42)
(43)
式(42)、(43)中:Zi表示目标函数值,Zimin表示第i个目标函数的解集中的最小值,Zimax表示第i个目标函数的解集中的最大值.定义双目标函数最大满意度λ为:
(44)
由图7和图9可知,确定性模型时选择制造中心数量为5,编号分别为1,2,3,5,6,选择配送中心数量为5,编号分别为1,4,5,6,9.例如,由制造中心1向配送中心1供应6 235件产品,由配送中心1为客户点4,16,11,13依次提供服务;并回收退换货产品612件,再由配送中心1将180件需要返修的产品运输至制造中心1.
图时网络规划方案求解结果
表6 单目标鲁棒优化模型的设施选址结果
网络总成本及二氧化碳排放节约量随不确定程度的变化趋势如图11和图12所示.
图11 不确定程度对网络总成本的影响
图12 不确定程度对CO2节约量的影响
由图11可知,随着不确定程度γ的增加,网络总成本大致线性增加.此时除了增加设施选址、生产、运输成本来维持网络鲁棒性,还需要增加减排技术投资成本,将网络中产生的二氧化碳维持在碳排放限制额度之内.
由图12可知,随着不确定程度γ的增加,CO2排放节约量增加,当γ=0.2时,二氧化碳排放节约量增幅明显,由表5可知,此时,供应链网络选择对制造中心4进行等级为1的减排技术投资方案,证明碳限额政策可以有效促使决策者做出更加低碳的决策.
制定合理的碳排放限制额度可以将企业二氧化碳排放量控制在一定水平,最大化减排效率的同时尽可能减少企业经济负担.为探究碳排放限制额度对不确定闭环供应链网络优化结果的影响,令碳排放限制额度Elim分别等于33 000、34 000、35 000、36 000、37 000、38 000,网络总成本与二氧化碳排放量随不确定程度变化趋势如图13和图14所示.
图13 碳限额对网络总成本的影响
图14 碳限额对网络CO2排放量的影响
由图13可知,碳限额的变化对不同不确定程度的网络总成本均有一定影响.当网络不确定程度γ≤0.6时,网络总成本随碳排放限额Elim的增加小幅减少;当γ>0.6时,网络总成本随Elim的增加减少明显,说明Elim对不确定程度较小的供应链成本影响不大.
根据图14可知,当网络不确定程度γ<0.8时,CO2排放量随碳排放限额Elim呈阶梯型递增的趋势;当γ≥0.8时,CO2排放量随Elim呈线性递增趋势.
当不确定程度较小,碳限额政策较宽松时,供应链总成本较低,碳排放量较高.此时,大幅度减少碳排放限制额度,可使得二氧化碳排放量降低,同时总成本增幅较小.当不确定程度较高,碳限额政策较严格时,供应链总成本较高,碳排放量较少,此时,只需要小幅度减少碳限额,就可使得二氧化碳排放量降低,但是需要付出的经济成本较大.
综上,不管不确定程度如何变化,碳限额政策对网络CO2排放将会起到较好的约束效果.随着Elim的增加,即碳政策的监管力度变宽松时,企业会改变减排技术投资方案,牺牲环境使总成本更小;当Elim大于某一额度时,总成本不再发生变化,碳限额政策不再起作用.
针对客户点需求量、退货量的不确定因素,考虑制造中心设施中断风险,本文建立了以最小化闭环供应链网络经济成本与碳排放的双目标鲁棒优化模型,采用Box和情景法处理不确定因素,设计了一种Prim-NSGAⅡ算法求解双目标鲁棒优化问题,采用Prim算法提高初始种群质量,并采用Prim-NSGAⅡ算法求解得到最优的网络规划方案,经验证可得以下结论:
(1)通过与传统NSGAⅡ算法的IGD、SM指标对比验证了算法的有效性,该算法达到了加快求解收敛速度、保持种群多样性的目的.
(2)不管网络不确定程度如何,碳限额政策对二氧化碳排放均会起到较好约束效果,当不确定程度较高时,碳限额政策将会提高经济成本.