阶梯式碳税下煤炭海铁联运网络调运优化研究

Optimization of Coal Allocation on Sea-rail Intermodal Transportation Network with Tiered Carbon Tax

摘要

在大规模海铁联运网络背景下，考虑铁路年运输能力和煤炭港口年吞吐量限制的基础上，引入阶梯式碳税量化碳排放成本，建立考虑碳排放的煤炭海铁联运网络调运方案非线性优化模型。模型以大规模海铁联运网络每条可行路径的年煤炭运输量为决策变量，根据模型特性设计了一种组合式求解算法(Combinatorial Algorithm，CA)，外层迭代算法采用增广拉格朗日乘子法，内层迭代算法采用可行方向法，将非线性规划模型转化为标准Hitchcock模型。采用国家能源集团煤炭海铁联运业务相关数据进行实例分析，并运用CA进行求解。在实例分析的基础上进一步对固定碳税和阶梯式碳税进行了敏感性分析，结果表明，阶梯式碳税上下限在合理度范围内时，采用其得到的最优调运方案在企业总效益和减排效果2个方面优于固定碳税，但当其上下限超出合理范围后，会对企业效益及企业减排积极性产生打击。

Abstract

In the context of a large-scale sea-rail intermodal transportation network, considering the annual transport capacity of railways and the annual throughput of coal ports, this paper introduced a tiered carbon tax to quantify carbon emission costs. A non-linear optimization model for the coal sea-rail intermodal transportation network was built. The model took the annual coal transportation volume along each feasible path in the large-scale sea-rail intermodal transportation network as decision variables. A combinatorial algorithm (CA) was designed based on the characteristics of the model, where the outer iteration algorithm employed the augmented LaGrange multiplier method, and the inner iteration algorithm employed the feasible direction method, transforming the non-linear programming model into a standard Hitchcock model. The instance analysis was conducted by relevant data from the coal sea-rail intermodal transportation business of the China Energy Investment Group Co., Ltd, and the CA was applied for the solution. Based on the instance analysis, sensitivity analyses were further performed on fixed carbon taxes and tiered carbon taxes. The results indicate that within a reasonable range of upper and lower limits for tiered carbon taxes, the optimal transportation scheme derived from it outperforms fixed carbon taxes in terms of overall enterprise benefits and emission reduction effects. However, when the upper and lower limits exceed a reasonable range, it adversely affects enterprise benefits and the motivation for emission reduction.

Graphical abstract

关键词

海铁联运 / 能源供应链 / 阶梯式碳税 / 非线性规划 / 增广拉格朗日乘子法

Key words

Sea-rail Intermodal Transportation / Energy Supply Chain / Tiered Carbon Tax / Nonlinear Programming / Augmented LaGrange Method

0　引言

我国煤炭资源总体呈现北多南少、西多东少的分布形式，而我国煤炭的主要需求单位火力发电厂，除部分分布于内陆铁路沿线的主要城市外，大部分分布于东南沿海等发达地区，这导致我国长距离煤炭运输需求旺盛。考虑到我国东南沿海地区航运系统发达，且长距离铁路运输与航运相比成本高且运量小，因此海铁联运是目前我国主要的煤炭运输方式之一。传统的海铁联运调运方案决策仅考虑运量及运输成本，而在国家深化交通运输供给侧结构性改革，以及实现“碳达峰”“碳中和”目标的背景下，低碳海铁联运煤炭运输的研究具有重要的理论和现实意义^[1]。近年来的文献中有较多低碳多式联运路径优化方面的研究。李珺等^[2]将运输时间、中转时间、客户需求和中转集拼货运量设定为四重混合不确定因素，通过分析受不确定因素影响下的多式联运成本，进行多式联运路径规划。邓学平、梅梦婷等^[3-4]均采用三角模糊数表示多式联运路径优化问题中的不确定需求和不确定时间，不同的是，邓学平采用机会约束规划理论将不确定模型确定化，梅梦婷则是构建考虑碳排放的多目标模型并采用多目标算法求解。袁旭梅等^[5]将不确定的货运需求表示为区间数形式，并基于此构建多式联运路径鲁棒优化模型，并且分析了不同的碳排放政策对模型的影响。Zhang等^[6]针对低碳多式联运路径优化问题，设计了混合鲁棒随机优化模型，并分别针对运输模式、货运需求和运输时间不确定场景下的模型及算法的有效性进行了分析。Jiang等^[7]研究了城市群发展背景下具有二氧化碳减排目标和需求不确定特点的区域多式联运物流网络设计问题。蒋琦玮等^[8]在时间不确定的条件下，针对集装箱多式联运路径优化问题，将碳排放成本内化为总成本，并研究了税率变化对成本的影响。张旭等^[9]将碳交易价格考虑为不确定条件，建立多式联运路径混合鲁棒随机优化模型。刘松、Zhao等^[10-11]分别研究了碳排放限制条件下的冷藏集装箱和普通集装箱多式联运路径优化问题。从以上研究可以看出，在进行多式联运路径优化研究时碳排放是一项非常重要的影响因素，有的研究是将碳排放内化为总成本的组成部分，有的分析对比了不同的碳排放政策对路径优化结果的影响。然而，对于碳排放价格的刻画，这些研究采用的都是固定费率的碳交易价格或者碳税。相关针对不同碳排放政策间替代效应的研究结果表明，就目前的实际情况来看，碳排放权交易价格往往过低，对促进碳减排的效果差强人意，而碳税的加入恰好能弥补这一不足。翁智雄等^[12]通过设定不同的税率测算出2020年的碳减排效果，发现税率越高，碳减排效果越好。王金南等^[13]基于可计算一般均衡模型(Computable General Equilibrium，CGE)开发的能源政策综合评价模型——能源经济模型，针对低、中、高3种碳税税率对我国国民经济、能源消费及碳减排的影响进行了情景模拟，提出宜在2012年按中等税率(20元/t)征收碳税，并认为税率应呈阶梯式上升。姚昕等^[14]采用动态DICE模型，研究发现我国最优碳税额度随时间推移而逐步上升，由2008年的7.31元/t上升至2020年的57.61元/t。因此可以看出，在合理阶梯高度范围内，阶梯式碳税相较于固定碳交易价格和碳税相比，对企业的碳减排有明显的促进作用。而一旦阶梯式碳税的阶梯过高，会有增加政策阻力风险^[15]。因此，为了探索更为合理的碳排放政策形式，在多式联运运输方案制定过程中引入阶梯式碳税具有重要的现实意义。

同时，针对煤炭运输优化问题，Gupta等^[16]基于层次分析法和数据包络分析技术构建多目标优化模型，同时优化采矿业车辆运输路线和类型的选择。Li等^[17]构建了一个涵盖我国30个省级行政区域的能源系统模型，分析了省际煤炭运输的碳排放量，同时预测了不同碳减排政策场景下各地区煤炭运输污染排放的变化趋势。Liu等^[18]以煤种、煤质、运输条件等因素为评价依据构建煤炭供需匹配度评价体系，进而构建以匹配程度最大和营业额最小为目标的煤炭运输分配优化模型。郭波等^[19]考虑包括铁路货运通道规划造价成本在内的区域全局运输成本，提出了一种煤炭铁水联运运输通道规划方法。以上针对多式联运低碳化转型和煤炭运输优化的研究大多是针对路径优化的研究，即多式联运运输方式的选择，其模型中会考虑单次运输的运量以计算运输成本及碳排放量，却没有考虑运输线路的运输能力和线路间的运力平衡。尤其是在铁路运输段，实际运营过程中各区间运输能力有着严格限制，若不考虑可行路径间的运力平衡，一旦超出限制将面临运输方案无法实施的后果。因此，在海铁联运网络的建模中，将着重考虑铁路区间运能限制和煤炭港口的吞吐量限制，以确保调运方案在实际运营过程中切实可行。

综上所述，研究在大规模海铁联运网络背景下，考虑煤矿与电厂间供需平衡约束和铁路与港口的运输能力限制约束，同时引入阶梯式碳税策略，构建海铁联运煤炭调运方案非线性优化模型，并针对模型特征设计一种组合式算法(Combinatorial Algorithm，CA)对模型进行求解，该算法分为内外2层迭代算法，外层迭代算法采用增广拉格朗日乘子法，内层迭代算法采用可行方向法。最后，通过算例分析验证模型和算法的有效性，并比较固定碳税和阶梯式碳税策略下的煤炭调运方案碳排放量变化趋势。

1　问题描述及符号定义

1.1　问题描述

煤炭由煤矿开采出后，经过短驳运输到达对应铁路装车站。经装车站装车后，煤炭正式进入海铁联运网络开始运输。装车后的煤炭分为2类去向，一是经过铁路运输至沿线各电厂卸车，完成煤炭运输，称为直达煤；二是运往各煤炭转运港，准备进行第二阶段的航运运输，称为下水煤。下水煤经由港口转运后，由货运船舶运输至东南沿海各电厂，进而完成煤炭海铁联运过程。直达煤的运输会占用铁路运输能力，从而影响下水煤的调运量。煤炭海铁联运网络示意图如图1所示。

因此，研究的问题是在考虑线路运能约束和阶梯式碳税的条件下，寻找满足煤矿和电厂煤炭供需平衡要求且煤炭生产企业总效益最大化的煤炭联运路径方案及各条路径上的年运输煤炭总量。

1.2　海铁联运网络符号定义

G\left(N,L\right)表示海铁联运网络，其中N=\left\{n_i,n_j,n_q,n_g\right\}表示网络中运输节点集合，由装车站、电厂、港口和其他节点组成，n_i=\left\{n_{i_{\mathrm{1}}},n_{i_{\mathrm{2}}},\cdots \right\}表示装车站集合，n_j=\left\{n_{j_{\mathrm{1}}},n_{j_{\mathrm{2}}},\cdots \right\}表示电厂集合，n_q=\left\{n_{q_{\mathrm{1}}},n_{q_{\mathrm{2}}},\cdots \right\}表示港口集合，n_g=\left\{n_{g_{\mathrm{1}}},n_{g_{\mathrm{2}}},\cdots \right\}表示除装车站、电厂、港口之外的网络节点。L=\left\{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\left|n_{\mathrm{1}},n_{\mathrm{2}}\in N\right.\right\}表示煤炭海铁联运网络路段集合(相邻节点)，其中l_{n_{\mathrm{1}},n_{\mathrm{2}}}表示相邻节点n_{\mathrm{1}}和n_{\mathrm{2}}之间的路段，因航运线路中间不经过其他港口，故其l_{n_{\mathrm{1}},n_{\mathrm{2}}}即表示由北方出发港到东南沿海到达港的整体路段。

1.3　阶梯式碳税

阶梯式碳税是指税率随着碳排放量的增加而呈现阶梯式增高^[15]。

阶梯式碳税的阶梯高度为相邻税率之间的差值，则阶梯高度\mathrm{\Delta }f\left(\rho \right)和碳税T\left(\rho \right)为

\mathrm{\Delta }f\left(\rho \right)=f\left({\rho }_a\right)-f\left({\rho }_{a-\mathrm{1}}\right)（1）

T\left(\rho \right)=\sum f\left({\rho }_a\right)\cdot {\rho }_a（2）

式中：\rho为除去初始碳配额外需要支付碳税的碳排放量；\alpha为能源企业初始碳配额，t。

考虑到以上阶梯式碳税为分段函数，导致在模型求解时需要先判断碳排放量所处区间。因此，为求解方便，这里将阶梯式碳税的区间取极限处理，即\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{\Delta }\rho \to \mathrm{0}，使得碳税的阶梯式趋于递进式，碳税的税率函数趋于连续的函数。

\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{\Delta }\rho \to \mathrm{0}条件下的碳税T\left(\rho \right)公式如下。

T\left(\rho \right)={\int }_{\mathrm{0}}^{\rho }f\left(\rho \right)\mathrm{d}\rho（3）

取极限处理的阶梯式碳税如图2所示。

如图2所示，求极限之前的原始阶梯型碳税呈离散形势分布，并且随着碳排放量的升高，\mathrm{\Delta }f\left(\rho \right)也会逐渐增加，而求极限之后的阶梯式碳税则呈连续分布，并随着碳排放量的增加，碳税曲线的斜率越来越大。当离散型阶梯式碳税转化为连续性阶梯式碳税后，会出现一定程度的碳税变化，如图2阴影部分所示，当碳排放越大并且越靠近离散型碳税区间边缘时，2种方式的碳税差异越大。但是，如果采用离散型碳税，就需要对每一条运输记录的碳排放量进行识别分类并累加碳税，这将极大地增加运算的时间和复杂度。因此，研究选择采用连续型阶梯式碳税，采用公式⑶计算碳税。

2　阶梯式碳税下煤炭海铁联运网络调运优化模型

2.1　模型假设条件

（1）模型优化范围由煤炭在装车站装车开始至船舶到达电厂结束。

（2）每个煤矿仅对应一个装车站，但一个装车站可对应多个煤矿。每个电厂仅对应一个港口，但一个港口可对应多个电厂。

（3）模型仅计算大宗煤种的年运输总量，不考虑价格差异较大的特种煤。

（4）只在港口进行运输方式的转换，且至多进行一次运输方式的转换。

（5）在运输过程中，煤炭的属性状态不发生变化，没有货损情况。

2.2　模型参数及变量定义

（1）海铁联运网络参数。A_i表示装车站n_i对应煤矿年可装车煤炭总量，t；B_j表示电厂n_j年煤炭需求，t；c_i表示装车站n_i对应煤矿年平均生产成本，元/t；p_j表示电厂n_j购买煤炭年平均售价，元/t；c_q表示港口n_q的平均转运成本，元/t；p_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}表示铁路区间或航线l_{n_{\mathrm{1}},n_{\mathrm{2}}}单位里程单位重量煤炭平均运价，包括铁路运价和船舶运价，元/t；s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}表示铁路区间或航线l_{n_{\mathrm{1}},n_{\mathrm{2}}}的长度，km；d_q表示港口n_q年煤炭吞吐量上限，t；d_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}表示铁路区间或航线l_{n_{\mathrm{1}},n_{\mathrm{2}}}年运能上限，t；\lambda表示煤矿生产煤炭总量低于电厂需求时，需求满足比例。

（2）可行路径属性参数。因网络中包含铁路沿线电厂和东南沿海电厂2类目的地，故模型中可行路径也包括2类：①a路径——煤矿→装车站→铁路沿线电厂；②b路径——煤矿→装车站→港口→东南沿海港口→东南沿海电厂。其中a路径虽然不采用海铁联运，但是其运输会占用海铁联运中铁路段的运能，因此需综合考虑。设集合R=\left\{r_{i,j}^k\left|i\in n_i,j\in n_j\right.\right\}为煤炭生产点到销售区域的全部可行路径集合，其中，r_{i,j}^k表示从n_i到n_j的第k条可行路径。设0-1变量y_{r_{i,j}^k,i}表示路径r_{i,j}^k是否以装车站n_i为起始点，是为1，否则为0。设0-1变量y_{r_{i,j}^k,j}表示路径r_{i,j}^k是否以电厂n_j为终到点，是为1，否则为0。设0-1变量y_{r_{i,j}^k,l_{n_{\mathrm{1}},n_{\mathrm{2}}}}表示路径r_{i,j}^k是否经过路段l_{n_{\mathrm{1}},n_{\mathrm{2}}}，是为1，否则为0；设0-1变量y_{r_{i,j}^k,q}表示路径r_{i,j}^k是否经过港口n_q，是为1，否则为0。

（3）决策变量。设决策变量x_{r_{i,j}^k}，表示路径r_{i,j}^k上的年煤炭运输总量，t。

2.3　目标函数

煤炭海铁联运调运模型是以能源企业煤炭业务总效益最大化为目标。其中，利润部分由煤矿向电厂售卖煤炭所得组成，能源企业煤炭业务成本包括煤矿生产成本、铁路运输成本、船舶运输成本、港口转运成本和碳排放成本。

Z_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{e}}=\sum _{r_{i,j}^k\in R}{p}_j{x}_{r_{i,j}^k}（4）

Z_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\sum _{r_{i,j}^k\in R}{c}_i{x}_{r_{i,j}^k}（5）

Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{p}_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}（6）

Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{f}\mathrm{e}\mathrm{r}}=\sum _{r_{i,j}^k\in R}{c}_q{y}_{r_{i,j}^k,q}{x}_{r_{i,j}^k}（7）

对于碳排放成本，以引入阶梯式碳税下的碳排放成本函数进行表示。此处令\alpha为能源企业初始碳配额，t；{\theta }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}为在铁路区间或航线l_{n_{\mathrm{1}},n_{\mathrm{2}}}上运输1 t煤炭的单位碳排放量，kg/(t·km)；\rho为除去初始碳配额外需要支付碳税的碳排放量。

\rho =\left(\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{\theta }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}\right)-\alpha（8）

为求解方便，这里将阶梯式碳税的区间取极限处理，即\mathrm{l}\mathrm{i}\mathrm{m}\mathrm{\Delta }\rho \to \mathrm{0}，则该条件下，碳税T\left(\rho \right)公式如下。

T\left(\rho \right)={\int }_{\mathrm{0}}^{\left(\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{\theta }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}\right)-\alpha }f\left(\rho \right)\mathrm{d}\rho（9）

公式⑷表示煤矿售煤所得总利润，元；公式⑸至公式⑺和公式⑼为能源企业煤炭业务成本，元；公式⑸为煤矿生产成本，元；公式⑹为铁路运输和船舶运输总成本；公式⑺为港口转运成本；公式⑼为阶梯式碳排放成本，元。

因此，该模型的总目标函数可以表示为

［M 1］

\mathrm{m}\mathrm{a}\mathrm{x}Z\left(x\right)=Z_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{e}}-Z_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}-Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}-Z_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{f}\mathrm{e}\mathrm{r}}-T\left(\rho \right)（10）

2.4　约束条件

煤炭海铁联运调运方案必须满足生产-需求约束和运输能力约束。

\sum _{r_{i,j}^k\in R}{y}_{r_{i,j}^k,i}{x}_{r_{i,j}^k}=A_i，\forall i\in n_i（11）

\sum _{r_{i,j}^k\in R}{y}_{r_{i,j}^k,j}{x}_{r_{i,j}^k}=\lambda B_j，\forall j\in n_j（12）

\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}\le d_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}，\forall l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L（13）

\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},q}{x}_{r_{i,j}^k}\le d_q，\forall q\in n_q（14）

x_{r_{i,j}^k}\ge \mathrm{0}，\forall r_{i,j}^k\in R（15）

公式⑾和公式⑿为生产-需求约束，即各煤矿运出的煤炭总量必须符合煤矿的煤炭供应能力，同时运到各电厂的煤炭也必须满足各电厂的需求。公式⒀和公式⒁为运输能力约束，即在铁路运输部分各区间的煤炭运输总量不能超过铁路区间的运输能力，到达港口的煤炭量也不能超过港口的最大年吞吐量。公式⒂为决策变量非负约束。

3　模型分析与内外双层迭代算法设计

研究所构建的模型引入了阶梯式碳税来量化碳排放成本，但由于阶梯式碳税在目标函数中为积分形式，因此模型从原不考虑碳税的线性规划模型转变为了带有路段和港口能力限制约束的非线性规划模型。因此，每条可行路径上分配的调运量由铁路区间通行能力约束、港口吞吐量约束和随碳排放量呈非线性递增的阶梯式碳税共同影响决定。

首先，[M 1]原模型转化为最小化形式

[M 1]

\mathrm{m}\mathrm{i}\mathrm{n}\underset{x\in \mathrm{\Omega }}{Z_{\mathrm{1}}\left(\mathit{x}\right)}=\sum _{r_{i,j}^k\in R}{c}_i{x}_{r_{i,j}^k}+\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{p}_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}+\sum _{r_{i,j}^k\in R}{c}_q{y}_{r_{i,j}^k,q}{x}_{r_{i,j}^k}+{\int }_{\mathrm{0}}^{\left(\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{\theta }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}\right)-\alpha }f\left(\rho \right)\mathrm{d}\rho -\sum _{r_{i,j}^k\in R}{p}_j{x}_{r_{i,j}^k}(16)

s.t. \sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}\le d_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}，\forall l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L（17）

\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},q}{x}_{r_{i,j}^k}\le d_q，\forall q\in n_q（18）

\Omega =\left\{x_{r_{i,j}^k}\left|\sum _{r_{i,j}^k\in R}{y}_{r_{i,j}^k,i}{x}_{r_{i,j}^k}=A_i,\sum _{r_{i,j}^k\in R}{y}_{r_{i,j}^k,j}{x}_{r_{i,j}^k}=\lambda B_j\right.,x_{r_{i,j}^k}\ge \mathrm{0}\right\}（19）

上述模型是一个非线性优化模型，针对该模型特性，本研究设计了一个组合式求解算法，该算法包含内外双层迭代算法，外层算法以增广拉格朗日乘子法为基础进行设计，将铁路区间通行能力和港口吞吐量限制约束作为拉格朗日乘子项，构建增广拉格朗日函数。然后，采用可行方向法作为内层迭代算法，将外层算法得到的非线性目标函数线性化。

3.1　外层迭代算法

外层迭代算法采用增广拉格朗日乘子法，将铁路区间通行能力约束和港口吞吐量约束作为拉格朗日乘子项，放入目标函数中，并在拉格朗日函数中增加了一项惩罚函数，从而构建海铁联运增广拉格朗日目标函数。外层迭代算法的求解思路如下。

首先，在公式⒄和公式⒅中，令

g_{\mathrm{1}}\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)=\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}-d_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}，\forall l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L（20）

g_{\mathrm{2}}\left(q\right)=\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},q}{x}_{r_{i,j}^k}-d_q，\forall q\in n_q（21）

再通过添加松弛变量，将不等公式⒄和公式⒅转化为等式形式，得到转化后的等式为

g_{\mathrm{1}}\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)+s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}=\mathrm{0}（22）

g_{\mathrm{2}}\left(q\right)+s_q=\mathrm{0}（23）

式中：s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}和s_q为松弛变量向量。

目标函数的增广拉格朗日函数表达式为

\begin{array}{l}{L}_{\gamma }\left(x,s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},s_q,{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},{\mu }_q\right)=Z_{\mathrm{1}}\left(x\right)-{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}\cdot \left(g_{\mathrm{1}}\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)+s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}\right)+\frac{\mathrm{1}}{\mathrm{2}}\gamma \cdot {\left(g_{\mathrm{1}}\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)+s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}\right)}^{\mathrm{T}}\cdot \left(g_{\mathrm{1}}\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)+s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}\right)-\\ {\mu }_q\cdot \left(g_{\mathrm{2}}\left(q\right)+s_q\right)+\frac{\mathrm{1}}{\mathrm{2}}\gamma \cdot {\left(g_{\mathrm{2}}\left(q\right)+s_q\right)}^{\mathrm{T}}\cdot \left(g_{\mathrm{2}}\left(q\right)+s_q\right)\end{array}（24）

式中：{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}，{\mu }_q为拉格朗日乘子向量，其维数分别与L元素个数和n_q元素个数相等；\gamma为惩罚参数。

从而，[M 1]可以转化为求解以下问题。

[M 2]

\mathrm{m}\mathrm{i}\mathrm{n}\underset{x\in \mathrm{\Omega }}{L_{\gamma }}\left(x,s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},s_q,{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},{\mu }_q\right)（25）

\mathrm{s}.\mathrm{t}.s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},s_q>\mathrm{0}（26）

由于[M 2]包含多个变量，因此可对[M 2]分步求最小化。对给定x，首先求s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}和s_q的最小值，然后再求x的最小值。

通过对s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}和s_q求最小值，可以得到

s_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}=\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0},-\frac{{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}}{\gamma }-g\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)\right]，\forall l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L（27）

s_q=\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0},-\frac{{\mu }_q}{\gamma }-g\left(q\right)\right]，\forall q\in n_q（28）

把公式27和公式28带入目标函数，则可以得到增广拉格朗日函数问题。

［M 3］

\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}\underset{x\in \Omega }{L_{\gamma }}\left(x,{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},{\mu }_q\right)=Z_{\mathrm{1}}\left(x\right)+\frac{\mathrm{1}}{\mathrm{2}\gamma }\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}\left[\mathrm{m}\mathrm{a}\mathrm{x}{\left(\mathrm{0},{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}+\gamma \cdot g_{\mathrm{1}}\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)\right)}^{\mathrm{2}}-{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\mathrm{2}}\right]+\\ \frac{\mathrm{1}}{\mathrm{2}\gamma }\sum _{q\in n_q}\left[\mathrm{m}\mathrm{a}\mathrm{x}{\left(\mathrm{0},{\mu }_q+\gamma \cdot g_{\mathrm{2}}\left(l_q\right)\right)}^{\mathrm{2}}-{\mu }_q^{\mathrm{2}}\right]\end{array}（29）

接下来，假设当前迭代点为{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\left(k\right)}和{\mu }_q^{\left(k\right)}，通过求解[M 3]，得到x^{\left(k\right)}。[M 3]的求解步骤见下一节内层迭代算法。此时惩罚参数和拉格朗日乘子可更新如下。

{\gamma }^{\left(k+\mathrm{1}\right)}\Leftarrow \kappa \cdot {\gamma }^{\left(k\right)}\left(\kappa \ge \mathrm{1}\right)（30）

{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\left(k+\mathrm{1}\right)}\Leftarrow \mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0},{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\left(k\right)}+{\gamma }^{\left(k\right)}\cdot g\left(l_{n_{\mathrm{1}},n_{\mathrm{2}}}\right)\right]（31）

{\mu }_q^{\left(k+\mathrm{1}\right)}\Leftarrow \mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0},{\mu }_q^{\left(k\right)}+{\gamma }^{\left(k\right)}\cdot g\left(q\right)\right]（32）

此处惩罚参数更新时的扩张因子取1<κ≤2，该范围可以避免内层算法在求解子问题时发生病态，无法得到最优解。

3.2　内层迭代算法

内层迭代算法是用以求解外层算法在{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\left(k\right)}和{\mu }_q^{\left(k\right)}条件下的子问题[M 3]，即在海铁联运调运优化问题中每一次外层迭代，都需要求解[M 3]。

[M 3]是一个仅含有等式及非负约束的非线性优化子问题，因此求解该模型可以通过将目标函数线性化，转化为求解线性规划子问题，从而得到目标函数的可行下降方向，这也正是可行方向法的基本思想。

内层迭代算法具体思路如下：假设当前迭代点为x^{\left(\tau \right)}，在该点将[M 3]目标函数线性化，得到以下子问题。

[M4]

\mathrm{m}\mathrm{i}\mathrm{n}\nabla \underset{x\in \Omega }{L_{\gamma }}{\left(x^{\left(\tau \right)},{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\left(k\right)},{\mu }_q^{\left(k\right)}\right)}^{\mathrm{T}}x（33）

假设y^{\left(\tau \right)}为上述子问题的解，则可行下降方向\mathrm{\Delta }{x}^{\left(\tau \right)}为

\mathrm{\Delta }{x}^{\left(\tau \right)}=y^{\left(\tau \right)}-x^{\left(\tau \right)}（34）

沿该方向进行一维搜索，即通过求解

\underset{\mathrm{0}\le \lambda \le \mathrm{1}}{\mathrm{m}\mathrm{i}\mathrm{n}}\underset{x\in \mathrm{\Omega }}{L_{\gamma }}\left(\left(x^{\left(\tau \right)}+\lambda \mathrm{\Delta }{x}^{\left(\tau \right)}\right),{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}^{\left(k\right)},{\mu }_q^{\left(k\right)}\right)（35）

得到最优步长\lambda，则下一个迭代点x^{\left(\tau +\mathrm{1}\right)}为

x^{\left(\tau +\mathrm{1}\right)}=x^{\left(\tau \right)}+\lambda \mathrm{\Delta }{x}^{\left(\tau \right)}（36）

为得到y^{\left(\tau \right)}，需要求解[M4]。[M4]可表示为

\mathrm{m}\mathrm{i}\mathrm{n}{Z}_{\mathrm{2}}\left(x\right)={\left[\frac{\partial L_{\gamma }\left(x_{r_{i,j}^k},{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},{\mu }_q\right)}{\partial x_{r_{i,j}^k}}\right]}^{\mathrm{T}}\cdot x（37）

\mathrm{s}.\mathrm{t}.\sum _{r_{i,j}^k\in R}{y}_{r_{i,j}^k,i}{x}_{r_{i,j}^k}=A_i，\forall i\in n_i（38）

\sum _{r_{i,j}^k\in R}{y}_{r_{i,j}^k,j}{x}_{r_{i,j}^k}=\lambda B_j，\forall j\in n_j（39）

x_{r_{i,j}^k}\ge \mathrm{0}，\forall r_{i,j}^k\in R（40）

其中，令\varpi =\frac{\partial L_{\gamma }\left(x_{r_{i,j}^k},{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},{\mu }_q\right)}{\partial x_{r_{i,j}^k}}，\varpi可视为每个OD对的广义路径阻抗，该阻抗具体推导过程如下。

\begin{array}{l}{\varpi }_{r_{i,j}^k}=\frac{\partial L_{\gamma }\left(x_{r_{i,j}^k},{\mu }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}},{\mu }_q\right)}{\partial x_{r_{i,j}^k}}=\sum _{r_{i,j}^k\in R}{c}_i+\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{p}_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}+\sum _{r_{i,j}^k\in R}{c}_q{y}_{r_{i,j}^k,q}+f\left(\left(\sum _{r_{i,j}^k\in R}\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}{\theta }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}\right)-\alpha \right)-\sum _{r_{i,j}^k\in R}{p}_j+\\ \mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{0},\sum _{l_{n_{\mathrm{1}},n_{\mathrm{2}}}\in L}\left({\beta }_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}+\gamma \cdot \left(\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}{x}_{r_{i,j}^k}-d_{l_{n_{\mathrm{1}},n_{\mathrm{2}}}}\right)\cdot y_{x_{r_{i,j}^k},l_{n_{\mathrm{1}},n_{\mathrm{2}}}}\right)\right)+\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{0},\sum _{q\in n_q}\left({\beta }_q+\gamma \cdot \left(\sum _{r_{i,j}^k\in R}{y}_{x_{r_{i,j}^k},q}{x}_{r_{i,j}^k}-d_q\right)\cdot y_{x_{r_{i,j}^k},q}\right)\right)\end{array}

通过上述推导过程可以看出，通过将[M3]的目标函数线性化，所得到的子问题[M4]，恰巧是一个标准的Hitchcock问题。因此其可运用现有Hitchcock问题求解算法^[20-23]快速求解。

CA流程图如图3所示。

4　案例研究

4.1　案例背景

以国家能源投资集团有限责任公司(以下简称“国家能源集团”)煤炭海铁联运网络为例，采用其煤炭业务相关数据作为参数，构建煤炭能源企业海铁联运调运方案优化模型，内层迭代算法Hitchcock模型采用Gurobi进行求解。国家能源集团海铁联运网络示意图如图4所示。

该网络中参与煤炭生产、装卸、运输的节点一共有188个，其中煤炭装车站47个，电厂33个，港口8个，其余为其他车站。煤炭矿区集中在西北地区，因此装车站也集中分布于西北地区，在铁路运输部分，煤炭由西向东进行运输，在航运部分，煤炭由北向南运输。根据节点的邻接矩阵可以得到网络中从任一生产节点到任一销售区域共有7 134条可行路径。将7 134条可行路径进行编号，编号为No.0~No.7133。初始化最短可行路径集R^{\left(\mathrm{0}\right)}选取85条路径，初始可行路径集中部分可行路径信息如表1所示。

4.2　计算过程及结果

对于阶梯式碳税价格区间的设定，张济建等^[15]对既有研究中针对我国碳税政策发展趋势和合理价格区间的研究进行了总结和探讨。其研究结果表明，将碳税税率设定在29~73元/t之间波动更符合我国的经济发展和碳排放现状。根据《中国能源统计年鉴2018》中能源消耗量折算为碳排放量的结果，铁路运输的单位碳排放按照0.035 kg/(t·km)，航运碳排放按照0.101 kg/(t·km)计算^[24]。R^{\left(\mathrm{0}\right)}条件下进行各OD对的调运量加载，生成初始方案，其总企业效益为-1.610 8×10¹¹元。再进行第二次迭代，煤炭调运方案的总企业效益为1.318 3×10¹⁰元。再进行第三次迭代，生成新调运方案总效益为1.704 4×10¹⁰元。然后，方案满足收敛性判断条件，因此算法终止，1.704 4×10¹⁰元即为该问题的最优解。

为验证设计的CA在求解质量和效率方面具有良好的性能，算例直接在包含7 134条的全可行路径集条件下，采用遗传算法(Genetic Algorithm，GA)对原始非线性规划模型进行求解，得到不同算法结果对比图如图5所示，不同算法的求解时间如表2所示。

由图5可以看出，GA随机生成的初始种群最优解是-2.975×10¹¹元，而CA在R^{\left(\mathrm{0}\right)}条件下生成的初始方案为-1.610 8×10¹¹元，相较于GA提升了45.9%。GA收敛后的最优解为8.765×10⁹元，而CA满足收敛条件后的最优解为1.704 4×10¹⁰元，相较于GA提升了94.5%。在收敛性方面，GA在迭代323次后收敛，而CA仅迭代3次即满足收敛性判断条件，但是由于所设计的CA为启发式算法，GA为元启发式算法，因此两者的收敛次数不具有可比性。

虽然收敛性不具有可比性，但是收敛时间仍能够表明算法的求解效率。由表2可以看出，CA相较于GA在求解效率方面提升了69.8%。

综上所述，所设计的CA针对煤炭能源企业海铁联运调运方案优化模型的求解，不论在初始解的生成、最优解的质量和求解效率方面都明显优于GA，证明CA在求解质量和效率方面均具有更好的性能。

4.3　不同碳排放政策下的敏感性分析

为对比固定碳税和阶梯式碳税对煤炭海铁联运方案企业效益和减排效果的影响，验证阶梯式碳税相较于固定碳税的优势，根据文献[15]对相关碳税研究的归纳总结，在0~200元/t的范围上分别对采用固定碳税政策和阶梯式碳税政策进行敏感性分析。在固定式碳税敏感性分析中，碳税的取值由0元/t向200元/t逐渐增加，在阶梯式碳税中，碳税取值为区间形式，上限由0元/t向200元/t逐渐增加，下限则根据上限取值，遍历所有合理范围内的取值。固定碳税敏感性分析结果如图6所示，阶梯式碳税敏感性分析结果(a)如图7所示，阶梯式碳税敏感性分析结果(b)如图8所示。

在图6中，固定式碳税在由0元/t向200元/t逐步增加时，最优方案的企业效益逐步降低，但最优方案的总碳排放量不变。因为当碳税的税率固定不变时，根据碳排放量的计算方法，总碳排放量与各方案的运输距离和货运量相关。这导致固定碳税在模型目标中表述为一次项，相当于在原决策变量的系数中增加了一项成本系数，不会对最优方案的路径选择产生影响。在目前政策实行过程中，大多使用固定式碳税，由于其计算方式简单，政策实行方便。但是固定式碳税的实施更多地会导致政策“一刀切”的后果，运输企业会在保证货运量不变的情况下涌向最短路径，这虽然能够减少总碳排放量，但是不利于海铁联运网络的整体均衡发展。相较之下，阶梯式碳税的灵活性能够在促进货运企业减少碳排放的情况下，促进海铁联运网络的发展。

图7为阶梯式碳税的上限和下限由0元/t向200元/t逐步增加时的最优方案企业效益，这一敏感性分析指标旨在探究阶梯式碳税区间变化对企业减排积极性的影响。可以看到，当碳税下限小于110、碳税上限大于90时，会出现最优企业效益小于0的情况。这与既有研究得到的不适宜的阶梯式碳税税率和梯度会使企业利益受损，打击企业减排积极性^[15]的结论相吻合。

图8为阶梯式碳税的上限和下限由0元/t向200元/t逐步增加时最优方案的总碳排放量，这一敏感性分析指标旨在探究阶梯式碳税区间变化对减排效果的影响。可以看到，当阶梯式碳税的上限逐步减小，下限逐步升高时，最优方案的总碳排放量在逐渐减少。这说明，在一定程度上减小阶梯式碳税上限和下限之间的差距，能够提升碳税的减排效果。

结合图7和图8可以得到，当阶梯式碳税下限逐渐上升，上限逐渐下降时，该碳税下的最优调运方案的企业效益和减排效果均呈上升趋势，其结果将优于相同条件下采用固定碳税得到的最优方案的企业效益和减排效果。

5　结论

通过对阶梯式碳税下煤炭海铁联运网络调运优化问题的研究，得到结论如下。

（1）煤炭海铁联运网络调运优化模型可以在考虑铁路与港口运输能力约束的基础上，有效地实现能源企业与电厂间的供需关系平衡，并在以企业总效益最大化为目标的基础上，促进和提升调运方案减排效果。

（2）所设计的CA能够满足大规模海铁联运网络的现实场景需求，能够达到较高的求解质量。

（3）阶梯式碳税相较于固定碳税更适用于我国煤炭海铁联运场景，能够在保证货运企业效益的前提下，促进货运企业的碳减排效果，同时保障海铁联运网络货运通道的均衡建设和发展。

（4）国家相关部门在制定碳税政策时，需考虑阶梯式碳税的合理上下限和阶梯高度，以求在保证能源企业生产效益的基础上实现节能减排效果的最大化。根据我国碳排放现状和相关敏感性分析结果，研究建议国家相关部门针对国家能源集团的阶梯式碳税区间定为29~73元/t。

参考文献

[1] 陈维亚，龚浩，方晓平. 考虑运输碳税与质量承诺的多式联运路径优化[J]. 铁道科学与工程学报，2022，19(1)：34-41.

[2] CHEN Weiya， GONG Hao， FANG Xiaoping. Multimodal Transportation Route Optimization Considering Transportation Carbon Tax and Quality Commitment[J]. Journal of Railway Science and Engineering，2022，19(1)：34-41.

[3] 李珺，杨斌，朱小林. 混合不确定条件下绿色多式联运路径优化[J]. 交通运输系统工程与信息，2019，19(4)：13-19，27.

[4] LI Jun， YANG Bin， ZHU Xiaolin. Path Optimization of Green Multimodal Transportation under Mixed Uncertainties[J]. Journal of Transportation Systems Engineering and Information Technology，2019，19(4)： 13-19，27.

[5] 邓学平，陈露，田帅辉. 不确定需求下考虑混合时间窗的多式联运路径优化[J]. 重庆邮电大学学报(自然科学版)，2021，33(4)：689-698.

[6] DENG Xueping， CHEN Lu， TIAN Shuaihui. Multimodal Transport Path Optimization Considering Mixed Time Windows under Uncertain Demand[J]. Journal of Chongqing University of Posts and Telecommunications (Natural Science Edition)，2021，33(4)：689-698.

[7] 梅梦婷，米小珍，郑晓军. 基于DE和NSGA-Ⅱ的集装箱多式联运的路径优化算法[J]. 现代电子技术，2019，42(15)： 144-149.

[8] MEI Mengting， MI Xiaozhen， ZHENG Xiaojun. Path Optimization Algorithm for Container Multimodal Transportation Based on DE And NSGA-‍Ⅱ[J]. Modern Electronics Technique，2019，42(15)：144-149.

[9] 袁旭梅，降亚迪，张旭. 低碳政策下基于区间的模糊多式联运路径鲁棒优化研究[J]. 工业工程与管理，2021，26(4)：134-141.

[10] YUAN Xumei， JIANG Yadi， ZHANG Xu. Research on Robust Optimization of Interval-based Fuzzy Multimodal Transport Paths under Low-Carbon Policies[J]. Industrial Engineering and Management，2021，26(4)：134-141.

[11] ZHANG X， JIN F Y， YUAN X M， et al. Low-Carbon Multimodal Transportation Path Optimization under Dual Uncertainty of Demand and Time[J]. Sustainability，2021，13(15)：8180-8198.

[12] JIANG J H， ZHANG D ZH， MENG Q， et al. Regional Multimodal Logistics Network Design Considering Demand Uncertainty and CO₂ Emission Reduction Target：A System-Optimization Approach[J]. Journal of Cleaner Production，2020(248)：119304-119324.

[13] 蒋琦玮，林艺，冯芬玲. 模糊时间下考虑碳税值变化的多式联运路径优化问题研究[J]. 工业技术经济， 2020， 39(4)： 81-88.

[14] JIANG Qiwei， LIN Yi， FENG Fenling. Research on Multimodal Transport Path Optimization Problem Considering Carbon Tax Value Changes Under Fuzzy Time[J]. Journal of Industrial Technological Economics，2020，39(4)：81-88.

[15] 张旭，袁旭梅，降亚迪. 需求与碳交易价格不确定下多式联运路径优化[J]. 系统工程理论与实践，2021，41(10)：2609-2620.

[16] ZHANG Xu， YUAN Xumei， JIANG Yadi. Optimization of Multimodal Transportation under Uncertain Demand and Stochastic Carbon Trading Price[J]. Systems Engineering-Theory & Practice，2021，41(10)：2609-2620.

[17] 刘松，邵毅明，彭勇. 碳排放限制下的冷藏集装箱多式联运路径优化[J]. 应用数学和力学，2020，41(2)：204-215.

[18] LIU Song， SHAO Yiming， PENG Yong. Optimization of Intermodal Transportation Path of Refrigerated Container under Carbon Emission Limitation[J]. Applied Mathematics and Mechanics，2020，41(2)：204-215.

[19] ZHAO S T， ASHRAF M A. Study on Optimization of Multimodal Transportation of Marine Container Considering Carbon Emission[J]. Desalination and Water Treatment， 2021(219)：84-89.

[20] 翁智雄，程翠云，葛察忠，等. 碳税政策视角下的中国碳减排政策研究[J]. 环境保护科学，2018，44(3)：1-7.

[21] WENG Zhixiong， CHENG Cuiyun， GE Chazhong， et al. Research on China's Carbon Emission Reduction Policy from The Perspective of Carbon Tax Policy[J]. Environmental Protection Science，2018，44(3)：1-7.

[22] 王金南，严刚，姜克隽，等. 应对气候变化的中国碳税政策研究[J]. 中国环境科学，2009，29(1)：101-105.

[23] WANG Jinnan， YAN Gang， JIANG Kejun， et al. The Study on China's Carbon Tax Policy to Mitigate Climate Change[J]. China Environmental Science，2009，29(1)：101-105.

[24] 姚昕，刘希颖. 基于增长视角的中国最优碳税研究[J]. 经济研究，2010，45(11)：48-58.

[25] YAO Xin， LIU Xiying. Optimal Carbon Tax in China with the Perspective of Economic Growth[J]. Economic Research Journal，2010，45(11)：48-58.

[26] 张济建，丁露露，孙立成. 考虑阶梯式碳税与碳交易替代效应的企业碳排放决策研究[J]. 中国人口·资源与环境，2019，29(11)：41-48.

[27] ZHANG Jijian， DING Lulu， SUN Licheng. Research on Enterprise Carbon Emission Decision Considering the Substitution Effect of Stepped Carbon Tax and Carbon Trading[J]. China Population，Resources and Environment，2019，29(11)：41-48.

[28] GUPTA P， MEHLAWAT M K， AGGARWAL U， et al. An Integrated AHP-DEA Multi-Objective Optimization Model for Sustainable Transportation in Mining Industry[J]. Resources Policy， 2021(74)：101180.

[29] LI N， CHEN W Y. Modeling China's Interprovincial Coal Transportation under Low Carbon Transition[J]. Applied Energy，2018(222)：267-279.

[30] LIU F， LV T， SAJID M， et al. Optimization for China's Coal Flow Based on Matching Supply and Demand Sides[J]. Resources， Conservation and Recycling，2018(129)：345-354.

[31] 郭波，陈频. 基于全局运输成本的煤炭铁水联运通道走向分析[J]. 交通运输系统工程与信息，2018，18(S1)：30-34.

[32] GUO Bo， CHEN Pin. Direction Analysis of Railway-waterway Multimodal Transport Corridor Based on the Global Transportation Cost[J]. Journal of Transportation Systems Engineering and Information Technology，2018，18(S1)：30-34.

[33] BRENNER U. A Faster Polynomial Algorithm for the Unbalanced Hitchcock Transportation Problem[J]. Operations Research Letters，2008，36(4)：408-413.

[34] MATSUI T， SCHEIFELE R. A Linear Time Algorithm for the Unbalanced Hitchcock Transportation Problem[J]. Networks，2016，67(2)：170-182.

[35] MORAIS R G， MOURELLE L M， NEDJAH N. Hitchcock Birds Inspired Algorithm[C]//Computational Collective Intelligence：10th International Conference，ICCCI 2018，Proceedings， Part II 10. Bristol：Springer International Publishing，2018：169-180.

[36] MORAIS R G， NEDJAH N， MOURELLE L M. A Novel Metaheuristic Inspired by Hitchcock Birds' Behavior for Efficient Optimization of Large Search Spaces of High Dimensionality[J]. Soft Computing，2020，24(8)：5633-5655.

[37] 周晴. 电力企业电煤多式联运路径优化问题研究[D]. 北京：华北电力大学，2020.

[38] ZHOU Qing. Study on the Optimization of the Multimodal Transport Path of Thermal Coal in Power Enterprises[D]. Beijing：North China Electric Power University，2020.

基金资助

国家重点研发计划项目(2016YFB1200105)

国家重点研发计划项目(2016YFB0900502)

国家电网电力科学研究院研发项目(FD-GF-DWGS-200001)