(中南大学交通运输工程学院,长沙 410075)
在 《车辆运输车治理工作方案》出台后,整车物流领域发生重大变革,自2018年7月1日起,不合规车辆运输车被全面禁止通行,整车公路运输成本急剧上升,企业开始采用图1所示的整车“库前移”物流网络代替原有的公路直达运输网络。该网络涉及配送中心选址、库存及运输等环节,需求量的变化会导致这3个环节的成本变化,且由于整车需求量具有明显的季节性规律,在对整车物流网络进行优化时需考虑不同时期的不同需求量要求,使得决策方案具有较强的鲁棒性。此外,由于环境问题日益严重,而物流配送作为碳排放的主要来源之一,考虑碳税政策对整车物流网络优化方案的影响具有非常重要的现实意义。
图1 整车物流网络
针对低碳视角下的物流网络优化问题,Yang等[1]提出了碳税约束型城市物流配送网络规划模型;Jin等[2]对不同碳减排政策下大型零售商的供应链网络优化及运输方式选择进行了研究;成耀荣和谭维[3]构建了基于低碳经济视角的多式联运路径优化模型;Saxena等[4]研究了碳税和减排奖励政策对印度轮胎再制造物流网络设计的影响。针对整车物流网络优化问题,王沛等[5]研究了多层级、多周期的整车铁路物流基地选址;秦璐等[6]利用正交试验方法对整车物流网络优化进行了研究;郑燕等[7]对长江流域整车公水联运网络进行了优化。在物流网络鲁棒性研究方面,Xi等[8]研究了非等区域动态设施布局,提出了考虑货物装卸点位置的混合鲁棒优化模型;周愉峰等[9]针对国家血液战略储备库选址,构建了考虑多情景、多血型的随机p-鲁棒选址-库存优化模型;孙华丽等[10]考虑应急需求的随机性,建立了多目标选址-路径鲁棒优化模型。
已有文献大多是单独考虑物流网络鲁棒性或碳排放而建立的优化模型,鲜有文献同时考虑这两者的共同影响并针对整车物流行业背景进行研究。因此,本文同时考虑碳税政策和多种需求量情景,建立考虑碳税政策的整车物流网络鲁棒优化模型,设计遗传退火算法进行求解,为整车物流网络优化问题提供更全面的解决思路。
由主机厂、配送中心和4s店构成的整车物流网络,如图1所示。由于不同时期整车的需求量具有较大差异,研究采用 “情景”处理方式,即设置若干个典型的需求情景,每种情景对应不同的需求量。要解决的问题是:在碳税政策下,确定配送中心的数量和选址位置、配送中心与4s店的服务关系、配送中心的库存控制参数以及主机厂至配送中心的运输方式,使得该方案在满足鲁棒性约束的要求下,在所有需求量情景下的总成本期望值最小化。
(1)配送中心容量能够满足需求;
(2)1个主机厂和若干4s店的位置已确定,整车从主机厂至4s店经过且只经过1个配送中心,不考虑从主机厂至4s店的直达运输;
(3)各情景发生的概率可以主观设定;
(4)各个4s店的需求相互独立,且单位时间 (天)内的需求量均服从正态分布;
(5)主机厂至配送中心的干线运输可选择公路、铁路和水路运输中的一种,但从配送中心至4s店只采用公路运输进行配送。
(6)配送中心采用连续检查的(Q,R)库存控制策略,而4s店只保有少量的库存,忽略4s店的库存成本。
I配送中心候选点集合,i∈I;
J4s店集合,j∈J;
O整车主机厂所在地点;
G需求情景集合,g∈G;
K运输方式集合 (k∈K,k=1,2,3分别表示公路运输、铁路运输和水路运输);
ck运输方式k的单位运输成本 (元/(km·辆));
ek运输方式k的碳排放因子(kg/(km·辆));
doi主机厂o与配送中心i的距离;
dij配送中心i与4s店j的距离 (km);
α库存服务水平;
Zα库存服务水平为α时的安全库存系数;
si配送中心i的安全库存量 (辆);
hi配送中心i的储存成本 (元/(辆·天));
Ri配送中心i的订货点 (辆);
Ai配送中心i的平均库存 (辆);
TPi配送中心i连续两次订购的间隔时间(天);
Q配送中心i的订购批量 (辆);
fi配送中心i的固定设施日折旧成本 (元/天);
从主机厂至配送中心i采用运输方式k的提前期 (天);
OC配送中心补货的固定成本 (元/次);
φ碳税 (元/kg);
qg情景g发生的概率;
p遗憾值限定系数。
决策变量:
Xi在节点i建立配送中心时为1,否则为0;Yijj由配送中心i服务时为1,否则为0;Uik从主机厂o至配送中心i采用运输方式k时为1,否则为0。
(1)库存成本分析
整车属于高价值商品,需要经常性的进行库存盘点,因此配送中心采用连续检查的(Q,R)库存控制策略。根据模型假设 (4),由数理统计知识可知配送中心的需求服从正态分布,其中:
根据库存控制相关理论,配送中心的日常库存控制参数如下:
①提前期:
②安全库存:
③订货点:
④平均库存:
⑤连续两次订购的间隔时间:
将式 (8)对Qi进行求导,并令其等于0,可得到配送中心i的最佳订购批量, 如式 (9)所示:
将式 (9)代入式 (7)中即可确定TPi的值。
(2)运输成本分析
运输成本包括主机厂至配送中心的长途干线运输成本、铁路及水路运输两端的接驳换装成本、以及从配送中心至4s店的配送成本,则单位时间内的平均运输成本如式 (10)所示:
(3)碳排放成本分析
物流网络中的总碳排放量包括干线运输、中转以及配送产生的碳排放成本,则单位时间内的平均碳排放成本如式 (11)所示:
(4)固定设施折旧成本分析
配送中心i的固定设施日折旧成本如式 (12)所示:
综上所述,在情景g下系统总成本Zg可表示为各项成本之和,如式 (13)所示:
(1)基于特定情景的确定性模型
对于某种特定情景g,其需求量参数是确定的,综合以上的分析,可建立如下确定性优化模型A:
式 (14)为最小化系统总成本,式 (15)表示每个4s店只能分配给1个配送中心;式 (16)表示干线运输只能选择1种运输方式;式 (17)确保4s店不会分配给未选中的配送中心,式 (18)表示至少有1个备选配送中心被选中,式 (19)~(21) 为0~1变量约束。
(2)鲁棒优化模型
由模型A虽然可求得每个情景所对应的优化方案及目标函数值Z∗g,但由于不同情景下的需求量有较大差异,该情景下的优化方案在其他情景中可能表现会很差。因此,在模型A的基础上建立鲁棒优化模型B。模型B力求得到一个在所有情景下的目标函数值与该情景下的最优目标函数值的差值在可接受范围内的满意解,其定义如下所示[11]: 设为情景g下的最优目标函数值,r为满足所有情景g∈G下的可行解,为r在情景g下的目标函数值,如果解r满足鲁棒约束条件,则称r为该问题的鲁棒解,其中参数p被称为遗憾值限定系数,反映了模型的鲁棒性水平,p值越小,则模型的鲁棒性越强。根据定义可知,满足约束条件的鲁棒解r有多个,其中必有一个最佳的解r∗使得总成本最小。模型B如下所示:
其他约束条件同式 (15)~(21), 式 (22) 为最小化鲁棒解在各情景下成本的期望值,式 (23)表示各情景发生的概率之和为1,式 (24)表示鲁棒约束条件。此外,当不考虑约束条件 (24)时,模型B即等同于随机规划模型C。
选址—库存问题已被证明为NP难问题[12]。而本文所研究的整车物流网络优化涉及选址、库存、运输方式选择等,相比选址-库存问题而言更为复杂,难以用精确算法进行求解。在求解该类问题的启发式算法中,遗传算法的全局寻优能力强,收敛速度快[13],但其局部搜索能力较差,易产生 “早熟”现象[14]。而模拟退火算法具有较强的局部搜索能力,因此,本文结合这两种算法优点,设计遗传退火算法对该问题进行求解,即在遗传操作陷入局部收敛后 (算法所得的历史最优解在一定代数内不发生变化时可视为发生局部收敛),再引入模拟退火操作使算法跳出局部最优解。先给出求解模型A的具体步骤,如下所示:
Step1:编码。采用实数编码,每条染色体的编码分为两段,第1段编码包含n位基因,对应4s店与配送中心的分配关系,取值范围为[1,m]之间的整数;第2段编码包含m位基因,对应主机厂至配送中心的运输方式选择变量,取值范围为[1,3]之间的整数。
Step2:随机生成规模为popsize的初始种群1,设置算法最大迭代次数为maxgen,当前迭代次数为gen,令gen=0。以目标函数作为适应度函数,记录种群1中最优个体的适应度值作为历史最优解。
Step3:选择。结合 “轮盘赌”与 “精英保留”策略对种群进行筛选。
Step4:交叉。采用单点交叉方式,如图2所示。
Step5:变异。采用交换变异方式,由于染色体的两段编码的取值范围不同,为了避免变异后产生非法个体,只有当随机选择的两个变异点位于同一段编码时才交换两点的值。如图2所示:
图2 交叉、变异过程
Step6:记交叉变异后的种群为种群2。用种群1中的最优个体替换种群2的最差个体,产生新种群3,计算适应度,更新历史最优解及其保持不变的代数。
Step7:判断是否满足模拟退火操作执行条件。当历史最优解保持不变的代数超过100代时为模拟退火操作执行的必要条件。若满足该条件,则从种群3中随机选择比例为∂的个体构成模拟退火操作的初始解集,执行Step8,否则,跳转至Step13。
Step8:初始化模拟退火操作。令当前解为w0,初温为T0,温度下降系数为low,每一温度下的最大迭代次数为maxnum,令当前温度T=T0,当前温度下已迭代次数num=0。
Step9:随机扰动生成新解w1。新解产生机制为随机选取w0中染色体π的某个位置γ的基因作为扰动点,当γ位于第1段编码时,扰动范围为[1,m]之间的整数,如式 (25);当γ位于第2段编码时,扰动范围为[1,3]之间的整数,如式(26)。其中ceil为MATLAB中朝正无穷方向取整函数,rand为(0,1)之间的随机数。
Step10:根据Metropolis准则判断是否接受新解。令Δw=w1-w0,当Δw<0,则接受w1作为新的当前解w0, 否则以概率 exp(-Δw/T)接受w1作为新的当前解w0,num=num+1。
Step11:判断当前迭代次数num是否小于当前温度下的最大迭代次数maxnum。当num<max-num时,返回Step9;否则,执行Step12。
Step12:用新解替换种群3中的原解。更新当前温度T,执行降温操作T=low∗T0,l=l+1,并令num=0。
Step13:判断是否满足算法终止条件。当gen=maxgen时,算法终止。否则,gen=gen+1,返回Step3。
以目标函数和罚函数之和作为Step3中的适应度函数,如式 (28)所示:
现设计由1个主机厂、12个备选配送中心及40个4s店构成的物流网络仿真算例,由于篇幅限制仅列出6个备选配送中心及20个4s店的数据,设置备选配送中心相关参数如表1所示,备选配送中心与4s店的距离如表2所示,其他各项参数如表3所示,设置3种典型情景1、2、3分别代表需求量低、中、高的情况,需求均值及方差如表4所示。
表1 备选配送中心相关参数
表2 备选配送中心和4s店的距离 (km)
表3 其他参数取值
表4 不同情景下4s店需求均值及标准差数据
续 表
表5 鲁棒优化模型与随机规划模型性能对比
为了比较鲁棒优化模型B和随机规划模型C的性能,以上述数据为基础设置3组不同节点规模的物流网络,分别用模型B和模型C进行求解。令情景g下的最优目标函数值为,鲁棒优化解在情景g下的目标函数值为,鲁棒优化解的总成本为Zr,的相对遗憾值,ωr的方差为Vr。 随机规化解在情景g下的目标函数值为,随机规划解的总成本为Zs,的相对遗憾值,ωs的方差为Vs。Zr与Zs的相对差值为Zr-s=(Zr-Zs)/Zs×100%,Vs与Vr的相对差值为Vs-r=(Vs-Vr)/Vr×100%。 取p=0.12, 结果如表5所示。
由表5可知,鲁棒优化解的相对遗憾值方差远小于随机规划解的相对遗憾值方差,两者的相对差值最大达到了533.33%,而方差越小说明遗憾值的波动越小,在各种情景下的表现越稳定,物流网络的鲁棒性越强[15]。且ωr的变动范围为(0,11.48%),说明即使最坏的情景发生,鲁棒优化解在该种情景下的遗憾值也能控制在p值内。而ωs的变动范围为(0.03%,27.31%),在不同情景下的表现存在较大差距,这表明当某种特定情景发生时 (如6∗8算例中的情景1),随机规划解会产生较大的遗憾值,极大地增加了决策风险。此外,虽然随机规化模型的总成本较鲁棒优化模型低,但Zr-s最大为4.55%,成本差距并不明显。因此,相比随机规划模型,鲁棒优化模型所求得的物流网络优化方案具有更强的鲁棒性,能够有效减少决策风险。
为了验证本文所设计遗传退火算法 (GASA)的有效性,以上述数据为基础,设计不同规模算例,对比GASA与遗传算法 (GA)及模拟退火算法 (SA)的求解结果,经过调参实验后,设置遗传算法的种群规模popsize为100,交叉概率Pc为0.9,变异概率Pm为0.06;模拟退火算法的初温T0为9000,降温系数low为0.99,每一温度下的最大迭代次数maxnum为400。对于规模较小的算例1、2、3,遗传算法和模拟退火算法的最大迭代次数均设置为1000,对于规模较大的算例4,设置最大迭代次数为5000。设置遗传退火算法的相关参数取值与上述两种算法保持一致,∂=0.05,θ=1000000,在MATLAB R2017b环境下测试,每一算例均运行10次,3种算法在情景1的运行结果对比如表6所示。
表6 算法运行结果对比
从表6可以看出,GASA在10次中的最优值、最差值及平均值均小于SA与GA,虽然在融合了两种算法的基本操作步骤后,GASA的平均运行时间大于SA与GA,但能够满足物流网络优化问题对时效性的要求。通过以上的对比分析可证明本文所设计遗传退火算法的有效性。
保持其他参数取值不变,在6∗10规模网络(即i1~i6、j1~j10) 下对遗憾值限定系数p及碳税φ分别进行敏感性分析。
(1)遗憾值限定系数p敏感性分析
p值对目标函数值的影响如图3所示。由图3可知,随着p值的减少,目标函数值逐渐增大,这意味着物流网络的鲁棒性与成本之间存在效益背反规律,若要增强网络的鲁棒性则需付出更高的成本。当p<0.1430时,算法未能找到可行解,这是因为p存在一个下限值pmin,且当p接近pmin时,目标函数值对p的变化较为敏感,此时增强网络的鲁棒性所付出的成本相对更高。
而p>0.1678时,目标函数值不再发生变化,这是因为p存在一个上限值pmax,当p>pmax时,约束条件 (24)将失效,此时即等同于随机规划模型。因此,决策者需要在鲁棒性和成本之间均衡选择,通过动态调整参数p调节整个物流网络的鲁棒性水平。
图3 遗憾值限定系数敏感性分析
(2)碳税值φ敏感性分析
调节碳税φ的大小,得到碳排放量和总成本随碳税φ的变化情况如图4所示。从图4可以看出,随着碳税值的增加,碳排放量呈现阶梯型递减特征,目标函数值则呈阶段性线性增长趋势——当碳税值在某一区间 (如φ=0~0.11,0.12~0.20,0.21~0.30时)内增加时,鲁棒优化方案及碳排放量均不发生变化,在该区间内,目标函数值呈线性增长,其斜率即对应鲁棒优化方案的碳排放量。但当碳税值增加到某一临界点 (φ=0.12,0.21)时,碳排放量突然下降。这是由于此时改变优化方案所减少的碳税成本大于所增加的其他项成本,企业将选择碳排放量更少的优化方案 (如表7所示),目标函数值在该阶段的斜率也随之减小。因此,碳税政策制定者需要根据这一规律,结合经济发展情况进行相应的调研分析,找到碳排放减少的临界点设置碳税值,使其在有效减少碳排放量的同时尽可能的降低企业的经济负担。
图4 碳税值敏感性分析
表7 不同碳税值下的优化方案
根据以上分析,本文得出如下结论:
(1)建立了考虑碳税和多种需求量情景的整车物流网络鲁棒优化模型,通过与随机规划模型对比,证明鲁棒优化模型能够有效减少决策风险,但提高网络的鲁棒性必须付出更高的成本。
(2)设计了遗传退火算法对模型进行求解,并通过不同规模的算例实验证明其较参数取值相同下的遗传算法和模拟退火算法性能更优。
(3)灵敏度分析实验表明,总成本随着碳税值的增加逐渐增加,碳排放量随着碳税值的增加呈现明显的阶梯状下降特征,碳税政策制定者可利用这一规律在碳排放量变化的临界点设置碳税值,使其既能达到减排目的又能最大限度地减轻企业的经济负担。