0 引言

随着人们应对全球变暖意识的提高，低碳、绿色相关问题受到广泛关注。中国提出“碳中和”、“碳达峰”的双碳目标，明确了各行业的碳减排要求。交通运输行业作为三大碳排放来源之一，是节能减排的重点领域。多式联运是一种低碳高效的货物运输方式，能够通过优化运输结构以提高运输效率，降低物流成本，减少碳排放。碳税是控制碳排放的市场化政策工具，通过内生化企业碳排放的外部性引导企业低碳投资，驱动低碳创新^[1]，被很多国家和地区所使用。中国2017年颁布的《“中国碳税税制框架设计”专题报告》提出了碳税实施框架、推出时间等方案^[2]，但并未明确征收机制、税率水平等具体问题。因此，通过路径优化研究分析不同碳税征收机制及其税率水平在多式联运中的减排效果与影响，能够为碳税制度设计者和具体决策者(货运企业)提供参考，对于发展低碳交通、推进双碳目标的实现具有重要的理论和实践意义。

国内外学者在多式联运路径优化方面取得了丰富的研究成果^[3-5]，并随着绿色低碳等发展理念的提出开始关注低碳多式联运相关问题。关于路径优化问题，现有研究可分为2个方面：一是基于实证研究分析多式联运路径优化在低碳环境下的优越性，如Santos等^[6]以比利时为例分析碳排放控制下铁路多式联运企业的竞争优势；金琳等^[7]以西南地区集装箱多式联运为例，对单一运输和多式联运方式进行SWOT分析，验证多式联运在环保性、经济性方面的优势；李晓东等^[8]以中国东北地区集装箱多式联运为例进行低碳运输的实证研究，确定运输方案并证明多式联运的优势；二是基于不同碳排放政策建立多式联运路径优化模型并分析其影响，如Hua等^[9]、张旭等^[10]研究碳交易政策下的路径选择问题，分析碳排放配额限制和碳交易价格对决策的影响；Bouchery等^[11]构建多目标动态多式联运网络模型，确定强制碳排放控制下的运输方案；蒋琦玮等^[12]基于统一碳税机制进行多式联运路径优化，前者以成本和时间最小构建双目标模型，后者分析不同碳税率变化对决策的影响；崔娥英等^[13]比较统一碳税机制和分段累进碳税机制对多式联运网络优化决策的影响；程兴群等^[14]、袁旭梅等^[15]通过模型分析探讨强制碳排放、碳税、碳交易和碳补偿4种政策对多式联运碳排放及成本的影响。受市场波动、订货提前期和天气环境等的影响，实际多式联运过程中的货物需求存在明显的不确定性，学术界对此采用不同方法进行研究：Sun等^[16]、张旭等^[17]针对多式联运过程中需求不确定的情形，分别构建混合整数非线性模型和具有遗憾值约束的随机优化模型并求解；陈汨梨等^[18]针对运输速度不确定问题，结合随机期望值模型与机会约束规划理论构建路径优化组合模型并求解；Ziaei等^[19]考虑需求的不确定性构建多目标路径优化模型，并采用鲁棒优化方法对随机模糊变量进行处理。

综上，现有研究在低碳多式联运路径优化方面已取得了一定的成果，但还有待进一步完善在不确定的决策环境下同时分析不同碳税机制影响的多目标路径优化研究。鉴于此，本研究考虑由于季节性变化或突发性补货等意外因素导致的需求难以预测的问题，在不确定需求下构建总成本和总碳排放量最小的双目标低碳多式联运路径优化模型，基于随机期望值进行模型转换，并设计改进的非支配排序遗传算法对模型求解。通过算例分析模型和算法的有效性，计算运输路径与方案，并进一步探讨不同碳税机制对总成本和总碳排放量的影响，为政策设计和多式联运企业的决策制定提供参考。

1 问题描述

假设某运输企业拟承运一批不确定需求量的货物，该任务需在规定时间内由供应地O运输至需求地D，其间经过若干节点进行中转，相邻两节点间可采用公路、铁路、水路等不同方式运输。已知各运输方式对应的运输能力、距离、速度、单位费用及碳排放量。以总成本和总碳排放量最小为目标确定运输路径与方式，明晰不同碳税机制的影响。模型中的总成本包括各节点间的直接运输成本、转运成本、时间成本及基于碳税的碳排放成本。同时考虑多式联运的实际情况，对模型做出假设为：(1)假设货物在运输过程中不可分割，即同一批货物在两相邻节点间只能选择一种方式进行运输；(2)假设同一节点处最多发生一次运输方式转变且节点具备足够的中转能力；(3)假设不同运输方式的货运能力均能满足客户的货运量要求；(4)不考虑运输过程中天气、货损、设备故障等突发事件造成的等待时间、费用等的影响。研究中涉及的主要参数及含义如表 1所示。

2 模型构建

本研究的低碳多式联运路径优化问题不仅要考虑如何表示需求的不确定性并加入目标函数，还要考虑不同碳税机制变化对总成本和碳排放量的影响。具体建模思路为：设需求量Q服从均值为μ、方差为σ²的正态分布，即Q~N(μ, σ²)。首先构建需求不确定下基于不同碳税机制的多式联运路径优化基础模型，其次借助随机优化理论将基础模型转换为随机期望值模型。

2.1 不同碳税机制下的模型构建

碳税是环境税的一种，无论是实践还是学术领域均认为其是控制碳排放的有效手段，但目前还未对具体的征收机制和税率水平达成共识。本研究分别在统一碳税机制和分段累进碳税机制下构建需求不确定的多式联运路径优化模型。

2.1.1 统一碳税机制下的路径优化模型(模型Ⅰ)

统一碳税机制下，按照统一的碳税率对运输中的碳排放量征收碳税，对应的模型Ⅰ如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

式(1)为总成本最小的目标，包括运输成本C₁，由货运量、节点间运输距离和单位运输成本的乘积表示；中转成本C₂，由货运量和某一节点运输方式转换对应的单位成本的乘积表示；时间成本C₃，由运输时间和转运时间构成，为提前到达的仓储成本或延迟到达的惩罚成本；碳排放成本C₄，由碳排放量和碳税率的乘积表示。式(2)为碳排放量最小的目标，包括运输过程碳排放量E₁^e，由货运量、节点间单位碳排放量和运输距离的乘积构成；转运碳排放量E₂^e，由货运量和某一节点运输方式转换时相对应的碳排放系数的乘积表示。式(3)表示货物需求量服从均值为μ、方差为σ²的正态分布。式(4)~ (6)为与假设条件相关的约束，其中式(4)表示节点i, j间至多允许采用一种运输方式；式(5)表示节点i处至多进行一次运输方式转换；式(6)约束了运输的连续性，以防止运输方式转换时各路径发生冲突。式(7)和式(8)限定决策变量为0-1变量。

2.1.2 分段累进碳税机制下的路径优化模型(模型Ⅱ)

分段累计碳税机制是按照不同碳排放量确定不同税率的碳税征收模式，具体做法是将碳排放量划分为若干等级并针对每一等级设计对应的碳税率。在现实中，电价、水价等均采用此种收费机制，以响应节能号召，抑制浪费，碳排放量越大对应的税率越高^[13]。本研究结合现有文献中的分段函数形式^[13]及个人所得税中速算扣除数的概念^[17]，将分段累进碳税机制下的碳排放成本C₄表示为：

(9)

式中β_m为速算扣除数，是为解决超额累进税率分级计算税额时的复杂技术问题，而预先计算出的数据为：

(10)

根据研究问题特点，设E₀=0，β₁=0。例如：当m=2时，碳排放量E^e处于第2阶段，即E^e∈ (E₁, E₂]，则C₄=α₁E₁+(E₂－E₁)α₂，表示不超过E₁部分按低碳税率α₁征收，超出E₁部分按较高碳税率α₂征收。

分段累进碳税机制下的模型Ⅱ中，碳排放量最小的目标函数不变(见式(2))，总成本最小的目标函数见式(11)。与式(1)相比，碳排放成本C₄的计算发生了变化。

(11)

(12)

同时，式(3)~ (8)，(10)成立。其中，式(12)表示m级碳税对应的碳排放阈值。

2.2 基于随机期望值的模型转换

研究结合随机优化理论，通过极大化不确定目标的期望值，将基础模型转换为随机期望值模型。以模型Ⅰ为例，当R→∞时，利用Monte Carlo模拟和大数定律得到的运输成本、转运成本、时间成本、运输碳排放量、转运碳排放量的期望值形式E (C₁)，E (C₂)，E (C₃)，E (E₁^e)，E (E₂^e)可以表示为式(13)~(17)形式。

(13)

(14)

(15)

(16)

(17)

转换后的模型Ⅰ为：

(18)

(19)

同时式(3)~(8)成立。其中，式(18)表示总成本的期望目标函数，式(19)表示总碳排放量的期望目标函数。类似地，模型Ⅱ可以基于上述方法转换为确定型模型。

3 改进的非支配排序遗传算法

研究所构建的模型Ⅰ和模型Ⅱ均为双目标0-1优化模型，各目标相互冲突，通常无法求得唯一最优解，而只能得到非劣解集，即帕累托最优解集^[20]。相比传统遗传算法，NSGA-Ⅱ算法引入了快速非支配排序、拥挤度计算和精英策略，降低了计算复杂度，提高了搜索能力，是目前较为成熟的求解多目标问题的算法。但其易陷入“早熟”，为避免该缺陷，本研究设计了改进的非支配排序遗传算法(Ⅰ-NSGA-Ⅱ)。该算法在生成初始化种群时加入拓扑排序规则，采用改进的循环拥挤排序方法进行选择操作，在精英策略中引入比例法。

3.1 染色体编码与种群初始化

(1) 染色体编码

编码是遗传算法的第1步，直接影响遗传算法的运算结果。多式联运路径优化问题需要确定运输的路径和方式，因此研究采用双层编码结构：第1层为运输路径x_ij^k的编码，第2层为运输方式y_i^kl的编码，如式(20)，二者均以实数方式编码。其中，染色体总长度为2N-1，即共N个节点采用N-1种运输方式：

(20)

(2) 基于拓扑排序的种群初始化

由于多式联运网络中的各节点是按照一定顺序排列的，为了保证解的可行性，且考虑到研究特点，采用拓扑排序规则生成初始化种群，即基于多式联运网络生成的有向无环图设计拓扑排序序列，并据此进行染色体的初始化，以抑制非法方案的产生。

3.2 适应度函数

考虑到遗传操作的特点，将总成本和碳排放量最小的目标函数进行最大化转换，构造适应度函数。设C (i)和E^e (i)分别为个体i对应的总成本和碳排放量，则适应度函数可相应地表示为和。

3.3 快速非支配排序与拥挤度计算 3.3.1 快速非支配排序

根据计算得到的适应度值，采用快速非支配排序对种群进行排序，具体步骤为：

(1) 为种群中的每一个解p(p∈R)定义2个变量N_p和S_p，分别为解p被支配的数量和解p支配的解的集合；

(2) 找出所有N_p为0的解，形成非支配解集F₁，令其非支配等级为1，为集合中的所有解设置相同的非支配排序；

(3) 令当前集合F₁中的每一个解所支配的解的集合为S_q，遍历S_q中的每个解，并执行N_q=N_q-1，若N_q-1=0则将个体q存入集合F₂中，以F₂作为当前集合；

(4) 重复步骤(2)和步骤(3)，确定每个解集的非支配等级及每个解的非支配排序；

(5) 将找出的所有非支配解集进行排序，{F₁>F₂>…>F_n}表示集合F₁中的个体最优，其他依次降低。

3.3.2 拥挤度计算

在快速非支配排序中，个体所在的非支配层越低则个体越优。然而，若个体处于同一非支配层，则需要进一步进行拥挤度比较。拥挤度表示空间中个体的密度值，可以用个体周围不包括其他个体的长方形表示(见图 1)。种群的拥挤度越大，个体n与周围个体的相对距离就越大，解的多样性越好。拥挤度n_d计算为：

(21)

式中f_h^max和f_h^min分别为个体集合中第h个目标函数的最大值和最小值。

3.4 遗传操作 3.4.1 基于循环拥挤排序方法的选择算子

传统NSGA-Ⅱ算法一般采用锦标赛选择法，随机选择2个个体比较其非支配排序等级，等级低的进入下一代种群，若排序等级相同，则对每个非支配层中的个体进行拥挤度计算。但事实上仅采用拥挤度大小来保持群体多样性并不能客观反映个体间的真实拥挤程度^[21]。因而引入改进的循环拥挤排序计算方法^[22]，以剔除的方式保留个体，并要求循环剔除排序后的个体数等于要选择的个体数，以使最终得到的Pareto解均匀分布。具体步骤为：

(1) 确定种群中非劣前沿等级上的总个体数即非支配解个数为n，要求从该等级中筛选i个个体；

(2) 计算n个个体的拥挤度，剔除拥挤度最小的个体；

(3) 判断剔除后的个体数是否为i，若是则结束，否则，进入下一步；

(4) 逐个剔除拥挤度最小的个体，每剔除一个个体，重新计算剩余个体的拥挤度，并重复上一步，直至结束。

3.4.2 交叉与变异

研究采用单点交叉法对染色体进行交叉操作，即按照固定的交叉概率对父代种群中的2条染色体的部分基因做随机交叉运算，形成新的编码串并根据拓扑排序规则和适应度函数值判断可行性。变异可以通过基因突变产生新的个体，同样地，采用单点变异法进行变异操作，即在一条染色体上选择2个变异点，按照一定的变异概率对变异点之间的区域做逆转操作，从而产生新个体以保证种群的多样性。

3.5 基于比例法的精英策略

为保证精英个体不被遗漏，同时降低计算复杂度，NSGA-Ⅱ中加入精英策略，但传统精英策略可能导致种群中子代个体被拒的概率随迭代次数的增加而增加，使算法陷入局部最优。因此，设计基于比例法的精英策略^[23]，即设Gen为当前迭代次数，最大迭代次数为NL max，通过在第NL max/3代时采用比例法并持续到NL max代，其余迭代中则采用传统精英策略，保证父代个体所占比例随迭代次数的增加而呈线性下降趋势，增加子代个体在生成新种群过程中被选择的概率。综上，研究所设计的改进的非支配排序遗传算法(Ⅰ-NSGA-Ⅱ)的流程如图 2所示。

4 算例分析 4.1 算例数据

某多式联运运输企业要将一批货物由包含15个节点的联合运输网络从出发点O运输至目的点D，客户期望的时间窗为[80 h，90 h]，货物提前到达的单位存储成本为18元/h·t，延迟到达的单位惩罚成本为30元/h·t。参考各国碳税征收标准，设不分段时的碳税率为0.2元/kg，分为两阶段时的碳税率为(0.2, 0.3] 元/kg，每阶段的碳排量下限为(0，4 000) kg。不同节点间可采用公路、铁路、水路3种方式运输，对应的运输距离见表 2。表中(a，b，c)分别表示公路、铁路、水路对应的运输距离，“-”表示两节点间没有对应的运输方式。

查阅货运服务的运价率表、《IPCC国家温室气体清单指南》和中国交通年鉴等资料，确定各运输方式的平均速度、单位运输成本、单位碳排放量(见表 3)。表中(a，b，c)分别表示运输节点间距离在“500 km以内”，“500~1 000 km”，“1 000 km以上”这3个范围时的单位运价。通过查阅相关文献资料和各港口、货运中转站等的具体情况，得到不同运输方式的单位转运成本、时间和碳排放系数(见表 4)。由于运输计划的超前性，需求企业的需求量难以确定。根据现实中货运量的分布特征及相关文献中的设计^{[10, 17]}，设货物运输量Q服从均值为100、方差为64的正态分布，即Q~N(100, 64)。

4.2 结果分析 4.2.1 算法参数与运行结果

根据所建立的模型，取随机次数R=100，采用改进的非支配排序遗传算法，以Matlab2019程序编程对模型求解。算法参数设置为：种群规模为100个，迭代次数为200次，交叉概率为0.1，变异概率为0.8。以固定碳税为例，目标函数Pareto最优解集的进化过程如图 3所示，逐步减小并在到达一定迭代次数后最终收敛，稳定在最优解集附近保持不变。

为分析不同碳税机制对总成本和碳排放量的影响，分别在固定碳税和分段累进碳税2种情况下计算得到Pareto最优解集，结果如表 5所示。每种机制的2个方案分别对应总成本和碳排放量最小。

在表 5中的2种碳税机制下，方案1的总成本均小于方案2，但方案1的碳排放量均大于方案2，产生这一“此消彼长”现象是由于双目标优化问题中并不存在唯一最优解使2个目标同时达到最优。通常一个目标函数达到最优时，另一个目标函数结果就会变差。具体来看，在分段累进碳税机制下，方案1和方案2的总成本和碳排放量相较固定碳税的结果均有所下降，特别是方案2，其总成本减少了4.75%，总碳排放量降低了16.13%。这说明分段累进碳税机制具有更加明显的减排效果。

4.2.2 算法有效性检验

在算法有效性方面，以固定碳税率为例，分别采用传统NSGA-Ⅱ算法和改进NSGA-Ⅱ算法求解模型，各运行10次，2种算法的运行时间与求解结果如表 6所示。其中，Ⅰ-NSGA-Ⅱ算法运行时间比传统NSGA-Ⅱ算法高出5.72 s，但其运行结果对应的平均最优成本相比传统算法减少了1.47%，平均最优碳排放量则降低了6.85%。这是由于相较于传统NSGA-Ⅱ算法，改进的NSGA-Ⅱ算法由于采用基于循环拥挤排序法的选择算子和基于比例法的精英策略，扩大了搜索空间与范围，能够获得更加优秀的个体与方案。

4.2.3 模型适用性分析

以固定碳税为例，采用传统遗传算法对只考虑总成本的单目标优化模型求解，以说明上述考虑总成本和总碳排放量的双目标模型相比单目标模型对于双碳背景的适用性，具体结果如表 7所示。

对比表 5和表 7，与单目标模型相比，双目标优化中的总成本由64 887.53元变为66 293.89元，增长2.17%；碳排放量由4 617.06 kg变为4 021.95 kg，降低9.48%。可以看出双目标函数下单位成本所产生的碳排放量下降更加明显，采取同时考虑总成本和总碳排放量的双目标策略，可以使决策者更加清晰地认识到企业仅需略微提高成本即可取得一定的减排效果，更加适合双碳背景下的运输场景。

4.3 不同碳税机制影响分析

由表 5可见，相较于固定碳税机制，分段累进碳税机制下碳排放量下降显著，这是由决策者的高阶碳税率厌恶促使的，为使碳排放量小于低阶碳税率对应的排放阈值而调整决策方案。由此，分段累进碳税机制是引导和鼓励企业控制碳排放的有效手段。

为进一步分析碳税率变化对决策的影响，使固定碳税率和分段累进碳税率分别在0~0.6之间取值，采用第3节所提出的改进的非支配排序遗传算法，基于模型Ⅰ的统一碳税机制和模型Ⅱ的分段累进碳税机制下的总成本和碳排放量计算公式，借助Matlab2019程序编程计算不同情景下对应的成本与碳排放量，结果如表 8所示。

由表 8可见，在2种碳税机制下，Ⅱ级碳税率与Ⅰ级碳税率(不考虑碳排放成本，税率为0)时的运输方案与碳排放量均无差异，碳税约束效应不明显，企业选择牺牲环境以节约成本；Ⅲ~Ⅳ级碳税率下，碳排放量降低，且在分段累进碳税机制下表现的更加显著；从Ⅴ级碳税率开始，高碳税率仅带来了成本的增加而碳排放量并未发生变化，减排效果不理想。综上，无论在固定碳税还是分段累进碳税机制下，企业的碳排放控制效果均会受到碳税率的影响，应设置合适的碳税率以调动企业减排积极性。

假定每阶段的碳税率仍为(0.2, 0.3]，使碳排放量上限由4 000 kg降低到2 000 kg，分析分段累进碳税机制中低阶碳税率对应的排放阈值变化对决策的影响，结果如图 4所示。由图可见，排放阈值由4 000 kg降低到2 000 kg时，总成本呈现持续上升趋势；排放阈值由4 000 kg降低到3 000 kg时，总碳排放量呈下降趋势，但排放阈值的进一步降低并未带来碳排放量的改变。这是由于企业的碳排放量会受到现有能力及减排技术的限制，在二者不变的条件下若已达到最低碳排放量，即使面对逐渐严苛的减排要求，企业亦束手无策，只能被动承受成本增加带来的压力。由此，应考虑企业现有能力与减排技术水平，确定合适的碳税率与排放阈值。

4.4 不同碳税机制下需求波动影响分析

考虑到实际运输任务中货运量的日常数据及对不同货运量的关注，同时参考相关文献中货运参数的设置及其计算方法^[24]，在期望值不变的前提下令需求方差分别为25，36，49，64，81，分别计算固定碳税和分段累进碳税机制下的运输方案、成本与碳排放量，探究需求波动幅度对决策的影响，结果如表 9所示。

由表 9可见，需求方差即波动性的增加意味着供需双方信息不确定性的增大，会对决策产生影响。无论在固定碳税还是分段累进碳税机制下，路径优化系统均为应对需求不确定性增加产生的影响而分配了更多的能力，导致成本和碳排放量的增长。当方差较小时，固定碳税和分段累进碳税机制下的总成本与碳排放量均差别不大，但随着方差的增加，分段累进碳税机制的优势逐渐凸显，总成本和碳排放量相对较小。这主要来源于决策者在高阶碳税率厌恶的利益驱使下对低碳排放量的主动追求。

5 结论

本研究在双碳背景下研究多式联运路径优化问题，基于改进的非支配排序遗传算法求解需求不确定的双目标路径优化模型，并分析不同碳税机制的影响。通过算例分析了改进的NSGA-Ⅱ算法的优越性及双目标优化模型对低碳运输场景的适用性，求得不同碳税机制下满足总成本和碳排放量最小的相对较优解。同时，比较了不同碳税机制下的路径优化方案，结果显示分段累进碳税机制具有更加显著的减排效果，特别在需求不确定性较高的情况下表现得更加明显。因此，相关部门应考虑企业现有能力与减排技术，确定合适的碳税率与排放阈值。

本研究可为相关部门和多式联运企业提供政策参考与决策依据，但其中仅考虑了单一运输任务，未来将对多任务低碳路径优化问题展开研究。此外，交通运输领域双碳目标的实现，不仅要依托相关部门和多式联运企业，客户的低碳选择同样具有重要作用，考虑客户低碳行为与偏好的路径优化问题将是下一步的研究重点。