输入不确定性下乙烯裂解炉低氮氧化物排放与高乙烯收率的动态保真度自适应鲁棒优化

《Journal of Industrial and Engineering Chemistry》：Dynamic-fidelity adaptive robust optimization for low-NOx emissions and high ethylene yield in ethylene cracking furnaces under input uncertainty

　　

• PSO-BP代理模型使用166个耦合工况数据和外部留一集进行评估。• DARO将自适应搜索与直接蒙特卡洛保真度校正相结合。• 两个决策变量揭示了稳健的NO—乙烯收率之间的折衷关系。• 三十组配对实验在ZDT1—ZDT3基准上比较了DARO与NSGA-II。• 耦合模型的重新计算支持了推荐的工业操作点。

引言
背景
乙烯是大宗化学品、精细化学品和新能源材料的基础平台化学品[1]。在石油化工路线中，乙烯裂解炉提供了高温环境，将原料（如石脑油）热分解为较轻的烯烃，如乙烯和丙烯。在炉子运行过程中，燃料气燃烧产生大量烟气，其中氮氧化物（NOx）占NOx总排放量的95%以上，是主要的大气污染物，通常以NO的形式出现[2]。针对NOx的环境影响，工业炉被设置了越来越严格的排放上限（如NO ≤ 100 ppm）[3]。为满足裂解炉的环保要求，有必要在乙烯生产过程中的燃烧阶段控制NOx的生成[4]。

工业裂解炉面临不可避免的输入不确定性——例如设定值与实际燃料流量之间的偏差、燃料组成的波动以及风门位置的变化——这些源于测量噪声、执行器限制和过程扰动。炉子将火箱中的湍流反应流、辐射/对流传热和复杂化学反应与盘管中的裂解反应耦合在一起。因此，不确定性引起的速度、温度和气相场的波动会改变管壁热通量和反应物停留时间，从而改变裂解深度和选择性。这既影响出口组成（如乙烯收率），也影响燃烧产生的NOx。在实际操作中，操作人员面临多目标折衷：在保持高乙烯收率和对扰动鲁棒性的同时，降低平均NOx排放及其变异性。这促使在不确定性下进行稳健优化，而非纯粹的确定性调节。

文献综述
由于乙烯裂解炉内复杂湍流和剧烈燃烧化学反应的综合作用，直接实验和设备上的优化耗时且成本高昂；在某些情况下，物理测试甚至不可行。因此，数值模拟变得不可或缺。过去十年间，计算流体力学（CFD）[5]提供了一种有效的数值工具，用于解决耦合的流动、燃烧、传热和组分输运过程。Hirsch [6]在其工作中将CFD定义为"一种使用计算机对流体流动进行数值模拟的技术方法"，它以相对较低的成本实现大量精确预测，从而减少了对昂贵试验的依赖。在裂解炉研究中，大多数CFD分析采用确定性模型，在固定边界和工况下指定热力学、动力学和输运参数。例如，Liu等人[7]在CFX上使用标准k-ε模型获得了详细的炉内速度和温度场，与工厂数据吻合良好。Ni等人[8]使用CFD数值模拟方法评估了基于标准k-ε模型、RNG k-ε模型和可实现k-ε模型构建的湍流模型。然而，一个关键局限仍然存在：真实的工业系统表现出多维不确定性（如输入、物性、模型参数），而确定性模型忽略了这些，可能导致模拟与工厂行为之间的偏差。这一差距促使在不确定性下研究CFD。

在这一背景下，模拟中的不确定性量化（UQ）方法可大致分为非统计（谱/代理）方法和统计（基于采样）方法。在非统计方面，Avdonin等人[9]应用非侵入式多项式混沌展开（NIPCE）对预混火焰的热声稳定性进行正向UQ和敏感性分析，并将结果与蒙特卡洛（MC）进行验证。类似地，Yousefian等人[10]在实际旋流稳定燃烧器中集成了NOx的不确定性、敏感性和可靠性分析，强调了UQ对燃烧工具链日益增长的重要性。在统计方面，MC仍是经典基准：在大样本情况下，MC量化不确定性而不简化模型假设，能处理高维度和非线性，并产生可信的变化范围以供控制和优化参考。其通用性在各领域均有体现[11]。Toli?等人[12]通过MC评估了光伏电池填充因子的不确定性。Tong [13]将MC与《测量不确定度表示指南》（GUM）进行了对比，揭示出不容忽视的差异，表明MC更为可靠。Zhu [14]通过避免基于导数的传播展示了MC的计算便利性。Hu等人[15]使用MC推导了极端风标准差的广义表达式。总之，虽然MC可能需要大量样本，但当模型复杂且解析处理困难时特别适用——这在裂解炉CFD中是常见情况——使得基于MC（或混合代理-MC）的UQ高度适用。

尽管不确定性量化取得了最新进展，但过程系统中的传统优化仍然经常通过确定性离散化来处理不确定性，即在几个预设场景下求解目标[16]。这种做法将随机变异性简化为静态点，忽略了随机扰动、参数漂移和环境噪声，无法显式优化稳健性[17][18]。相比之下，稳健优化纳入了不确定性集合或分布，寻求平衡性能与韧性的解决方案（例如，在最大化乙烯收率的同时最小化NO均值及其分散度）。这种折衷在裂解炉的辐射段尤为尖锐，其中排放量和收率在不确定性下呈拮抗响应。因此，多目标进化算法（MOEAs）具有吸引力：作为黑箱、无梯度优化器，它们近似Pareto前沿并揭示决策折衷[19][20]。其中，非支配排序遗传算法II（NSGA-II）[21]因快速非支配排序、拥挤距离多样性及已验证的工业影响力而备受瞩目——例如，Safikhani等人[22]使用CFD耦合NSGA-II优化了旋风分离器几何形状；Zarifi等人[23]通过多目标NSGA-II自动调谐PVT数据；Liu等人[24]通过在岛模型并行框架中集成SPEA2+NSGA-II提高了可扩展性。

近期研究已将预测优化应用于动态控制系统。Hao等人[25]将模型预测控制与序列凸规划结合用于车辆轨迹规划，而Shi等人[26]开发了带有步补偿的预测功能控制用于非线性与扰动下的电机转速调节。与这些闭环方法不同，本研究在概率扰动下对稳态炉操作条件进行离线多目标优化。所得的稳健操作点可作为未来炉控制的参考目标。

对于这个稳态稳健设计问题，传统NSGA-II提供了有用的基准。然而，嵌套蒙特卡洛公式需要显式的可行性管理、近似误差校正以及搜索算子间的自适应分配，这促使了3.3节中介绍的DARO方法的发展。

动机与主要贡献
为了解决不确定性工况下确定性优化的局限性，本研究开发了PSO-BP-MC稳健设计工作流程和一个重新设计的动态保真度自适应稳健优化器（DARO）。DARO是搜索算法的名称；它不是可调节的两阶段稳健优化模型，也不意味着时变的工厂动态。主要贡献如下：

1. 使用包含166个真实工况的耦合CFD-Coilsim1D数据库来建模五个炉输入条件下的NO排放和乙烯收率，包括新添加的配对温度工况。
2. 随机划分的PSO-BP代理模型与100,000样本的蒙特卡洛传播和Sobol分析相结合，以量化响应不确定性和有影响的燃料流量输入。
3. 工业稳健设计问题扩展到两个决策变量，并明确写为均值-标准差多目标问题，约束条件为E(NO)+3σ(NO)≤100 ppm。
4. DARO结合了可行性感知的存档、基于在线信用的算子组合、直接嵌套MC保真度校正和自适应学习的约束安全裕度。

本文结构
本文其余部分组织如下。第2节建立了概率稳健设计问题。第3节介绍了耦合代理-MC工作流程和重新设计的DARO。第4节描述了166工况数据集、不确定性和实验设置。第5节报告了代理、敏感性、基准和工业结果。第6节对论文进行总结。

问题描述
乙烯裂解炉为烃类热解提供了高温环境，与产物形成和污染物生成密切相关。在典型的炉系统中，辐射段是燃料燃烧、传热和裂解反应强耦合的主要区域。因此，其运行状态对燃料供应、进料条件和入口温度的变化高度敏感。

乙烯裂解是一个高度精密的过程，需要严格控制多个操作参数，包括原料和燃料组成、燃料质量流量和预热温度。然而在实际中，设备退化、测量误差和控制偏差等因素不可避免地导致校准值与实际输入值之间的差异。这些偏差表现为随机波动，通常称为工业输入不确定性。

这类不确定性的传播如图1所示。系统输入的波动不可避免地平移为输出的变异性。对于乙烯裂解炉，这些波动破坏了操作稳健性，使得同时实现关键目标变得困难：减少NOx排放、保持高乙烯收率和增强抗扰动能力。因此，不确定性下的优化变得不可或缺。在本研究中，稳健优化问题通过将不确定性量化获得的统计响应作为优化目标来建立。具体而言，NO排放均值、NO排放标准差和乙烯收率分别用于表示排放水平、操作稳健性和经济性能。

下载：下载高分辨率图片（183KB）下载：下载全尺寸图片
图1. 不确定性传播过程。

然而，作为基准采用的传统NSGA-II不会显式地根据种群状态的变化来调整搜索算子，也不会校正快速稳健响应近似引入的保真度误差。为了解决这些局限性，本研究开发了用于输入不确定性下炉优化的DARO。DARO结合了：(i) 分层种群初始化与可行性感知的选择和存档；(ii) 在线信用分配机制，自适应地在SBX-多项式变异、差分变异、停滞响应Beta变异和目标空间Pareto间隙修复之间分配后代；以及(iii) 动态保真度验证，其中选定候选者通过CFD训练的代理模型进行直接蒙特卡洛传播重新评估，以更新局部残差校正和数据驱动的约束安全裕度。这些机制共同设计用于平衡收敛性、Pareto前沿多样性、约束可靠性和计算成本。

方法论
乙烯裂解炉模拟
基本控制方程
在本研究中，我们使用CFD进行炉模拟，使用Coilsim1D [27]进行反应管模拟。在炉模拟中，采用了标准k-ε湍流模型；有限速率/涡耗散模型（FR/ED）用作燃烧模型[28]；DO模型用作辐射模型，NO计算使用NOx后处理模型进行[29]。炉出口的NO浓度以干基体积分数（ppm）报告。未对参考氧浓度进行校正。CFD预测和工厂测量采用相同基准处理。反应管模拟中使用的方程选择自先前研究[29][30]。本实验模拟中使用的公式见表1。除非另有说明，方程中使用的所有符号都在术语表中定义。

表1. CFD模拟公式。

标题 | 方程
炉腔模拟 |
连续性方程 | ?/?xi(ρUi)=0 (1)
能量方程 | ?/?xi[Ui(ρE+p)]=?/?xj(keff?T/?xj-Σj hj→Jj+Uiμeff[(?Uj/?xi+?Ui/?xj)-2/3?Ul/?xlδij])+Sh (2)
动量方程 | ?/?xj(ρUiUj)=?peff/?xi+?/?xj[μeff(?Ui/?xj+?Uj/?xi-2/3δij?Ul/?xl)] (3)
湍流模型 | ?/?xi(ρkUi)=?/?xj[(μ+μt/σk)?k/?xj]+Gk-ρε (4)
| ?/?xi(ρεUi)=?/?xj[(μ+μt/σ?)??/?xj]+ε/k(C1εGk-C2ερε) (5)
组分输运方程 | ?ρUjYi/?xj=?/?xj[(ρDi,m+μt/Sc t)?Yi/?xj]+Ri (6)
燃烧模型 | Ri=Mw,iΣ NRr=1 R?i,r (7)
| R?l,r=Γ(ν″i,r-ν′i,r)(kf,r Π Nj=1[Cj,r](η′j,r+η″j,r)) (8)
| Ri,k=ν′i,kMw,iAβεk min(YR/ν′R,k Mw,R) (9)
| Ri,k=ν′i,kMw,iABρεk ΣPYP Σ Nj ν″j,k Mw,j (10)
辐射模型 | ?·(I(r→,s→)s→)+(α+σs)I(r→,s→)=αn2σT?/π+σs/4π ∫4π0 I(r→,s'→) Φ(s→·s'→) dΩ' (11)
反应管模拟 |
连续性方程 | dFj/dz=πdt2/4 Rj (12)
能量方程 | dT/dz=πdt Σj Fj cpj q(z)+πdt2/4 Σj Fj cpj Σj(-ΔH0fj)Rj (13)
机械能方程 | dpt/dz=(2fdt+ζπrb)ρg u2+ρg u du/dz (14)

炉和反应管耦合模型
炉模型和反应管模型以耦合方式求解，以描述外部传热与内部裂解反应之间的相互作用。在耦合过程中，信息反复交换，直到管壁温度在连续迭代之间不再显著变化。CFD模型提供炉侧场，包括气体流动、燃烧行为和辐射热通量，而Coilsim1D用于计算反应管内工艺气体的裂解反应。在本研究中，管壁温度被采用为耦合炉-反应器计算的收敛指标[31]。耦合计算方法如图2所示。

下载：下载高分辨率图片（88KB）下载：下载全尺寸图片
图2. 炉和反应管之间的耦合计算方法。

PSO-BP模型与蒙特卡洛不确定性传播
尽管BP神经网络可以表示复杂的非线性关系，但其梯度下降训练对初始权重和偏置敏感，可能收敛缓慢或陷入局部极小值。为缓解这一局限性，PSO用于识别有利的初始化，之后BP进一步细化网络参数。图3总结了数据划分、PSO-BP模型开发和验证过程。外部留一集与预处理和模型开发保持分离，仅用于最终评估。一旦验证通过，在蒙特卡洛不确定性传播期间，代理模型取代耦合的CFD和Coilsim1D模拟器，避免了每种操作条件需要约10—12小时的重复模拟。详细模型设置见第4节。

下载：下载高分辨率图片（197KB）下载：下载全尺寸图片
图3.

## PSO-BP代理模型工作流程

PSO是一种基于群体智能的优化方法，其中一组候选解（粒子）通过共享信息并调整轨迹来迭代搜索解空间。每个粒子根据两个引导项更新其速度和位置：

- 粒子自身发现的个体最优位置（pbest）
- 整个群体发现的全局最优位置（gbest）

通过这种机制，PSO通过平衡探索与开发来加速收敛并避免陷入局部最优。在BP训练之前，采用粒子群优化来搜索适当的初始连接权重和阈值参数，从而减少了对随机初始化的依赖。这种混合方法显著提高了收敛稳定性和预测精度。

PSO的速度和位置更新方程如下：

(15) $v_i^{k+1} = w \times v_i^k + c_1 \times r_1 \times (pbest_i - x_i^k) + c_2 \times r_2 \times (gbest_i - x_i^k)$

(16) $x_i^{k+1} = clip(x_i^k + \alpha \cdot v_i^{k+1}, x_{min}, x_{max})$

每个粒子的质量通过在随机留出的验证子集上使用均方根误差进行评估：

(17) $J(\theta) = \sqrt{\frac{1}{n_{val}} \sum_{i=1}^{n_{val}} [y_i - \hat{y}(x_i; \theta)]^2}$

其中 $\theta$ 表示粒子所代表的完整连接权重和偏差向量，$n_{val}$ 为验证样本数量，$y_i$ 和 $\hat{y}(x_i; \theta)$ 分别为样本 $i$ 的仿真目标值和代理模型预测值。PSO搜索使验证均方根误差最小化的参数向量。PSO获得的最佳参数向量随后用于初始化BP微调，在微调过程中监测验证误差以提高训练稳定性并控制过拟合。

训练完成后，MC仿真生成的不确定输入样本使用PSO-BP代理模型进行评估，而无需反复调用CFD和Coilsim1D耦合求解器。得到的氮氧化物排放量和乙烯收率样本用于估计它们的均值和标准差。这些统计量为后续的敏感性分析和鲁棒多目标优化提供了基础。

## 动态保真度自适应鲁棒优化器

DARO被重新设计为用于嵌套不确定性评估的约束多目标优化器。其名称代表两个算法特性：算子概率根据观察到的后代性能自适应调整，响应保真度通过按计划的嵌套MC评估进行动态修正。它并不代表两阶段可调鲁棒优化或时域动态优化。

### 可行性感知排名与双归档

初始种群通过分层随机采样[32]生成。给定种群规模为 $N$，个体 $i$ 在第 $j$ 维上的归一化坐标 $z_{ij}$ 和决策值 $x_{ij}$ 生成如下：

(18) $z_{ij} = \frac{\pi_j(i) + r_{ij}}{N}$

(19) $x_{ij} = l_j + z_{ij} \cdot (u_j - l_j)$

其中 $\pi_j$ 是 $\{0, \cdots, N-1\}$ 的独立随机排列，$r_{ij} \in [0, 1]$。该初始化在每个决策维度的每个分层中放置一个候选解。

约束处理遵循Deb[33]的可行性原则。设 $g^t(x)$ 为当前快速/修正模型预测的约束值，$s^t$ 为自适应安全裕度。第 $t$ 代候选解 $x$ 的修正约束违反度计算如下：

(20) $v^t(x) = \max\{0, g^t(x) + s^t\}$

违反度较小的候选解优先保留。如果两个候选解具有相同的可行性状态，则依次比较非支配等级和拥挤距离。工业约束 $g(x)$ 在第4.1.3节中定义。

DARO维护一个全局非支配归档和一个鲁棒可行归档。全局归档保持目标空间探索，而鲁棒归档仅包含满足 $v^t(x) = 0$ 的唯一候选解。当归档超过容量时，使用拥挤距离进行截断。鲁棒归档为最终直接蒙特卡洛验证提供候选解。

### 反馈驱动的算子组合

使用四个后代生成器：带多项式变异的模拟二进制交叉（SBX-PM）、差分变异[34]、停滞响应Beta变异和目标空间Pareto间隙修复。前两者提供标准重组和方向搜索，后两者响应停滞和稀疏Pareto前沿覆盖。

设 $q^t$ 为连续未增加第一非支配前沿规模的代数。搜索尺度和Beta变异因子为：

(21) $a^t = \min(1, 0.2 + \frac{q^t}{q_{max}}), \quad \beta^t = \max(1, \frac{5}{a^t}), \quad B^t \sim Beta(1, \beta^t)$

使用独立方向变量 $r \sim U(0,1)$，决策维度 $j$ 的变异值如下：

(22) $x'_j = \begin{cases} x_j + (u_j - x_j) \cdot B^t, & r < 0.5 \\ x_j - (x_j - l_j) \cdot B^t, & r \geq 0.5 \end{cases}$

随着停滞的积累，$a^t$ 增加且 $\beta^t$ 减小，从而增加预期变异步长。因此，该算子响应停滞扩展探索，而不是随代数单调减小变异半径。

对于Pareto间隙修复，第一前沿上的目标向量被归一化并按第一个目标排序。相邻对中欧式目标空间间隙最大的对 $(a, b)$ 识别如下：

(23) $(a, b) = \arg\max_{(i, i+1)} \| \tilde{F}_{i+1} - \tilde{F}_i \|_2$

识别该对后，在其决策空间中点附近生成新候选解：

(24) $x' = clip[\frac{x_a + x_b}{2} + \varepsilon, l, u], \quad \varepsilon_j \sim \mathcal{N}(0, [0.05 \cdot a^t \cdot (u_j - l_j)]^2)$

算子概率通过在线信用分配[35]更新。在第 $t$ 代，算子 $k$ 的奖励从后代存活、非支配性、鲁棒可行性和基于拥挤度的多样性计算如下：

(25) $R_k^{(t)} = 0.40 \cdot S_k^{(t)} + 0.25 \cdot N_k^{(t)} + 0.20 \cdot F_k^{(t)} + 0.15 \cdot D_k^{(t)}$

为减少代间波动，瞬时奖励通过指数平滑处理：

(26) $M_k^{(t)} = \rho \cdot M_k^{(t-1)} + (1 - \rho) \cdot R_k^{(t)}$

平滑信用首先转换为归一化softmax分数：

(27) $q_k^{(t)} = \frac{\exp\{4[M_k^{(t)} - M_{max}^{(t)}]\}}{\sum_{h=1}^{K} \exp\{4[M_h^{(t)} - M_{max}^{(t)}]\}}$

然后引入非零概率下界 $p_{min}$ 获得选择概率：

(28) $p_k^{(t)} = p_{min} + (1 - K \cdot p_{min}) \cdot q_k^{(t)}$

因此，成功的算子获得更多后代，而 $p_{min}$ 防止任何算子被永久淘汰。

### 动态保真度修正与自适应安全裕度

普通搜索评估使用快速响应模型。在预定代数或搜索停滞之后，选定的非支配候选解通过直接蒙特卡洛传播经PSO-BP代理模型重新评估。此验证绕过了快速鲁棒响应模型，但不需要为每个蒙特卡洛样本反复执行CFD-Coilsim1D仿真。

对于已验证候选解 $x_i$，直接值减快速值的目标残差和修正目标预测定义如下：

(29) $r_i = F^{direct}(x_i) - F^{fast}(x_i), \quad F^{corr}(x) = F^{fast}(x) + \hat{r}(x)$

为估计 $\hat{r}(x)$，在归一化决策空间中应用高斯核加权[36]，结果修正在已验证邻域外进行衰减：

(30) $d_i(x) = \| \frac{x - x_i}{u - l} \|_2, \quad w_i(x) = \exp\left[-\frac{d_i^2(x)}{2h^2}\right], \quad \hat{r}(x) = \exp\left[-\frac{d_{min}^2(x)}{2h^2}\right] \cdot \frac{\sum_i w_i(x) r_i}{\sum_i w_i(x)}$

相同的直接评估用于计算约束预测误差 $e_i$ 和约束低估的安全裕度 $s^t$：

(31) $e_i = g^{direct}(x_i) - g^{pred}(x_i), \quad s^t = Q_{0.95}(\{\max(0, e_i)\}_{i=1}^{n^t})$

仅正误差有贡献，因为它们代表被低估的约束风险。每次保真度更新后，当前种群重新评估并使用修正后的目标和安全裕度更新排名。保真度间隔、停滞阈值、批大小、带宽和直接评估预算在实验设置中报告。优化后，鲁棒归档候选解再次通过直接蒙特卡洛检查，仅报告直接验证可行的解。

## 实验实施

### 实验设计

#### 乙烯裂解炉的几何与边界条件

以一工业乙烯裂解炉为研究对象。炉体尺寸为11.544 m（长）× 2.693 m（宽）× 19.821 m（高），配备24个多级底烧器和24个侧壁烧器。对炉体的六分之一对称截面进行建模以进行数值仿真，以减少计算成本。该截面包含48根垂直反应管，分为六组，每组八根。主要几何和操作参数列于表2。边界条件设置见图4。

**表2. 六分之一乙烯裂解炉模型的固定参数。**

| 参数 | 值 |
|------|------|
| **炉体边界条件** | |
| 过量空气系数 | 1.18 |
| 炉体长度（x方向）(m) | 1.924 |
| 炉体宽度（y方向）(m) | 2.693 |
| 炉体高度（z方向）(m) | 19.821 |
| 燃料组分比（mol%）H?/CH?/C?H?/CO | 5/93.8/0.6/0.6 |
| **反应管条件** | |
| 反应管数量 | 8 |
| 程数 | 2 |
| 反应管总长度 (m) | 26.59 |
| 入口管内径 (10?3m) | 52 |
| 出口管内径 (10?3m) | 60 |
| 入口管壁厚 (10?3m) | 7 |
| 出口管壁厚 (10?3m) | 8 |
| 反应管入口温度 (K) | 853 |
| 反应管入口压力 (Pa) | 300,000 |
| 反应管出口压力 (Pa) | 90,000 |
| 裂解原料组分（wt%）：正构烷烃/环烷烃/烯烃/芳香烃 | 58.5/31.5/0.0/10.0 |

#### 不确定性量化

修订后的工业优化将可控决策变量与不确定操作输入区分开来。过量空气系数和空气预热温度被视为决策变量，而主底烧器、副底烧器和侧壁烧器的燃料流量被视为不确定输入。决策边界、数据集覆盖范围、燃料流量范围和相应单位总结于表3。

**表3. 输入变量分布表。**

| 参数类型 | 参数名称 | 参数分布 | 均值 | 标准差 |
|----------|----------|----------|------|--------|
| 不确定参数 | 底烧器主燃料流量（kg/s） | 高斯 | 1.14e-02 | 1.14e-04 |
| | 底烧器副燃料流量（kg/s） | 高斯 | 6.86e-02 | 6.86e-04 |
| | 侧壁烧器燃料流量（kg/s） | 高斯 | 6.51e-02 | 6.51e-04 |
| 决策变量 | 过量空气系数 | 上下界 [1.0, 1.3] | — | — |
| | 空气预热温度（K） | 上下界 [293, 310] | — | — |

三个燃料流量输入的均值取自相应工厂操作设定点。遵循GUM框架[37]，通过基于过程知识和操作经验的B类评定分配1%的相对标准不确定度，并假设服从高斯分布。由于三个燃料流量在正常操作中分别调节，其短期波动被视为相互独立。具体数学表示为第 $i$ 个不确定参数服从均值为 $\mu$、标准差为 $\sigma = 0.01\mu$ 的高斯分布：

(32) $X \sim N(\mu, (0.01\mu)^2)$

确定不确定参数后，使用拉丁超立方采样（LHS）方法生成训练样本，以确保在整个参数空间上均匀且有代表性的分布。这些样本随后使用PSO-BP神经网络进行建模和训练，建立输入不确定性与关键炉体输出之间的非线性映射。本文使用的PSO-BP超参数如表4所示。

**表4. PSO-BP超参数。**

| 超参数 | 值 |
|--------|------|
| 网络架构 | 5-5-2，含一个隐藏层 |
| 激活函数 | 隐藏层Tanh，输出层线性 |
| PSO惯性权重 | 1.0 |
| PSO认知系数 | 4.494 |
| PSO社会系数 | 4.494 |
| PSO群体/迭代 | 5/50 |
| PSO LM阻尼 | 初始值10?3；接受步后除以3，否则乘以10 |
| PSO L2系数 | 1×10?? |
| PSO目标函数 | 权重和偏差的验证RMSE |
| BP优化器 | Levenberg-Marquardt，最多1000轮 |
| 学习率 | 0.01 |

PSO-BP代理模型训练完成后，蒙特卡洛仿真生成的输入样本直接由代理模型评估，而无需反复进行CFD计算。得到的响应样本随后用于获取氮氧化物排放和乙烯收率的统计特征，如均值、标准差和变异系数。这样，输入波动引起的响应不确定性可以被高效地量化。为进一步识别响应变异的来源，对不确定输入变量进行Sobol敏感性分析。该方法将输出方差归因于各单独输入及其耦合效应。对于具有 $d$ 维不确定输入向量的模型，方差分解可写为：

(33) $V(Y) = \sum_{i=1}^{d} V_i + \sum_{1 \leq i \leq j \leq d} V_{ij} + \cdots + V_{1,2,\cdots,d}$

(34) $V_i = V[E(Y|\xi_i)]$

(35) $V_{ij} = V[E(Y|\xi_i, \xi_j)] - V_i - V_j$

基于此方差分解，一阶Sobol'指数 $S_i$ 定义为：

(36) $S_i = \frac{V_i}{V(y)}$

式(36)表示单个参数独立贡献的输出方差比例。总Sobol'指数为：

(37) $S_T^i = \frac{V(Y) - V_{-i}}{V(Y)} = 1 - \frac{V[E(Y|\xi_{-i})]}{V(Y)}$

### 基于PSO-BP-MC的DARO框架

第5.2节的不确定性量化结果表明，不确定燃料流量输入引起的乙烯收率变异性远小于氮氧化物排放的变异性。因此，乙烯收率的波动不作为额外的鲁棒性目标引入。相反，乙烯收率本身作为经济性能目标保留，而氮氧化物排放的均值和标准差均被考虑，分别描述排放水平和其对输入不确定性的敏感度。

设决策向量定义为：

(38) $x = (\lambda, T_{air})^T$

其中 $\lambda$ 和 $T_{air}$ 分别表示过量空气系数和空气预热温度。三个不确定烧器燃料流量收集于 $\xi$ 中，其概率模型在第4.1.2节中定义。对于给定决策向量 $x$，氮氧化物排放的均值和标准差表示为：

(39) $\mu_{NO}(x) = E_\xi[Y_{NO}(x, \xi)], \quad \sigma_{NO}(x) = \sqrt{Var_\xi[Y_{NO}(x, \xi)]}$

由于不确定性分析表明乙烯收率波动相对于其生产水平可忽略不计，乙烯目标由在名义不确定输入处评估的收率表示：

(40) $Y_{C_2H_4}(x) = Y_{C_2H_4}(x, \bar{\xi})$

基于这些定义，工业鲁棒多目标优化问题表述为：

(41) $\min_x F(x) = [\mu_{NO}(x), \sigma_{NO}(x), -Y_{C_2H_4}(x)]^T$

s.t. $\mu_{NO}(x) + 3\sigma_{NO}(x) \leq 100 \text{ ppm}$

$1.0 \leq \lambda \leq 1.3$

$293 \leq T_{air} \leq 310 \text{ K}$

$Y_{C_2H_4}$ 前的负号将乙烯收率最大化转换为多目标优化器使用的常见最小化形式。约束 $\mu_{NO} + 3\sigma_{NO} \leq 100 \text{ ppm}$ 引入了氮氧化物排放的保守均值-标准差裕度。它作为鲁棒工程约束使用，而非被解释为精确的概率保证。

为将鲁棒优化应用于本研究中的乙烯裂解炉，必须从不确定性量化结果构建高效的目标评估模块，以决策向量为输入。该模块允许优化算法在可行决策空间上进行全局搜索，而无需对每个候选解反复执行完整的不确定性传播程序。为构建训练数据集，过量空气系数和空气预热温度在可行范围内以均匀间隔水平采样，形成二维决策网格。在每个采样决策点，PSO-BP代理模型嵌入蒙特卡洛程序中，计算鲁棒优化模型所需的统计响应。所得的决策-响应数据用于训练基于PSO-BP的鲁棒响应评估器，该评估器封装了完整的PSO-BP-MC不确定性量化程序。在DARO优化过程中，该封装模块提供快速的目标评估，而选定的候选解定期使用原始PSO-BP-MC程序重新计算，所得差异用于动态保真度修正。最终Pareto解在报告前也须通过直接PSO-BP-MC验证。不确定性量化与鲁棒优化之间的耦合如图5所示。

**图5. 鲁棒优化模型。**

表5总结了NSGA-II和DARO使用的主要算法参数。

两种算法均采用 30 组配对随机种子、200 个个体的种群规模以及 500 代进化，每次运行产生 100,200 次候选解评估。表 5 仅保留两种算法之间存在差异且对可复现性具有实质影响的参数配置。

表 5. NSGA-II 与 DARO 所采用的参数配置。

| 配置项 | 标准 NSGA-II | 动态保真度自适应鲁棒优化器（DARO） |
|---|---|---|
| 初始化 | 均匀随机采样 | 分层拉丁超立方体采样 |
| 选择 | 基于排名和拥挤距离的二元锦标赛 | 约束感知锦标赛，综合考虑约束违反程度、排名和拥挤距离 |
| 变异算子 | SBX 和多项式变异 | SBX-PM、差分进化、目标差值修复以及停滞响应型 Beta 变异的自适应组合 |
| 主要变异参数 | pc=1.0, ηc=20; pm=1/nvar, ηm=20 | SBX (ηc=20)，PM (ηm=20)，DE (F~U(0.45,0.85), CR=0.90) |
| 算子自适应 | p0=0.2, pmin=0.08 | 奖励记忆系数=0.70 |
| 动态保真度控制 | — | 每 25 代或 15 代停滞之后更新一次；每次更新使用 6 个直接计算点；最大使用 120 次 |
| 鲁棒约束处理 | — | Deb 风格约束处理，并附加 95 分位预测误差安全余量 |
| 外部存档 | — | 全局存档和鲁棒可行 Pareto 存档，各自最多 1,000 个解 |

总体而言，所提出的框架包含六个主要阶段：耦合 CFD–Coilsim1D 数据库的构建、PSO-BP 代理模型的训练与验证、指定燃料流量不确定性的传播、双决策鲁棒响应模型的构建、基于 DARO 的 Pareto 搜索，以及对最终候选解的直接 PSO-BP-MC 验证。相应的代理模型、不确定性和鲁棒优化结果将在第 5 节中展示。

结果与讨论

CFD 模拟分析

CFD 的主要产品收率数值和炉膛信息如表 6 所示，两类数据均与实际工业数据进行了对比。模拟值与实测值之间的微小偏差表明，耦合炉-反应器模型能够可靠地描述乙烯裂解炉内部的热状态。

表 6. 三种代表性工况的运行条件和验证结果。

| 类别 | 参数 | 数据来源 | 工况 1 | 工况 2 | 工况 3 |
|---|---|---|---|---|---|
| 运行条件 | 进料流量 (kg/s) | 模型输入 | 1.2188 | 1.14153 | 1.12856 |
| | 水碳比 (kg/kg) | 模型输入 | 0.50 | 0.5586 | 0.5589 |
| | 空气入口流量 (kg/s) | 模型输入 | 1.60968 | 1.41142 | 1.49774 |
| | 侧壁燃料流量 (kg/s) | 模型输入 | 0.06509 | 0.05444 | 0.06023 |
| | 一次底部燃料流量 (kg/s) | 模型输入 | 0.01140 | 0.01037 | 0.01065 |
| | 二次底部燃料流量 (kg/s) | 模型输入 | 0.06860 | 0.06241 | 0.06411 |
| 验证结果 | 反应器管出口 C2H4 (wt%) | 模拟/实测 | 26.65/28.03 | 26.70/28.08 | 26.31/27.67 |
| | 反应器管出口 C3H6 (wt%) | 模拟/实测 | 11.40/12.34 | 11.42/12.36 | 11.25/12.18 |
| | 反应器管出口温度 (K) | 模拟/实测 | 1157/1161 | 1082/1086 | 1092/1096 |
| | 炉膛出口 NO (ppm) | 模拟/实测 | 60.028/57.687 | 59.79/55.635 | 55.24/51.43 |
| | 炉膛出口 O2 (mol%) | 模拟/实测 | 2.51/2.58 | 1.53/1.58 | 1.28/1.32 |
| | 炉膛出口温度 (K) | 模拟/实测 | 1402/1399 | 1311.6/1308.8 | 1324.1/1321.3 |

上述验证结果充分表明，该耦合模拟方法的模拟精度能够满足工程应用的要求，为后续相关研究和工业优化提供了可靠的模型支持。

代理模型验证、不确定性量化与敏感性分析

数据集构建与 PSO-BP 验证

首先，使用拉丁超立方采样生成了 150 个耦合 CFD-Coilsim1D 工况。为了增强新增的进气预热温度决策维度方向上的数据覆盖度，额外获取了 16 个实际耦合模拟工况。最终数据库包含 166 个实际工况，覆盖规定的输入范围，并用于开发 PSO-BP 代理模型。

在模型拟合之前，使用固定随机种子将 166 个耦合模拟工况随机划分为 133 个开发集和 33 个留出集。在开发集中，106 个工况用于 BP 训练过程中的网络权重更新，其余 27 个工况作为内部调优子集，用于评估 PSO 候选权重并在后续 BP 微调过程中控制早停。图 6 对比了开发集和留出集中耦合模拟响应与相应的 PSO-BP 预测值。

图 6. PSO-BP 代理模型的预测结果。

表 7 中的留出结果显示，最终的 PSO-BP 模型在未被用于拟合或内部调优的样本上仍保持满意的预测精度。对于 NO 排放，留出集的 R2、RMSE 和 MAE 分别为 0.9920、7.994 ppm 和 2.528 ppm。NO 的 RMSE 与 MAE 之间的差异主要归因于一个保留了较大预测误差的留出样本。对于乙烯收率，相应数值分别为 0.9930、0.04480 wt% 和 0.03159 wt%。这些结果在所采用的固定随机划分方式下，提供了对拟合代理模型的外部评估。

表 7. PSO-BP 炉膛代理模型的预测性能。

| 响应变量 | 开发集 RMSE | 开发集 R2 | 留出集 RMSE | 留出集 MAE | 留出集 R2 |
|---|---|---|---|---|---|
| NO 排放 (ppm) | 2.309 | 0.9993 | 7.994 | 2.528 | 0.9920 |
| 乙烯收率 (wt%) | 0.07592 | 0.9766 | 0.04480 | 0.03159 | 0.9930 |

图 7、图 8 使用五折交叉验证对比了 BP、PSO-BP、SVR、GPR 和 RF 五种模型。PSO-BP 在 NO 排放和乙烯收率两个方面均获得了最高的 R2 以及最低的 RMSE 和 MAE，表明其在五种模型中具有最佳的总体预测性能。PSO 搜索仅用于离线网络初始化；训练完成后，PSO-BP 作为一个紧凑的前馈模型运行，可为 100,000 样本的蒙特卡洛计算和重复的 DARO 目标评估提供平滑的响应映射和向量化评估。因此，PSO-BP 的选择同时基于其预测精度和其在大规模代理评估重复使用中的适用性。

图 7. NO 排放预测候选模型的性能评估。
图 8. 乙烯收率预测候选模型的性能评估。

蒙特卡洛收敛与响应不确定性

图 9、图 10 展示了随着蒙特卡洛样本量从 103 增加到 10?，NO 和乙烯收率估计量的归一化均值、标准差和 95 分位数的收敛过程。这三个量在 10? 时均收敛至 1，超出该点后波动可忽略不计。这表明 10? 的样本量足以支持后续不确定性分析中的可靠统计估计。

图 9. NO 排放蒙特卡洛估计量的样本量收敛性。
图 10. 乙烯收率估计量的样本量收敛性。

基于经过验证的代理模型和已收敛的蒙特卡洛采样结果，定量不确定性分析结果汇总于表 8。在三个关键不确定性参数的扰动下，NO 排放指标表现出显著波动特征。其排放浓度的期望值达到 58.8135 ppm，标准差为 2.7876 ppm，相应的变异系数为 4.739%。相比之下，乙烯收率在相同运行条件下表现出很强的稳定性。其平均产量 26.4263 wt% 对应的标准差和变异系数分别为 0.0294 wt% 和 0.111%。这一对比结果表明，裂解炉运行期间的 NO 排放参数对输入条件扰动更为敏感。在不确定性参数影响下，NO 排放性能波动显著，其变异系数比乙烯收率高出一个数量级以上。

表 8. 不确定性定量结果。

| 响应特征 | 均值 | 标准差 | 变异系数 |
|---|---|---|---|
| NO (ppm) | 58.8135 | 2.7876 | 4.739% |
| 乙烯 (wt%) | 26.4263 | 0.0294 | 0.111% |

Sobol 敏感性分析

使用公式 (33) 至 (37) 计算了每个不确定性参数对 NO 和乙烯的敏感性系数，结果如图 11 所示。结果表明，NO 排放和乙烯收率对输入参数扰动的响应特征截然不同。对于 NO 排放，一次燃料流量贡献 79.64%，侧壁燃料流量贡献 14.40%，而二次燃料流量仅占 1.37%。相比之下，乙烯收率的波动主要由二次燃料流量驱动，其一阶 Sobol 指数为 49.99%，其次分别是侧壁和一次燃料流量，为 16.14% 和 6.04%。一阶指数与总效应指数之间的差异表明存在交互效应，尤其是在乙烯收率方面。然而，乙烯收率的整体标准差仅为 0.0294 wt%，其变异系数为 0.111%，表明在所考虑的输入不确定性下，其绝对波动较小。

图 11. NO 排放和乙烯收率的 Sobol 系数。

鲁棒优化结果

DARO 在 ZDT 测试函数上的应用

由于无法获得乙烯裂解炉优化问题的真实 Pareto 前沿，首先采用 ZDT（Zitzler-Deb-Thiele）基准测试函数来检验所提出的 DARO 算法的优化能力。标准 NSGA-II 和 DARO 使用相同的 30 组配对随机种子、200 的种群规模、500 代进化以及每次运行 100,200 次候选解评估。为了定量比较获得的 Pareto 解集，引入了两个常用指标，即世代距离（GD）和分布度（SP）。GD 用于评估所获得解与参考 Pareto 前沿的接近程度，而 SP 反映非支配解的分布均匀性。这两个指标的定义如下：

(42) GD = (1/n)·√(Σ???? d?2)
(43) SP = (Σ???^|Q|-1 |d? - d?|) / ((|Q|-1)·d?)

如图 12 和表 9 所示，DARO 在全部三个测试函数上获得了更低的平均 SP 值。其在 ZDT1 上的平均 GD 与 NSGA-II 相当，而在 ZDT2 和 ZDT3 上更低。最明显的改进出现在 ZDT2 上，其 GD 和 SP 均有下降。

表 9. 不确定性定量结果。（注：原文表 9 内容与表 8 相似，此处按原文保留）

| 响应特征 | 均值 | 标准差 | 变异系数 |
|---|---|---|---|
| NO (ppm) | 58.826 | 2.320 | 3.944% |
| 乙烯 (wt%) | 26.424 | 0.0248 | 0.0975% |

图 12. NSGA-II 与 DARO 在 ZDT1-ZDT3 上 30 次配对运行中 GD 和 SP 的分布。

表 10. 经直接 MC 验证的推荐工业运行点。

| 过量空气系数 | 进气预热温度 | NO 均值 (ppm) | NO 标准差 (ppm) | NO 均值+3σ (ppm) | 乙烯收率 (wt%) |
|---|---|---|---|---|---|
| 1.17549 | 293.04 | 56.2894 | 1.21891 | 59.9461 | 26.5449 |

基于 PSO-BP-MC 的 DARO 框架在乙烯裂解炉上的应用

图 13 展示了在过量空气系数和进气预热温度上，PSO-BP-MC 的响应面。增加过量空气系数通常会降低平均 NO 排放，但同时也会降低乙烯收率，揭示出明显的环境-经济权衡关系。进气预热温度进一步影响炉膛的局部响应和鲁棒性，而 μ_NO+3σ_NO 表面则识别出满足鲁棒排放约束的区域。

图 13. 二维决策空间中 PSO-BP-MC 的响应面。

基于这些响应特征，DARO 被应用于搜索在最小化平均 NO 排放、最小化 NO 波动和最大化乙烯收率三个目标下的鲁棒解。响应模型使用直接 PSO-BP-MC 评估进行动态修正。平均 NO、NO 标准差和乙烯收率的 RMSE 值分别为 0.09029 ppm、0.11594 ppm 和 0.0001356 wt%。

由于三个目标之间存在相互竞争，DARO 产生的是鲁棒可行的 Pareto 解集，而非单一最优解。图 14 展示了一次代表性运行中直接验证候选解和 Pareto 解的分布。为避免主观选择，推荐运行点采用以下规则确定：

(44) d* = argmax_d∈P_MC Y_C2H4(d)，s.t. μ_NO(d) + 3σ_NO(d) ≤ C_base
其中 C_base = 60.028 ppm

P_MC 表示直接验证的 Pareto 解集。该规则防止相对于基线排放水平出现恶化，然后选择乙烯收率最高的解。所得条件为过量空气系数 = 1.17549，进气预热温度 = 293.04 K。使用 100,000 个蒙特卡洛样本验证得到平均 NO 浓度为 56.2894 ppm，标准差为 1.21891 ppm，三倍标准差值为 59.9461 ppm。相应乙烯收率为 26.5449 wt%。预测的 NO 均值接近工业实测值 57.687 ppm，同时乙烯收率保持在 26.5 wt% 以上。因此所选点代表的是保持基线的工业折中方案，而非被推向 100 ppm 边界的解。

图 14. DARO 框架中的 Pareto 解和推荐设计方案。

在所选决策变量处进行了额外的耦合模拟验证，将不确定输入固定在其规定的均值。计算得到的 NO 浓度和乙烯收率分别为 56.7984 ppm 和 26.6358 wt%。这些结果与不确定性平均预测值 56.2894 ppm 和 26.5449 wt% 一致，支持了推荐运行点的可靠性。

结论

本研究开发了一个面向乙烯裂解炉的感知不确定性鲁棒优化框架，集成了耦合 CFD-Coilsim1D 模型、PSO-BP 代理模型、蒙特卡洛不确定性传播和 DARO。代理模型基于 166 个耦合模拟工况建立，NO 的留出 RMSE 为 7.994 ppm，乙烯收率为 0.04480 wt%。不确定性与敏感性分析表明，NO 排放对燃烧器流量扰动更为敏感，而乙烯收率的绝对波动仍然较小。

DARO 的性能在 ZDT1-ZDT3 上通过 30 次配对运行进行了评估。与标准 NSGA-II 相比，DARO 在全部三个测试函数上获得了更低的平均 SP 值，其平均 GD 在 ZDT1 上相当，在 ZDT2 和 ZDT3 上更低。在炉膛应用中，全部 30 次 DARO 运行在直接 PSO-BP-MC 验证后均产生了满足 100 ppm 全局可行性限制的 Pareto 解集，表明其在二维运行空间上具有稳定的搜索性能。

100 ppm 限值被视为全局可行性边界，而非期望运行水平。因此最终运行点选择采用了更严格的保持基线准则：在要求 μ_NO+3σ_NO ≤ 60.028 ppm 的条件下最大化乙烯收率。所选条件为过量空气系数 = 1.17549，进气预热温度 = 293.04 K，平均 NO 浓度为 56.2894 ppm，NO 标准差为 1.21891 ppm，三倍标准差值为 59.9461 ppm，乙烯收率为 26.5449 wt%。

额外的一次耦合仿真验证将不确定的输入固定在其规定的均值处，得到56.7984 ppm NO和26.6358 wt%乙烯，与相应的代理模型预测结果基本一致。所得NO水平接近工业测量值57.687 ppm，而乙烯收率仍保持在26.5 wt%以上。本不确定性分析侧重于规定的操作输入扰动，并以拟合的PSO-BP代理模型为条件。由有限训练数据集和代理模型近似误差引起的认知不确定性未通过鲁棒优化框架进行显式传播，将在未来工作中使用集成或概率代理方法进行探讨。因此，所提出的框架提供了一种透明的手段来识别鲁棒操作折中方案，而不会将解推向排放上限。CRediT作者贡献声明：胡桂花：原创草稿撰写、可视化、软件、方法论、调研、资金获取、概念构思。毛俊明：原创草稿撰写、可视化、软件、方法论、形式分析、数据整理、概念构思。龙健：审稿与编辑、监督、概念构思。胡桂花|毛俊明|龙健