全球气候变化日益加剧，氢能源由于不含碳且终端使用效率高，被认为是降低温室气体排放的解决方案之一^[1]。工业副产氢是一种成本低廉、无额外碳排放的氢气来源，经过相应技术净化可达到国际标准^[2]。副产氢可用于基于燃料电池的热电联产、燃料电池汽车生产和家庭供暖等方面^[3]，但目前存储和供应链基础设施缺乏，商业模式有待开发。储氢可使氢气产销规模化，并解决产销不匹配的问题。在大规模储氢方面，由于地下盐穴具有密封能力强、材料惰性高和注采能力高等优势，盐穴储氢被认为是最有前景的储氢技术之一^[4]。

氢市场是氢气供应链的核心，包括制造、运输、储存和销售等环节。尽管氢市场尚不成熟，但已展现巨大潜力。液氢和气态氢储运与加注方法近年来被学界广泛研究，文[5-8]对不同氢气供应链发展模式下的成本及成本优化进行了对比和分析。然而，针对产(储)氢的地理位置、制氢厂的产氢量、氢气处理装置容量选择、运输方式和路径选择等问题，仍需展开深入研究。

现阶段盐穴仍不具备储存和销售副产氢的业务能力，由此, 本研究设想了一种由盐穴和化工厂构成的副产氢供应链，研究并优化该供应链中盐穴和化工厂的动态决策过程，从而实现包括盐穴和上游化工厂在内的供应链整体利益最大化。其中，化工厂将剩余副产氢运输至盐穴进行存储，由盐穴将副产氢出售给天然气公司等用户赚取利润。在该供应链中，考虑副产氢市场需求的同时，如何将多家化工厂的副产氢以最低成本进行二次处理并运输至盐穴进行存储，构成了副产氢集散的关键环节。本研究将针对副产氢的收储环节，重点研究由多家化工厂和盐穴构成的局部供应链。具体而言，提出了一种盐穴收储化工厂副产氢、化工厂联合运输以降低运输成本的商业模式；建模并优化了该商业模式下副产氢供应链中盐穴和化工厂的决策过程。

本研究的目的是探索所提商业模式可能存在的最大利润空间，所建立的优化模型以实现经济利益最大化为目标，且化工厂和盐穴都参与决策。进一步探究所提商业模式的利润空间发现，当利润空间较小时，商业模式可能不存在，因为这会降低最大经济利益的预期并进一步降低利润空间；如果利润空间较大，未来可通过进一步设计副产氢市场创造环境，鼓励各方合作以实现总经济利益最大化。

1 商业模式 1.1 商品流和收益流

本研究提出的商业模式由多家化工厂、盐穴和用户组成，相应的副产氢供应链流程如图 1所示。其中，化工厂是副产氢的供应商，将副产氢运输至盐穴进行储存；盐穴是副产氢的供应商，将副产氢出售给终端用户；用户是副产氢的消费者，从盐穴收储副产氢。

本研究所提商业模式的基本假定如下：

(1) 对于终端用户的纯度要求，各化工厂需要根据本厂的副产氢纯度和压缩机(液化机)容量确定纯化仪器的容量。

(2) 考虑到各化工厂的副产氢存储量有限，对于副产氢的运输方式，本文只考虑长管拖车(气氢)和液罐车(液氢)运输2种陆上运输方式，不考虑地下管道输送。

(3) 化工厂之间或化工厂与盐穴之间的合作关系不变，即各化工厂的运输目的地不随时间发生变化，但各时段交易量可能随氢气处理成本变化而变化。

(4) 本研究主要针对未来氢气市场进行供应链优化决策设计。根据文[9]，未来氢气市场中氢气价格由零售商(本研究中为盐穴)决定，用户需求可通过与零售价格呈负相关的线性需求函数描述。

1.2 商业模式存在的问题及解决方案

在1.1节提出的盐穴大规模储氢商业模式下，各个副产氢供应商(上游化工厂)和零售商(盐穴)构成副产氢供应链。该供应链中，各参与方根据市场需求，通过合作协调各方决策以实现共同收益。上述商业模式形成的副产氢供应链与现有氢气供应链存在差异，该差异主要体现在以下两个方面：

1) 主要成本构成不同。

与现有氢气供应链中供应商和零售商需承担制氢、运氢和储氢成本不同，作为化工厂主产品生产过程中的副产物，副产氢的制氢成本可忽略，副产氢供应链需支付的成本主要来自二次处理成本和运输成本。二次处理包括纯化和压缩(液化)2个环节，其成本主要包括成本Ⅰ和成本Ⅱ，成本Ⅰ为氢气处理设备(纯化仪器、压缩机/液化机)的初始投资成本，成本Ⅱ为氢气处理设备的运行成本，本研究主要考虑购电成本。运输成本主要包括成本Ⅲ和成本Ⅳ，成本Ⅲ为氢气运输设备(长管拖车和液罐车)的初始投资成本，成本Ⅳ为氢气运输设备的运行成本。

在建立副产氢供应链的数学模型时，针对成本Ⅰ，应考虑各化工厂副产氢的生产规模及各时段氢气二次处理量，以便合理规划纯化仪器、压缩机/液化机的容量；针对成本Ⅱ，由于电价随时段变化而波动，因此各氢气处理设备的运行策略应按小时进行建模；成本Ⅲ和Ⅳ主要由运输方式、运输量和运输距离决定。该数学模型应当能够优化各化工厂的运输方式，从而降低初始投资成本和运行成本。

2) 规划和调度的可协调性不同。

在现有的氢气供应链研究中，氢气的供应商和零售商在规划阶段已对制氢规模和技术与运输距离和方式进行了合理的协调。文[5-8]根据预测年的市场需求给定生产规模，进一步估计生产和运输设备的投资成本，并通过假定的某个固定电价粗略估计氢气生产的年运行成本。

然而，副产氢供应链中化工厂和盐穴的相对位置固定，且副产氢的生产规模易受钢铁、烧碱、氯等化工厂主产品生产计划的限制，进而限制了氢供应链的规划和调度灵活性，这将带来以下问题：

(1) 对于距离盐穴较远、日副产氢量在特定范围的化工厂，若采用气氢方式运输，则大量长管拖车远距离运输会产生较大的运行成本。

(2) 若采用液氢方式运输，则液罐车加注时间过长会造成较大的挥发损失。

因此，本研究设计了一个多家化工厂联合运输的场景。为降低副产氢远距离运输至盐穴的成本，在多家相距较近的化工厂中选择一家作为转运中心，先采用气氢方式集中运输至转运中心，再由转运中心进行液化并远距离运输至盐穴。图 2以3家化工厂和一家盐穴所构成的副产氢供应链为例，展示了2种可能的运输路径。图 2中，d_{x, y}(x, y=1, 2, 3, 4)表示化工厂x和y之间的转运关系，当化工厂x以化工厂y作为运输目的地时，d_{x, y}=1，否则d_{x, y}=0。蓝线表示3家化工厂以化工厂2为运输中心，化工厂1和3将副产氢以气氢方式运输至化工厂2，再由化工厂2将3家化工厂的副产氢进行液化并统一运输至盐穴进行储存；紫线表示化工厂2将副产氢以气氢方式运输至化工厂1，再由化工厂1将两家化工厂的副产氢进行液化并运输至盐穴，而化工厂3则单独将副产氢以气氢方式运输至盐穴。

此外，不同于文[5-8]中以季或年为时间尺度进行供应链优化决策建模，为了详细计及副产氢二次处理成本和运输成本，需要选择一个更详细的时间尺度，即在该时间尺度上不仅可以计及不同化工厂与盐穴运输距离间的差距，还可以计及电价波动和液氢挥发损失，从而为多家化工厂联合运输方式规划及氢气处理量调度奠定基础。

2 数学模型 2.1 问题建模

首先，对副产氢的收益流进行建模。盐穴和化工厂根据市场的副产氢需求出售氢气，即

$m_t^{\mathrm{dis}}=D\left(\rho_t\right), $

(1)

$\rho_t^{\min } \leqslant \rho_t \leqslant \rho_t^{\max }, \quad \forall t.$

(2)

其中：t为全优化过程的时段，ρ_t为t时段盐穴单位质量氢气的收购价，D为副产氢需求函数，ρ^min、ρ^max分别为盐穴从化工厂收购单位质量氢气的价格上下限，m_t^dis表示盐穴和化工厂根据市场的副产氢需求出售的氢气量。

根据文[9]，市场的副产氢需求可通过线性需求函数描述，具体可表示为

$D\left(\rho_t\right)=D_0-\alpha \rho_t.$

(3)

其中：D₀为ρ_t=0时市场副产氢的需求量，α表示副产氢需求随副产氢价格变化的斜率。

根据以上各时段市场的单位售氢价格和市场需求，上游供应链的售氢收入为$\sum\limits_{t=1}^T \rho_t m_t^{\text {dis }}$，T为时段总数。

1) 盐穴荷气状态变化。

假设运氢车(长管拖车或液罐车)到达盐穴后直接将氢气充入盐穴，不会在盐穴停留至下一时间段，则t时段内盐穴充气量m_t^ch和T_tran^{i, j}时段前从化工厂i运出的氢气量$m_{i, t-T_{\text {tran }}^{i, j}}^{\operatorname{tran}}$满足如下关系式：

$m_t^{\mathrm{ch}}=\sum\limits_{i=1}^I d_{i, I+1} m_{i, t-T_{\operatorname{tran}}^{i, j}}^{\operatorname{tran}}.$

(4)

其中：i、j为化工厂索引编号，i, j=1, 2, …, I，I为化工厂总数，I+1为盐穴的索引编号；d_{i, I+1}为0/1变量，表示化工厂i是否以盐穴为运输目的地；m_{i, t}^tran为t时段末从i运出的氢气总质量; T_tran^{i, j}表示从i到j的运输时长。

盐穴储气状态(state of charge，SOC)满足如下关系：

$A_t^{\mathrm{SOC}}+1=A_t^{\mathrm{SOC}}+\frac{1}{M_{\max }}\left(m_t^{\mathrm{ch}}-m_t^{\mathrm{dis}}\right).$

(5)

其中：M_max为盐穴储气最大质量，m_t^dis为t时段内盐穴出气量，A_t^SOC为t时段盐穴SOC状态。

2) 氢气加工及运输的动态过程。

随着化工厂生产的实时进行，氢气作为副产品实时产生并纯化后存在3种去向：

(1) 通过压缩机/液化机处理后直接充装至长管拖车/液罐车。在t时段末，仅有充装质量达到单车容量上限的车辆方可驶离，而未装满的车辆仍停留在化工厂，直至充装至单车容量上限。因此，i的氢气交易量实际由氢气加工质量和单车单程最大运输质量决定。

对于采用气氢方式运输的化工厂，其氢气运输的动态过程可表示如下：

$n_{i, t}^{\mathrm{cars}} \leqslant \frac{m_{i, t}^{\mathrm{pr}}+m_{i, t-1}^{\mathrm{st}}}{M_{\mathrm{com}}^{\mathrm{car}}} \leqslant n_{i, t}^{\mathrm{cars}}+1, \quad \forall t ;$

(6)

$m_{i, t}^{\text {tran }}=n_{i, t}^{\text {cars }} M_{\text {com }}^{\text {car }}, \quad \forall t ;$

(7)

$\begin{gathered}m_{i, t}^{\mathrm{st}}=m_{i, t-1}^{\mathrm{st}}+m_{i, t}^{\mathrm{pr}}-m_{i, t}^{\mathrm{tran}}, \\ \forall t \in\{2, 3, \cdots, T\} .\end{gathered}$

(8)

其中：M_com^car为单辆长管拖车的单程最大运输质量，kg；n_{i, t}^cars为t时段离开i的车辆数量; m_{i, t}^pr为t时段内i加工的氢气质量；m_{i, t}^st为t时段末停留在i的未装满车辆的车内氢气质量。式(6)和(7)表示t时段末已加注至容量上限的长管拖车从化工厂运出氢气，因此从i运出的氢气质量应为长管拖车最大运输质量的整数倍。式(8)表示各时段末因未装满停留在化工厂的长管拖车中气氢的变化情况。

对于采用液氢运输的化工厂，氢气运输的动态过程可表示如下：

$n_{i, t}^{\mathrm{cars}} \leqslant \frac{m_{i, t}^{\mathrm{pr}}+\gamma_{\mathrm{L} 1} m_{i, t-1}^{\mathrm{st}}}{M_{\mathrm{liq}}^{\mathrm{clar}}} \leqslant n_{i, t}^{\mathrm{cars}}+1, \quad \forall t ;$

(9)

$m_{i, t}^{\mathrm{tran}}=n_{i, t}^{\mathrm{cars}} M_{\mathrm{liq}}^{\mathrm{car}} \sum\limits_{j=1}^{I+1} d_{i, j} \gamma_{\mathrm{L} 2}^{i, j}, \quad \forall t ;$

(10)

$\begin{gathered}m_{i, t}^{\mathrm{st}}=\gamma_{\mathrm{L} 1} m_{i, t-1}^{\mathrm{st}}+m_{i, t}^{\mathrm{pr}}-m_{i, t}^{\mathrm{tran}} / \sum\limits_{j=1}^{I+1} d_{i, j} \gamma_{\mathrm{L} 2}^{i, j}, \\ \forall t \in\{2, 3, \cdots, T\} .\end{gathered}$

(11)

其中：γ_L2^{i, j}表示液氢运输过程的单程液氢挥发剩余比例，γ_L2^{i, j}=γ_L2T_tran^{i, j}，γ_L2为液氢运输过程的单时段液氢挥发剩余比例; γ_L1为液化装车过程的单时段液氢挥发剩余比例；M_liq^car为单辆液罐车的单程最大运输质量；d_{i, j}为0/1变量，表示i是否以j为运输目的地。式(9)和(10)表示t时段末从i运出的氢气质量是液罐车最大运输质量的整数倍。式(11)表示各时段末因未装满停留在化工厂的液罐车中液氢的变化情况。

(2) 未经处理的副产氢可在本地利用低压罐缓存，且将在未来某一时段进行处理并运输。

(3) 直接废弃。如果用户某一时段的收购价过低，而本地低压罐容量有限，化工厂可能会选择将副产氢进行废弃。

由副产氢去向(2)和(3)得到各化工厂本地低压罐中的氢气量变化情况，具体表示如下：

$\begin{gathered}m_{i, t}^{\mathrm{unpr}} \leqslant m_{i, t-1}^{\mathrm{unpr}}+M_{i, t}-m_{i, t}^{\mathrm{pr}}+ \\ \sum\limits_{j=1}^I d_{j, i} m_{i, t-T_{\text {tran }}^{i, j}}^{\mathrm{tran}}, \\ \forall t \in\left\{\max \left(1, T_{\text {tran }}^{i, j}\right), \max \left(1, T_d^{i, j}\right)+1, \cdots, T\right\} .\end{gathered}$

(12)

其中：M_{i, t}为t时段内i的副产氢量，m_{i, t}^unpr为t时段末由于未经适当处理(足够的压缩比/液化)暂存在低压罐中的氢气质量。式(12)表明i可以选择性地将副产氢就地存储或废弃。

2.2 供应链决策优化

对副产氢供应链中盐穴和化工厂的决策问题进行建模，从而实现副产氢供应链的优化。该优化问题包含如下约束：

1) 盐穴运行约束。

考虑盐穴单时段进/出气管道输送容量上限和稳定运行下限，时段t内盐穴进气量m_t^ch和出气量m_t^dis需满足

$M_{\mathrm{ch}}^{\min } \leqslant m_t^{\mathrm{ch}} \leqslant M_{\mathrm{ch}}^{\max }, $

(13)

$M_{\mathrm{dis}}^{\min } \leqslant m_t^{\mathrm{dis}} \leqslant M_{\mathrm{dis}}^{\max }.$

(14)

其中：M_ch^min、M_ch^max分别为盐穴单时段进气量上下限，M_dis^min、M_dis^max分别为盐穴单时段出气量上下限。

此外，盐穴储气SOC约束为

$A_{\min }^{\mathrm{SOC}} \leqslant A_t^{\mathrm{SOC}} \leqslant A_{\max }^{\mathrm{SOC}}, $

(15)

$A_{t=1}^{\mathrm{SOC}}=A_{t=T}^{\mathrm{SOC}}.$

(16)

其中A_min^SOC、A_max^SOC分别为盐穴SOC上下限。

式(13)、(15)和(16)为含整数变量的非线性约束，理由如下：由于d_{i, I+1}(0/1变量)和$m_{i, t-T_{\text {tran }}^{i, j}}^{\mathrm{tr}}$为决策变量，式(4)中m_t^ch包含决策变量的相乘项$d_{i, I+1} m_{i, t-T_{\text {tran }}^{i, j}}^{\operatorname{tran}}$，因此m_t^ch为非线性项，式(13)为m_t^ch的上下限约束，含整数变量的非线性约束；根据式(5)，A_t^SOC的取值由m_t^ch决定，而式(15)为A_t^SOC的上下限约束，式(16)表示A_t^SOC在开始时段与结束时段取值相同，因此式(15)为含整数变量的非线性约束。

2) 化工厂规划及运行约束。

化工厂规划主要包含化工厂运输路径规划、氢气处理装置投资规划，化工厂运行主要为受规划结果约束的氢气加工和运输的动态过程。

长管拖车和液罐车的实际装车过程将为化工厂的规划和调度决策提供多种选择。比如，对于产量为100 kg/h的化工厂，可选择购买容量为100 kg/h的氢气处理装置对副产氢实时加工，也可选择购买容量为1 000 kg/h的氢气处理装置，利用低压罐对10 h内的氢气进行缓存，并在第10 h末一次性加工。上游供应链的决策问题是在多种规划和调度方式中寻找最优决策，从而降低副产氢二次处理(压缩或液化)和运输的成本。

化工厂规划及运行的决策变量包括氢气处理装置的类型、各时段的氢气加工质量及化工厂间的转运关系。

(1) 运输目的地约束。i的运输目的地约束如下：

$\sum\limits_{i=1}^{I+1} d_{i, j}=1, \quad \forall i ;$

(17)

$d_{i, j}+d_{j, i} \leqslant 1, \quad \forall i, j.$

(18)

式(17)表明各化工厂运输目的地唯一；式(18)表明两个化工厂之间不会同时以对方作为运输目的地。

(2) 氢气处理装置选型约束。化工厂i的氢气处理装置选型约束可表示如下：

$\sum\limits_{\xi=1}^{E_{\text {com }}+E_{\text {liq }}} x_i^{\xi}=1.$

(19)

其中：ξ为氢气处理设备的索引编号；E_com为可选的压缩机类型数量，E_liq为可选的液化机类型数量(不同类型仅考虑容量不同)；x_i为由i的压缩机和液化机选型的0/1变量所构成的矩阵，x_i={x_i^ξ}, ∀ξ。式(19)表示各化工厂选取氢气处理装置唯一。

3) 氢气运输量约束。

对于采用气氢方式运输的化工厂，以式(6)—(8)所示的氢气运输动态过程作为约束。对于采用液氢方式运输的化工厂，以式(9)—(11)所示的氢气运输动态过程作为约束，由于n_{i, t}^cars(整数变量)、d_{i, j}(0/1变量)和m_{i, t}^tran均为决策变量，式(10)存在决策变量的相乘项n_{i, t}^carsd_{i, j}，式(11)存在决策变量的相除项$m_{i, t}^{\operatorname{tran}} / \sum\limits_{j=1}^{I+1} d_{i, { }_j} \gamma_{\mathrm{L} 2}^{i, }$，因此式(9)—(11)成为含整数变量的非线性约束。

为计及i本地低压罐中的氢气量变化情况，式(12)也作为氢气运输量的约束。由于d_{j, i}和$m_{i, t-T_{\operatorname{tran}}^{i, j}}^{\operatorname{tran}}$, j_tran^tran均为决策变量，且d_{j, i}为0/1变量，式(12)因存在决策变量的相乘项$d_{j, i} m_{i, t-T_{\mathrm{tran}}^{i, j}}^{\operatorname{tran}}$而成为含整数变量的非线性约束。

记$t_{\mathrm{dr}}^i=\sum\limits_{j=1}^{I+1} d_{i, j} T_{\mathrm{tran}}^{i, j}$，表示i到达目的地的时间。与氢气运输状态无关的约束表示如下：

$\sum\limits_{t=2 t_{\mathrm{dr}}^i}^t n_{i, t}^{\mathrm{cars}} \leqslant N_i^{\mathrm{cars}}, \quad \forall t \in\left\{2 t_{\mathrm{dr}}^i, 2 t_{\mathrm{dr}}^i+1, \cdots, T\right\} ; $

(20)

$m_{i, t}^{\mathrm{pr}} \leqslant \sum\limits_{\xi=1}^{E_{\text {com }}+E_{\mathrm{liq}}} x_i^{\xi} M_{\mathrm{pr}}^{\xi}, \quad \forall t ;$

(21)

$m_{i, t}^{\text {unpr }} \leqslant \sum\limits_{\xi=1}^{E_{\text {com }}+E_{\text {liq }}} x_i^{\xi} M_{\mathrm{pr}}^{\hat{\xi}}, \quad \forall t.$

(22)

其中：M_pr^ξ表示第ξ台氢气处理装置的容量，kg/h；N_i^cars为i的初始购置车辆数。式(20)表示考虑单车往返时间时，i的在途车辆总数不应超过初始购置车辆数。式(21)和(22)分别表示各时段氢气加工质量、本地缓存总量不超过氢气处理装置的容量。

本研究以上游供应链的净收入最大化为副产氢供应链的目标。上游供应链的售氢收入为$\sum\limits_{t=1}^T \rho_t m_t^{\mathrm{dis}}$。上游供应链的成本包括氢气处理设备投资建设费用C_INV1ⁱ、运氢车投资建设费用C_INV2ⁱ、氢气二次处理单时段运行成本C_opeⁱ和氢气二次处理单时段运输成本C_tranⁱ。以$x_i^{\mathrm{com}}=\sum\limits_{\xi=1}^{E_{\mathrm{com}}} x_i^{\xi}$和$x_i^{\operatorname{liq}}=\sum\limits_{\xi=E_{\text {com }}}^{E_{\text {com }}+E_{\text {liq }}} x_i^{\xi}$分别表示化工厂采用气氢运输方式和液氢运输方式的0/1变量，以上4种成本分别表示如下：

$C_{\mathrm{ope}}^i=m_{i, t}^{\mathrm{pr}} \sum\limits_{\xi=1}^{E_{\text {com }}+E_{\mathrm{liq}}} x_i^{\xi} w_t\left(\beta_{\mathrm{p}}+\beta_{\mathrm{com}} x_i^{\mathrm{com}}+\beta_{\mathrm{liq}} x_i^{\mathrm{liq}}\right), $

(23)

$C_{\text {tran }}^i=n_{i, t}^{\text {car }} \sum\limits_{j=1}^{I+1} d_{i, j} K_3 T_{\text {tran }}^{i, j}, $

(24)

$C_{\mathrm{INV} 1}^i=\sum\limits_{\xi=1}^{E_{\text {com }}+E_{\text {liq }}} x_i^{\xi} K_1^{\xi}, $

(25)

$C_{\mathrm{INV} 2}^i=N_i^{\mathrm{cars}}\left(x_i^{\mathrm{com}} K_2^{\mathrm{com}}+x_i^{\mathrm{liq}} K_2^{\mathrm{liq}}\right).$

(26)

其中：w_t为t时段的电价；β_p、β_com和β_liq分别为纯化仪器、压缩机和液化机的单位能耗，kW·h·kg^-1；K₁^ε为氢气处理装置ξ的投资建设费用, K₁={K₁^ξ}, ∀ξ；K₂^com、K₂^liq分别为单辆长管拖车、液罐车的投资建设费用；K₃为单车单时段运行成本，长管拖车和液罐车取值相同。

综合售氢收入和式(23)—(26)的成本，副产氢供应链决策问题的优化目标可表示为

$\begin{gathered}\max \theta=\sum\limits_i^I \sum\limits_{t=1}^T\left(\rho_t m_{i, t}^{\mathrm{tran}} d_{i, I+1}-C_{\mathrm{ope}}^i-\right. \\ \left.C_{\mathrm{tran}}^i\right)-C_{\mathrm{INV} 1}^i-C_{\mathrm{INV} 2}^i .\end{gathered}$

(27)

此外，由于m_{i, t}^pr、x_i^ξ、x_i^com、x_i^liq和N_i^cars均为决策变量，根据式(23)，C_opeⁱ包含m_{i, t}^pr、x_i^ξ、x_i^com(或x_i^liq)的相乘项，其中x_i^ξ、x_i^com和x_i^liq均为整数变量，使优化问题的目标存在含整数变量的非线性项。同样，根据式(26)，C_INV2ⁱ存在N_i^cars、x_i^com(或x_i^liq)的相乘项，而N_i^cars、x_i^com和x_i^liq均为整数变量，使优化问题的目标存在含整数变量的非线性项。

至此，上游供应链的决策问题被建模成一个混合整数非线性规划(mixed integer nonlinear programming，MINLP)问题。该问题中，为使上游供应链的利益最大化，盐穴和化工厂一方面需要确定氢气价格以提高其售氢收益，另一方面需要在副产氢量受化工厂主要产品生产计划限制的情况下，尽可能降低副产氢二次处理(压缩或液化)和运输的成本。

2.3 模型求解

本研究借助遗传算法求解前文构建的MINLP问题。其中，将整数变量矩阵x_i和d_i作为遗传算法的优化变量，x_i={x_iⁿ}, ∀n, d_i={d_{i, j}}, ∀j; 然后针对每个个体，求解以式(27)为目标函数、式(6)—(22)为约束的MINLP问题，并进行个体适应度评估。该过程具体如下。

第1步  初始化。设置进化代数计数器为0，设置最大进化代数为50，随机生成30个个体(x_i, d_i), ∀i∈I作为初始群体P(0)。

第2步  个体评价。给定(x_i, d_i)时，调用MINLP求解器，根据(x_i, d_i)和MINLP所得最优解计算θ，并将其作为群体P(t)中各个个体的适应度。

第3步  群体进化及终止。群体P(t)经过选择、交叉和变异运算之后得到下一代群体P(t+1)。终止条件为达到最大进化代数50，或两代群体适应度变化小于阈值10^-10。

上述方法实现了MINLP问题的求解。相比于完全利用遗传算法，本文方法降低了遗传算法中的个体计算复杂度，并利用线性规划商业求解器保证了下层解的全局最优性，以及问题的求解效率。

3 实例分析

考虑到化工厂数量主要影响副产氢运输路径的规划结果，而3家化工厂可以体现完全独立运输、部分联合运输和完全联合运输这3种运输路径，因此本文设计了一个由1个盐穴和3个化工厂构成的副产氢供应链，并对其进行优化设计。仿真计算在Intel Xeon W-2255处理器、3.70 GHz主频、128 GB内存的电脑上进行。MINLP问题的商业求解器为CPLEX 12.9。

3.1 系统参数

考虑由1个盐穴和3家化工厂构成的副产氢供应链存在转运决策和氢气定价问题, 故设定时间尺度为1 d(分为12个时段)。副产氢供应链场景参数如表 1所示。

气氢和液氢运输相关的经济技术参数参考文[9-10]，氢气净化技术有关成本参考文[4]，氢气运输和净化的相关经济技术参数如表 2所示。其中，单车单时段运行成本综合考虑汽油费、工作人员薪资和车辆维护费用。由于本文设定的时间尺度为1 d，在氢气的加工和运输成本中，各类设备的一次性投资建设费用按一定折扣率折算为日投资费用。

3.2 结果分析

1) 优化结果。

采用Python建模并调用CPLEX对模型进行求解，得到的化工厂运输路径策略、运输方式和氢气处理装置容量规划结果如图 3所示。其中，化工厂1采用容量为2 000 kg/h的压缩机和气氢运输方式，将副产氢压缩后运输至化工厂2。化工厂2采用容量为8 000 kg/h的液化机对化工厂1和2的副产氢进行液化，然后运输至盐穴。化工厂3采用容量为8 000 kg/h的液化机对副产氢进行液化并直接运输至盐穴。

盐穴-化工厂(上游供应链)氢气定价及与终端用户间的交易量如图 4所示。当市场副产氢需求由式(3)描述时，氢气的交易价格为14.17美元/kg，并在1 d的结束时段达到最低价格10.00美元/kg。此时，上游供应链的日收益为253 620.00美元/d，售氢量为33 988.79 kg，占日副产总氢量(65 000 kg)的52.3%。

3家化工厂的总交易量与分时电价如图 5所示。由图 5可知，化工厂倾向于在电价较低时进行氢气加工，以降低氢气加工成本。

以上结果表明，在表 2所示的相关经济技术参数下，该商业模式能够使所设计的上游供应链获得经济利益，从而证明本文所提商业模式成立。

此外，根据图 4和图 5的结果可知，利用所提优化模型进行供应链决策，化工厂的氢气加工量能够响应电价，电价波动最终反映氢气的定价与市场需求量。

2) 供应链管理决策优化带来的经济效益。

文[5-8]中，通常固定生产商(化工厂)转运路线(如生产商直接将氢气运输至存储点)，或固定运行决策(如不考虑氢气市场需求及电价波动，对各化工厂的所有副产氢进行最大限度加工处理)，因此本文将通过固定不同决策变量，展示第2章中将氢气处理方式、转运决策、初始购置车辆数、各时段氢气加工量和运出量均作为决策变量的必要性。由于将上述变量作为决策变量将不可避免地使供应链优化决策问题变成MINLP问题，因此效益提升结果也能够体现将供应链优化决策问题建模成MINLP问题的必要性。

固定不同决策变量的优化结果如图 6所示。其中，固定转运决策表示3家化工厂均采用单独运输至盐穴的方式(d_{1, 4}=d_{2, 4}=d_{3, 4}=1)；固定氢处理装置容量决策1表示采用液化运输方式的化工厂均使用容量为8 000 kg/h的液化机，采用气氢运输方式的化工厂均使用容量为2 000 kg/h的压缩机；固定氢处理装置容量决策2表示采用液化运输方式的化工厂均使用容量为4 000 kg/h的液化机，采用气氢运输方式的化工厂均使用容量为1 200 kg/h的压缩机；固定运行决策表示对各化工厂的所有副产氢进行最大限度加工处理，即各时段副产氢处理量达副产氢生产量上限或氢处理装置容量上限。

由图 6可知，相比于化工厂均进行独立运输，本文所提联合运输方式可使盐穴-化工厂收益提升7.77%。此外，在市场需求一定时，将副产氢进行最大限度加工处理并运输至盐穴进行储存，这并不能使上游供应链利益最大化，仍需要通过合理的调度实现利益最大化，而合理的调度可使利润增加146.6%。

综上所述，本文案例中化工厂间的联合运输、适当的氢处理装置容量选择，以及对氢气处理量进行优化，可使副产氢供应链收益提升，这说明本文所提优化模型在规划阶段优化运输路径、运输方式和氢处理装置容量，以及在运行阶段根据电价波动优化氢处理量十分必要。

3) 灵敏度分析。

为进一步证明本文方法具有一般性和合理性，对单车单时段运行成本、盐穴单时段进/出气管道容量输送上限、市场副产氢线性需求函数参数等进行了敏感度分析，研究其对优化结果的影响。

考虑到启发式算法结果的偏差可能超过参数不同取值得到的结果变化量，本文对整数变量矩阵x_i和d_i(i=1, 2, …, I)的所有可能取值进行遍历，使优化问题转化为MINLP问题，并采用CPLEX得到全局最优解，以获得最优的x_i、d_i和最优x_i、d_i下的其他最优决策。

单车单时段运行成本影响化工厂间的运输目的地选择。不同单车单时段运行成本下的优化结果如图 7所示。蓝色虚线和橙色虚线分别表示考虑化工厂联合运输的优化结果和各家化工厂均单独运输至盐穴的优化结果，紫色实线为化工厂优化运输策略后的结果。分析发现，随着单车单时段运行成本降低，盐穴-化工厂的收益升高。其中，单车单时段运行成本高于310美元/h时，化工厂间进行联合运输；单车单时段运行成本低于310美元/h时，化工厂采用气氢运输方式导致的长管拖车远距离运输成本升高，但低于液氢运输造成的挥发损失，此时产氢量较少的化工厂无需进行联合运输，以便降低成本。

不同盐穴单时段进/出气管道输送容量上限和市场需求下的优化结果如图 8所示。由图 8可知，在管道输送容量上限相同时，盐穴-化工厂收益随着市场需求提高而升高；在相同的市场需求下，盐穴-化工厂收益随着管道输送容量上限提升却单调不减。此外，在某一给定的市场需求下，存在一个最佳管道输送容量，在该最佳值以下，盐穴-化工厂收益随输送容量增加而增加；达到该最佳值后，由于市场需求一定，增加管道输送容量将不会继续提高收益，因此盐穴应根据市场需求规划输送容量。

4) 多场景分析。

(1) 不同单车单时段运行成本下化工厂生产规模和运输距离对优化结果的影响。

为得到副产氢供应链规划及调度的共性指导，本文在由1个盐穴和3个化工厂构成的副产氢供应链的基础场景上，批量修改化工厂副产氢的生产规模、化工厂至盐穴的距离、单车单时段运行成本，以产生更多场景，并求解所提优化模型，从而获得不同产量和运输距离下化工厂间运输路径及方式的一般性结论。

涉及的场景参数在表 2的副产氢供应链场景参数基础上进行修改。具体而言，记产氢规模倍数为ε，运输距离倍数为η，场景集合$\varPsi=\left\{\psi_{\varepsilon}, \eta, K_3\varepsilon \in \mathcal{E}, \eta \in \mathcal{H}, K_3 \in \mathcal{T}\right\}, $。该集合包括$\operatorname{card}(\mathcal{E}) \times \operatorname{card}(\mathcal{H}) \times \operatorname{card}(\mathcal{T}), $个场景，在场景ψ_{ε, η, K3}中，i的生产规模M_{i, t}^ε=εM_{i, t}，化工厂的运输距离矩阵T_tran^η=ηT_tran，单车单时段运行成本为K₃。

当中包括0.5、0.6、0.7、0.8、0.9、1.0、1.2、1.4、1.6、1.8，中包括1、2、3、4、5、6、7、8、9、10，中包括200、300、400、500、600时，将产生500个不同单车单时段运行成本、化工厂产量和运输距离场景。针对上述500个场景，分别求解所提优化模型，得到3家化工厂的最优运输路径、运输方式及联合运输带来的收益，结果如图 9所示。其中，联合运输带来的收益为最优运输路径下的收益与完全独立运输时的收益之差。

在图 9a—9e中，字母空缺处表示3家化工厂完全独立运输，B、C、D、E分别代表一种联合运输路径，其中，B表示化工厂1、2联合运输，化工厂3独立运输；C表示化工厂2、3联合运输，化工厂1独立运输；D表示化工厂1、3联合运输，化工厂2独立运输；E表示化工厂1、2、3联合运输。3种颜色分别代表以盐穴为运输目的地的化工厂采用方式，其中，绿色表示气氢运输，蓝色表示液氢运输，紫色表示3家化工厂中同时存在以气氢和液氢方式运输至盐穴的化工厂。在图 9c—9e中，由斜实线覆盖的场景表示化工厂1的氢气处理量为0，即在该场景下化工厂1处理并运输副产氢的成本将高于出售该部分副产氢所带来的收入，化工厂1将不再处理并运输副产氢；同样，由圆点覆盖的场景表示化工厂1和2每一时刻的氢气处理量均为0，即在该场景下，仅有化工厂3处理并运输氢气；由矩形覆盖的场景表示3家化工厂每一时刻的氢气处理量均为0。

由图 9a—9e可知，当3家化工厂产氢规模不变(ε不变)，运输至盐穴的距离增加(η从1增加至10)时，化工厂将由完全独立运输模式转变为联合运输模式，以降低运输成本；当运输至盐穴的距离不变(η不变)，化工厂产氢规模增加(ε从0.5增加至1.8时)时，采取的联合运输路径可能发生变化，且在大多数场景中，B为最优运输路径。由此可知，当不同化工厂间副产氢生产规模存在一定差距时，中小产量的化工厂在多数场景下更倾向于联合运输，而高产量的化工厂倾向于独立运输。此外，由图 9d和9e可知，当单车单时段运行成本较高(≥500美元/h)时，随着生产规模扩大，运输距离增加，以盐穴为运输目的地的化工厂从气氢运输逐渐转变为液氢运输。由图 9a和9c可知，当单车单时段运行成本较低(＜500美元/h)时，在大多数产量-运输距离场景中，随着生产规模扩大，运输距离增加，以盐穴为运输目的地的化工厂从气氢运输逐渐转变为液氢运输，少量场景(K₃=200, ε=0.7, η=10时)不满足上述规律。这说明，随着单车单时段运行成本降低，气氢运输由于无挥发损失而具有一定优势，可能抵消较多数量小容量长管拖车远距离运输带来的劣势。

由图 9c—9e可知，当3家化工厂产氢规模不变，运输至盐穴的距离增加时，可能会出现产氢规模较小的化工厂不再处理并运输副产氢的情况。当运输距离进一步增加，不同化工厂可能根据产氢规模从小到大依次出现不再处理并运输副产氢的情况。该结果表明：距离盐穴较远且副产氢产量较小的化工厂不会成为副产氢供应链中的供应商，因为其处理并运输副产氢的成本高于出售副产氢带来的收入。

对比图 9a—9e可知，随着单车单时段运行成本升高，图中字母空缺处越来越少，即在越来越多的场景下，化工厂之间会进行联合运输。由此可知，单车单时段运行成本升高会促进化工厂间进行联合运输，以此降低运输成本，验证了本文所设想的多家化工厂联合运输具有经济性。此外，随着单车单时段运行成本升高，由于化工厂运输同等质量氢气时所需液罐车数量少于长管拖车，因此化工厂在越来越多的场景下采用液氢运输方式以降低运输成本。

由图 9f—9j可知，在同一化工厂产氢规模和运输距离场景下，随着单车单时段运行成本升高，联合运输带来的收益增加，该结果同样验证了本文所设想的多家化工厂联合运输具有经济性。但是需要指出的是，化工厂产氢规模、运输至盐穴的距离和联合运输带来的收益间的关系是非单调的，这说明影响联合运输收益的因素复杂，应当在具体的产氢规模和运输距离场景下进行分析。

(2) 不同单车单时段运行成本下市场需求对优化结果的影响。

为得到副产氢供应链规划及调度的共性指导，本文在所提副产氢供应链的基础场景上，批量修改市场副产氢线性需求函数参数D₀和K₀，产生了更多场景，进而可研究市场需求变化对售氢收益影响的一般性结论。D₀和K₀分别表示市场副产氢供应量为0时的副产氢需求量和需求曲线斜率。其中，D₀分别设置为120、240、360、480、600、720、840、960、1 080和1 200，K分别设置为0.24、0.48、0.72、0.96、1.20、1.44、1.68、1.92、2.16和2.40， $\mathcal{H}$中包括100、200、300、400、500、600时，将产生600个不同的副产氢需求量场景。针对上述600个场景，分别求解所提优化模型，得到上游供应链售氢收益，结果如图 10所示。

根据图 10a—10f中任意一张图可知，当D₀越大，K₀越小，单车单时段运行成本K₃越小时，由优化问题得到的最优决策下的售氢收益越大。由此可知，当副产氢需求量越大，或副产氢需求量对单位氢气零售价格的灵敏度越小时，副产氢上游供应链收益越大。此外，为了提升售氢收益，上游供应链应改进氢气运输技术，减小单车单时段运行成本。

4 结论

首先，本研究提出了盐穴收储化工厂副产氢的商业模式，并对该商业模式下盐穴荷气状态变化、化工厂氢气加工及运输过程中的动态过程进行详细建模。其次，将该商业模式下的副产氢供应链决策问题建模成一个优化问题，并给出该优化问题的求解方法。最后，利用所提模型对1个盐穴和3家化工厂构成的副产氢供应链进行分析。分析结果表明：

1) 在不同的供应链场景参数和氢气运输与净化的相关经济技术参数下，盐穴-化工厂供应链可以实现净收入，存在一定利润空间，因此本文所提商业模式成立。

2) 本文所提优化模型在规划阶段优化运输路径、运输方式和氢处理装置容量，以及在运行阶段根据电价波动优化氢处理量，提升了副产氢上游供应链的收益。

3) 通过灵敏度分析展示了本文所提联合运输在单车单时段运行成本变化下带来的收益。此外，在某一给定的市场需求下，存在一个最佳管道输送容量，在该最佳值以下，盐穴-化工厂收益随着输送容量提高而增加；达到该最佳值后，增加管道输送容量将不再继续提高收益。

4) 通过批量修改化工厂副产氢的生产规模、化工厂至盐穴的距离及单车单时段运行成本，得出化工厂间运输路径及运输方式的一般性结论。当不同化工厂间副产氢生产规模存在一定差距时，中小产量的化工厂在多数场景下更倾向于联合运输，而高产量的化工厂倾向于独立运输；当单车单时段运行成本较高时，随着生产规模扩大和运输距离增加，以盐穴为运输目的地的化工厂从气氢运输方式逐渐转变为液氢运输方式，但单车单时段运行成本较低时该规律不一定满足；随着单车单时段运行成本升高，距离盐穴较远且副产氢产量不大的化工厂不会成为副产氢供应链中的供应商，因为其处理并运输副产氢的成本高于出售副产氢带来的收入；单车单时段运行成本升高会促进化工厂间进行联合运输，以此降低运输成本。

5) 在由1个盐穴和3个化工厂构成的副产氢供应链的基础场景上，批量修改市场副产氢线性需求函数参数D₀和K，得出市场需求变化对售氢收益影响的一般性结论，即当副产氢需求量越大，或副产氢需求量对单位氢气零售价格的灵敏度越小时，副产氢上游供应链收益越大。上游供应链应改进氢气运输技术，减少单车单时段运行成本以提升售氢收益。

本文所提商业模式为化工厂和盐穴提供了新的收入来源，通过减少工业废气排放提高了资源利用率，降低了废气排放对环境的影响，有利于盐穴自然资源的合理利用。