基于正交实验和熵值法耦合人工神经网络的PEMFC多目标优化

Multi-Objective Optimization of PEMFC Based on Orthogonal Experiment, Entropy Method, and Artificial Neural Network

摘要: 质子交换膜燃料电池(PEMFC)性能指标多样，通过多目标优化可以提升PEMFC的综合性能。本文综合运用正交实验(OEM)、方差分析(ANOVA)以及熵值法(EWM)与人工神经网络(ANN)，提出了一种综合评价PEMFC性能的多目标优化方法。首先，选择几个常用的PEMFC性能指标，针对PEMFC的几何参数和操作参数初步选择对PEMFC性能有影响的参数设计正交实验，随后使用COMSOL软件模拟计算获取对应工况下各个性能指标结果。通过方差分析筛选出对目标有显著影响的变量，可以达到减小优化空间、提高优化效率的目的。其次，根据显著变量再次进行正交实验，获取优化空间中具有全局代表性的数据集，随后提出了代表PEMFC综合性能(CP)的表达式，通过熵值法分析了CP中多目标的权重，进一步计算出CP的数学表达式并作为适应度函数。最后将ANN作为计算CP的代理模型，对显著变量进行优化，继而得到一种最优的参数组合使得对应的CP值达到最大，从而最大限度地提升PEMFC的综合性能。利用该方法对梯形挡板流道模型进行了多目标优化分析，在0.4 V电压下，优化后流道模型相比较基础流道模型而言，功率密度和净功率分别提升了53.71%和53.70%，然而阴极压降却仅仅只上升了7.93 Pa，氧气不均匀度降至0.19。上述方法可有效降低优化计算时间，对今后PEMFC的结构设计与操作运行具有一定的指导意义。

Abstract: Proton exchange membrane fuel cell (PEMFC) has various performance indicators, and the comprehensive performance of PEMFC can be improved by multi-objective optimization. In this paper, a multi-objective optimization method for comprehensively evaluating the performance of PEMFC is proposed by comprehensively applying orthogonal experiment (OEM), analysis of variance (ANOVA), and entropy value method (EWM) with artificial neural network (ANN). Firstly, several commonly used PEMFC performance indexes are selected, and orthogonal experiments are designed for the parameters that have an impact on the PEMFC performance by initially selecting the geometrical and operational parameters of the PEMFC, and then simulations are performed using COMSOL software to obtain the results of the performance indexes under the corresponding working conditions. Through the analysis of variance to screen out the variables that have a significant impact on the target, it can achieve the purpose of reducing the optimization space and improving the optimization efficiency. Secondly, orthogonal experiments were conducted again according to the significant variables to obtain a globally representative data set in the optimization space. Subsequently, an expression representing the comprehensive performance (CP) of PEMFC was proposed, and the weights of multiple objectives in the CP were analyzed by entropy method, and the mathematical expression of the CP was further computed and used as the fitness function. Finally, the ANN is used as a proxy model for calculating the CP to optimize the significant variables, and then an optimal combination of parameters is obtained to maximize the corresponding CP value, thus maximizing the comprehensive performance of the PEMFC. A multi-objective optimization analysis of the trapezoidal baffle runner model is carried out using this method, and the power density and net power of the optimized runner model are increased by 53.71% and 53.70%, respectively, compared with the base runner model at 0.4 V, while the cathode pressure drop is increased by only 7.93 Pa, and the oxygen inhomogeneity is reduced to 0.19. The above method can effectively reduce the optimization computation time, which is very useful for the design of future PEMFC structures and the development of the PEMFC.

1. 引言

质子交换膜燃料电池(PEMFC)因其高效能和低排放特性而备受瞩目，被视为未来能源转换的重要设备[1]-[3]。在汽车和电子领域，PEMFC展现出广阔的应用前景，为实现可持续发展提供了重要支持。然而，同时也面临着催化剂成本昂贵、耐久性不足等挑战，这些障碍限制了其在市场上的广泛应用[4]。因此，迫切需要加强对质子交换膜燃料电池性能的研究与改进，以推动其更广泛的应用和发展。

改善流道几何结构对提升PEMFC性能至关重要。例如在流道中添加挡板可以有效强化反应气体的传输效果，提升PEMFC输出性能，许多学者已在这一领域进行了大量研究。Yin等[5]定量研究了梯形挡板的前斜角和后斜角对PEMFC性能的影响。结果表明，梯形挡板的前斜角和后斜角均为45˚时，PEMFC的性能最佳。Perng等人[2]研究了不同几何参数(倾斜角和高度)的梯形挡板对PEMFC性能的影响。研究结果表明，当高度为1.125 mm、倾角为60˚时，PEMFC的净功率比不带挡板的PEMFC提高了约90%。Wang等人[6]研究了具有两种不同结构的PEMFC，发现交错梯形挡板(STBP)流道显著提高了PEMFC的最大净功率，与传统的平行流道和平行梯形障板(PTBP)流道相比，分别提高了6.39%和2.54%。

此外，为了实现质子交换膜燃料电池商业化的目标，还必须确保PEMFC在适当的操作条件下稳定运行。否则会降低燃料电池的输出性能，加速燃料电池的退化，从而大大缩短燃料电池的使用寿命[7]。针对这一挑战，研究人员对这一领域进行了探索，旨在提升质子交换膜燃料电池的性能表现。Yan [8]等实验研究了阴极入口气体流速、阴极入口加湿温度、电池温度等各种操作条件对采用传统流场和交错流场的PEMFC电流密度的影响。实验结果表明，随着阴极入口气体流速、阴极加湿温度和电池温度的升高，电池的性能会有所提高。Kim [9]等研究了相对湿度和反应物的化学计量比对PEMFC水饱和度的影响。结果表明，较高的氢化学计量比降低了阴极处的水饱和度，并增强了水从阴极到阳极的扩散，从而提高了PEMFC的性能。然而，上述研究都使用单一目标来优化PEMFC的性能，这限制了其对不同场景的适应性。

与单目标优化相比，PEMFC的多目标优化具有更强的实用性。由于在PEMFC中，不同的性能指标经常相互冲突，因此，进行多目标优化是一项有意义且具有挑战性的工作。Liu等人[10]以能效和输出功率为目标，通过多目标遗传算法对PEMFC的运行条件和流道结构进行了优化。Kwan等人[11]进行了一项基于NSGA-Ⅱ的优化，以优化PEMFC系统热电发电机设计的系统输出功率和效率。然而，采用多目标优化算法需要不断地去计算和评估适应度函数的值，往往比较耗时。因此，本文提出一种综合评估PEMFC性能的多目标优化方法，可快速、系统地实现电流密度、净功率、压降和氧气不均匀度等多指标的综合优化。基于该方法，本文以具有梯形挡板的PEMFC流道模型为例进行实施，在0.4 V电压下，优化后流道模型相比较基础流道模型而言，功率密度和净功率分别提升了53.71%和53.70%，然而阴极压降却仅仅只上升了7.93 Pa，氧气不均匀度降至0.19。上述方法可有效降低优化计算时间，对今后PEMFC的结构设计与操作运行具有一定的指导意义。

2. 数值模拟

2.1. 多目标优化方法

多目标优化流程如图1所示。首先进行多目标选择，选择几个常用的PEMFC性能指标。针对PEMFC的几何参数和操作参数初步选择对目标有影响的参数设计正交实验，随后通过COMSOL模拟计算来获取对应工况下各个性能指标结果。通过对每个性能指标进行方差分析，筛选出对各个性能指标有显著影响的变量，可以达到减小优化空间、提高优化效率的目的。针对显著变量再次设计正交实验以获取优化空间中具有全局代表性的数据集，并通过COMSOL模拟计算每组工况不同运行电压如(0.4~0.95 V，间隔0.05取值)下的四个性能指标值。随后提出了代表PEMFC综合性能的表达式CP，根据所获得的数据集通过熵值法分析了CP中多目标的权重，进一步计算出CP的数学表达式作为之后多目标优化的适应度函数并计算出每一组工况模拟对应的CP值以供后续ANN模型训练/测试使用，其中每组工况参数及其对应的CP值构成一组数据集。紧接着，在建立ANN模型之后，按照一定比例将获得的所有数据集随机分为训练集和测试集，采用训练集进行训练，通过测试集评估预测性能。最后在确定ANN模型预测的准确性以后重新定义显著变量取值范围，将ANN模型作为计算CP的代理模型获取预测取值范围内所有CP值，继而从中筛选出得到一种最优的参数组合使得对应的CP值达到最大，最大限度的提升了考虑多目标的PEMFC综合性能。

Figure 1. Multi-objective optimization flowchart

图1. 多目标优化流程图

2.2. 几何模型

Figure 2. Schematic of PEMFC runner geometry with trapezoidal baffle structure

图2. 带梯形挡板结构的PEMFC流道几何示意图

采用上述多目标优化方法，本文针对带梯形挡板的PEMFC流道进行多目标优化分析，其模型示意图见图2，具体模型的几何和物理参数见表1。该模型包括双极板(BP)、气体扩散层(GDL)、催化剂层(CL)和质子交换膜(PEM)。一个循环长度(E)包括一条与梯形挡板相连的线性路径。梯形挡板底边的长度等于一个循环长度(E/2)的一半。梯形挡板的两个侧壁角度相等。直线路径和梯形挡板构成一个完整的循环，延伸至整个流道。流道中梯形障板的高度(H)、侧壁角度(θ)和循环次数(N)被定义为可变几何参数。

Table 1. PEMFC model geometry and physical parameters

表1. PEMFC模型几何及物理参数

参数	数值
流道长度/mm	20
流道高度/mm	1
流道宽度/mm	1
肋板宽度/mm	1
扩散层厚度/mm	0.3
多孔电极厚度/mm	1.29 × 10^−^2
质子交换膜厚度/mm	0.108
膜的电导率/(S·m^−^1)	9.825
多孔电极渗透率/m^2	2.36 × 10^−^12
阴极参考交换电流密度/(A·m^−^2)	0.001
阳极交换电流密度/(A·m^−^2)	100
扩散层渗透率/m^2	1.18 × 10^−^11
扩散层电导率/(S·m^−^1)	222

2.3. 数学模型

本研究开发的模型为三维、稳态、单相、非等温的PEMFC模型。为简化数值计算，采用如下假设[12]：

1) 气体流动为层流且是不可压缩的理想气体。

2) 多孔介质均各向同性。

3) 反应气体不能穿透PEM。

4) 忽略重力作用。

基于以上假设，控制方程可以表述如下[13]：

质量守恒方程：

$\nabla \cdot \left(\rho u\right)=S_m$ (1)

式中：ρ为流体密度，kg/m³；u为速度矢量，m/s；S_m为质量源项。

动量守恒方程：

$\nabla \cdot \left(\epsilon \rho uu\right)=-\epsilon \nabla P+\nabla \cdot \left(\epsilon {\mu }_{eff}\nabla u\right)+S_u$ (2)

式中：ε为孔隙率；P为流体压力，Pa；μ_eff为有效动力粘度，Pa·s；S_u为动量源项。

能量守恒方程：

$\rho C_{p}u\cdot \nabla T=\nabla \cdot \left(k_{eff}\nabla T\right)+S_e$ (3)

式中：C_p为气体定压比热容，J/kg·K；T为温度，K；k_eff为有效导热系数，W/m·K；S_e为能量源项。

组分守恒方程：

$\rho u\cdot \nabla {\omega }_k=\nabla \cdot \left(\rho D_k^{eff}\nabla {\omega }_k\right)+S_k$ (5)

式中：ω_k为不同组分的质量分数；$D_k^{eff}$ 为不同组分有效扩散系数；S_k为组分源项。

电荷守恒方程：

$\nabla \cdot \left(-{\sigma }_s\nabla \cdot {\varphi }_s\right)=S_{sol}$ (6)

$\nabla \cdot \left(-{\sigma }_m\nabla \cdot {\varphi }_m\right)=S_{mem}$ (7)

式中：σ_s、σ_m为分别为固相和膜的电导率，S/m；ϕ_s、ϕ_m为固相和膜相电势，V；S_sol、S_mem为电子和质子电流源项；

Butler-Volmer方程：

$i=i_0{\left(\frac{C_k}{C_{k,ref}}\right)}^{{\gamma }_i}\left[\exp \left(\frac{{\alpha }_{R}F}{RT}\eta \right)-\exp \left(\frac{-{\alpha }_{0}F}{RT}\eta \right)\right]$ (8)

式中：i为电流密度，A/m²；i₀为交换电流密度，A/m²；α_R、α₀为传递系数；γ为浓度修正系数；C_k、C_k,ref为反应物浓度和参考反应物浓度；T为工作温度，K；η为活化过电势，V；F为法拉第常数，96,485 C/mol。

2.4. 边界条件

除阳极和阴极的入口和出口外，所有边界均采用零流量边界条件。入口处的温度、压力和相对湿度可用于计算各组分的质量分数。阳极和阴极的入口速度由以下公式提供：

$u_{in,a}={\xi }_a\frac{i_0}{2F}{A}_m\frac{1}{{\omega }_{H_2}^{in}}\frac{R{T}_a^{in}}{p_a^{in}}\frac{1}{A_{ch}}$ (9)

$u_{in,c}={\xi }_c\frac{i_0}{4F}{A}_m\frac{1}{{\omega }_{O_2}^{in}}\frac{R{T}_c^{in}}{p_c^{in}}\frac{1}{A_{ch}}$ (10)

式中：ξ_a、ξ_c分别为阳极和阴极的化学计量比；i₀为参考工作电流密度，A/m²；A_m、A_ch分别为活化区域面积和流道入口截面积，m²；ω为不同组分气体的质量分数。

在阳极和阴极的出口边界，边界条件设置为压力出口，出口压力为PEMFC的工作压力。

2.5. 网格无关性验证与模型验证

Figure 3. Grid-independent verification

图3. 网格无关性验证

本文选择了56,160、67,392、80,870、97,044和116,453五种不同的网格数来对PEMFC模型进行网格独立性验证。如图3所示，在0.4 V电压下，网格数为80,870和97,044的模型之间的电流密度相对误差在0.5%以内。因此，可以认为80,870的网格数模型具有足够的计算精度，并在后续模拟计算中使用该网格数。

Figure 4. Comparison of simulation results with experimental data

图4. 模拟结果与实验数据对比

极化曲线是评估PEMFC运行期间综合性能的关键标准。为了验证本文数值模型的准确性，根据Sezgin等[14]模型的几何条件和工作条件，对文献中的直流道进行对比验证，结果如图4所示。从图中可以看出，模拟数据与实验数据非常吻合，从而证实了上述数学模型的可靠性。

3. 结果与讨论

3.1. 多目标的选择和首次正交实验设计

Table 2. Selected 10 variables and their respective level values

表2. 所选10个变量及各自水平值

序号	参数	缩写	水平1	水平2	水平3
1	挡板高度/mm	H	0.2	0.4	0.6
2	挡板角度/˚	θ	30	45	60
3	挡板周期	N	2	4	6
4	扩散层孔隙率	ɛ	0.4	0.5	0.6
5	操作压力/atm	P	1	2	3
6	操作温度/K	T	333	343	353
7	阳极化学计量比	ζa	1.5	3.5	5.5
8	阴极化学计量比	ζc	2	4	6
9	阳极相对湿度/%	RHa	80	90	100
10	阴极相对湿度/%	RHc	80	90	100

选择四个常用的性能指标：1. 电流密度I；2. 阴极流道入口与出口之间的压降ΔP；3. GDL/CL界面处氧气不均匀度N；净功率J。在传统的优化研究中，显著变量的选择往往基于经验，具有一定的主观性。为了在最少的显著变量下获得较好的优化结果，采用方差分析方法在优化前进行显著变量识别。选择对目标有显著影响的变量，达到减小优化空间、提高优化效率的目的。针对梯形挡板流道模型选择10个常用的PEMFC参数进行分析，所选10个变量及各自水平值见表2。几何参数和操作参数都被考虑在内。几何参数为挡板高度、挡板侧壁角度、挡板周期数、扩散层孔隙率；操作参数为压力、温度、阴/阳极化学计量比、阴/阳极相对湿度。

为了进行方差分析，首先进行正交实验设计[15]。每个参数都有三个不同的水平数，如表2所示。需要强调的是，这10个参数的三个水平值的选择是基于PEMFC参数的一般范围和一些经验[16] [17]。因此，根据正交实验法，一个L₂₇(3¹⁰)正交阵列被选择，方案见表3。对于每行的工况组合，由于PEMFC通常在0.4 V电位左右达到最大功率密度，所以通过COMSOL模拟计算出0.4 V电位下对应的四个性能指标的值。因此一共需要进行27组模拟仿真，得出不同工况组合下对应的各个性能指标。计算结果如表3所示。

Table 3. Arrangement design of orthogonal table L₂₇(3¹⁰) and its simulation results

表3. 正交表L₂₇(3¹⁰)的安排设计及其模拟结果

序号	H	θ	N	ɛ	P	T	ζa	ζc	RHa	RHc	电流密度I/(W·cm^−^2)	压降P/Pa	氧气不均匀度N	净功率J/W
1	0.2	30˚	2	0.4	1	333	1.5	2	80	80	1.0727	2.2670	0.66	0.171631
2	0.2	30˚	2	0.4	2	343	2.5	4	90	90	1.1449	5.0318	0.33	0.183181
3	0.2	30˚	2	0.4	3	353	3.5	6	100	100	1.1596	8.0176	0.26	0.185529
4	0.2	45˚	4	0.5	1	333	1.5	4	90	90	1.1466	4.9899	0.25	0.183453
5	0.2	45˚	4	0.5	2	343	2.5	6	100	100	1.1636	7.9809	0.19	0.186169
6	0.2	45˚	4	0.5	3	353	3.5	2	80	80	1.1056	2.4889	0.51	0.176895
7	0.2	60˚	6	0.6	1	333	1.5	6	100	100	1.1603	7.9856	0.15	0.185641
8	0.2	60˚	6	0.6	2	343	2.5	2	80	80	1.1078	2.4826	0.46	0.177247
9	0.2	60˚	6	0.6	3	353	3.5	4	90	90	1.1570	5.4890	0.20	0.185117
10	0.4	30˚	4	0.6	1	343	3.5	2	90	100	1.1105	3.5188	0.46	0.177679
11	0.4	30˚	4	0.6	2	353	1.5	4	100	80	1.1552	7.8436	0.19	0.184827
12	0.4	30˚	4	0.6	3	333	2.5	6	80	90	1.1713	11.6720	0.15	0.187398
13	0.4	45˚	6	0.4	1	343	3.5	4	100	80	1.1328	8.5660	0.30	0.181243
14	0.4	45˚	6	0.4	2	353	1.5	6	80	90	1.1458	14.3200	0.24	0.183315
15	0.4	45˚	6	0.4	3	333	2.5	2	90	100	1.0627	3.7275	0.62	0.170031
16	0.4	60˚	2	0.5	1	343	3.5	6	80	90	1.1784	14.9900	0.18	0.188531
17	0.4	60˚	2	0.5	2	353	1.5	2	90	100	1.1114	4.5586	0.50	0.177823
18	0.4	60˚	2	0.5	3	333	2.5	4	100	80	1.1617	9.2141	0.24	0.185867
19	0.6	30˚	6	0.5	1	353	2.5	2	100	90	1.1303	6.1930	0.55	0.180846
20	0.6	30˚	6	0.5	2	333	3.5	4	80	100	1.1883	14.1060	0.27	0.190120
21	0.6	30˚	6	0.5	3	343	1.5	6	90	80	1.2027	25.5730	0.19	0.192410
22	0.6	45˚	2	0.6	1	353	2.5	4	80	100	1.1719	23.7330	0.19	0.187490

续表

23	0.6	45˚	2	0.6	2	333	3.5	6	90	80	1.1867	35.1280	0.15	0.189843
24	0.6	45˚	2	0.6	3	343	1.5	2	100	90	1.1157	10.5000	0.45	0.178509
25	0.6	60˚	4	0.4	1	353	2.5	6	90	80	1.1752	39.3790	0.24	0.187997
26	0.6	60˚	4	0.4	2	333	3.5	2	100	90	1.0781	9.6117	0.64	0.172493
27	0.6	60˚	4	0.4	3	343	1.5	4	80	100	1.1508	22.9710	0.31	0.184115

3.2. 显著变量识别

采用方差分析方法判断各变量对结果的显著性[18]。方差分析是一种统计方法，可用于检验变量对实验结果的显著性。在这项工作中，方差分析的计算过程是通过MINITAB软件实现的。根据表3的模拟结果，分别对4个性能指标进行方差分析，统计结果见表4。

Table 4. ANOVA results for the four performance indicators

表4. 四个性能指标的方差分析结果

电流密度	变量	离差平方和	自由度	均方	F值	P值	贡献率	排名
	挡板高度	0.002296	2	0.001148	23.38	0.001	6.65%	3
	挡板侧壁角度	0.000603	2	0.000301	6.14	0.035	1.75%	4
	挡板周期数	0.000123	2	0.000061	1.25	0.352	0.36%
	扩散层孔隙率	0.004414	2	0.002207	44.96	0	12.78%	2
	操作压力	0.000004	2	0.000002	0.04	0.96	0.01%
	操作温度	0.00049	2	0.000245	4.99	0.053	1.42%
	阳极化学计量比	0.000079	2	0.00004	0.81	0.49	0.23%
	阴极化学计量比	0.02606	2	0.01303	265.42	0	75.47%	1
	阳极相对湿度	0.000108	2	0.000054	1.1	0.393	0.31%
	阴极相对湿度	0.00006	2	0.00003	0.61	0.574	0.17%
	误差	0.000295	6	0.000049	/	/	0.85%
	总和	0.03453	26	/	/	/	100.00%
压降	变量	离差平方和	自由度	均方	F值	P值	贡献率	排名
	挡板高度	1206.18	2	603.09	515.17	0	48.66%	1
	挡板侧壁角度	149.1	2	74.552	63.68	0	6.01%	4
	挡板周期数	41.42	2	20.71	17.69	0.003	1.67%
	扩散层孔隙率	34.46	2	17.229	14.72	0.005	1.39%
	操作压力	9.51	2	4.755	4.06	0.077	0.38%
	操作温度	10.9	2	5.449	4.65	0.06	0.44%
	阳极化学计量比	4.73	2	2.365	2.02	0.213	0.19%
	阴极化学计量比	796.76	2	398.381	340.31	0	32.14%	2
	阳极相对湿度	151.28	2	75.641	64.61	0	6.10%	3
	阴极相对湿度	67.47	2	33.733	28.82	0.001	2.72%
	误差	7.02	6	1.171	/	/	0.28%
	总和	2478.83	26	/	/	/	100.00%

续表

氧气不均匀度	变量	离差平方和	自由度	均方	F值	P值	贡献率	排名
	挡板高度	0.001089	2	0.000544	1.41	0.314	0.16%
	挡板侧壁角度	0.001689	2	0.000844	2.19	0.193	0.24%	4
	挡板周期数	0.000089	2	0.000044	0.12	0.893	0.01%
	扩散层孔隙率	0.081067	2	0.040533	105.23	0	11.54%	2
	操作压力	0.000156	2	0.000078	0.2	0.822	0.02%
	操作温度	0.004822	2	0.002411	6.26	0.034	0.69%	3
	阳极化学计量比	0.000067	2	0.000033	0.09	0.918	0.01%
	阴极化学计量比	0.610956	2	0.305478	793.07	0	86.97%	1
	阳极相对湿度	0.000067	2	0.000033	0.09	0.918	0.01%
	阴极相对湿度	0.000156	2	0.000078	0.2	0.822	0.02%
	误差	0.002311	6	0.000385	/	/	0.33%
	总和	0.702467	26	/	/	/	100.00%
净功率	变量	离差平方和	自由度	均方	F值	P值	贡献率	排名
	挡板高度	0.000058	2	0.000029	23.26	0.001	6.58%	3
	挡板侧壁角度	0.000015	2	0.000008	6.16	0.035	1.70%	4
	挡板周期数	0.000003	2	0.000002	1.25	0.352	0.34%
	扩散层孔隙率	0.000113	2	0.000057	45.05	0	12.81%	2
	操作压力	0	2	0	0.04	0.959	0.00%
	操作温度	0.000013	2	0.000006	4.99	0.053	1.47%
	阳极化学计量比	0.000002	2	0.000001	0.8	0.49	0.23%
	阴极化学计量比	0.000666	2	0.000333	265.16	0	75.51%	1
	阳极相对湿度	0.000003	2	0.000001	1.08	0.397	0.34%
	阴极相对湿度	0.000002	2	0.000001	0.6	0.579	0.23%
	误差	0.000008	6	0.000001	/	/	0.91%
	总和	0.000882	26	/	/	/	100.00%

在表4中，离差平方和是一组数据中每个数据点与该组数据的均值之差的平方的总和，均方是离差平方和除以自由度的结果，F的高值和P的低值表示变量更显著。在电流密度方差分析中，按照各个参数对电流密度的贡献度进行排名，排名前四的分别是阴极化学计量比、扩散层孔隙率、挡板高度、挡板侧壁角度，四者对电流密度的总贡献率为96.65%；同理，在压降方差分析中，排名前四的分别是挡板高度、阴极化学计量比、阳极相对湿度、挡板侧壁角度，四者对压降的总贡献率为92.91%；在氧气不均匀度方差分析中，排名前四的分别是阴极化学计量比、扩散层孔隙率、操作温度、挡板侧壁角度，四者对氧气不均匀度的总贡献率为99.44%；在净功率方差分析中，排名前四的分别是阴极化学计量比、扩散层孔隙率、挡板高度、挡板侧壁角度，四者对净功率的总贡献率为96.60%；因此，选择阴极化学计量比、挡板高度、挡板侧壁角度、扩散层孔隙率、操作温度、阳极相对湿度作为后续工作的显著变量。这些变量的贡献如图5所示。对于其他对性能指标没有相对显著影响的变量，其值固定在一定的常数，以降低后续优化模型的复杂度。即操作压力2 atm，阳极化学计量比3.5，挡板周期数4，阴极相对湿度100%。通过本小节，优化空间中的变量数量从10个减少到6个。

Figure 5. Contribution of decision variables to the four performance indicators

图5. 决策变量对四个性能指标的贡献度

3.3. 二次正交实验设计和数据获取

第二次正交试验旨在获取优化空间中具有全局代表性的点，以供后续ANN模型的训练和测试。为了增加数据的多样性和代表性，对六个显著变量各自选取了5个水平值进行正交实验设计。具体的水平值和正交表见表5和表6。每种工况组合下，本文计算了12个不同运行电压(0.4~0.95 V，间隔0.05取值)下的四个性能指标值，共得到25 × 12 = 300组数据。

Table 5. Six variables and their levels in the second orthogonal experiment

表5. 第二次正交实验的六个变量及其水平

	变量	水平1	水平2	水平3	水平4	水平5
1	ɛ	0.2	0.3	0.4	0.5	0.6
2	ζc	1.5	2.5	3.5	4.5	5.5
3	θ	30˚	45˚	60˚	75˚	90˚
4	H/mm	0.2	0.3	0.4	0.5	0.6
5	RHa	80	85	90	95	100
6	T/K	333	338	343	348	353

Table 6. L₂₅(5⁶) orthogonal table with six variables and five levels

表6. 具有六个变量及五个水平的L₂₅(5⁶)正交表

序号	ɛ	θ	ζc	H	RHa	T
1	0.2	30˚	1.5	0.2	80	333
2	0.2	45˚	2.5	0.3	85	338
3	0.2	60˚	3.5	0.4	90	343
4	0.2	75˚	4.5	0.5	95	348

续表

5	0.2	90˚	5.5	0.6	100	353
6	0.3	30˚	2.5	0.4	95	353
7	0.3	45˚	3.5	0.5	100	333
8	0.3	60˚	4.5	0.6	80	338
9	0.3	75˚	5.5	0.2	85	343
10	0.3	90˚	1.5	0.3	90	348
11	0.4	30˚	3.5	0.6	85	348
12	0.4	45˚	4.5	0.2	90	353
13	0.4	60˚	5.5	0.3	95	333
14	0.4	75˚	1.5	0.4	100	338
15	0.4	90˚	2.5	0.5	80	343
16	0.5	30˚	4.5	0.3	100	343
17	0.5	45˚	5.5	0.4	80	348
18	0.5	60˚	1.5	0.5	85	353
19	0.5	75˚	2.5	0.6	90	333
20	0.5	90˚	3.5	0.2	95	338
21	0.6	30˚	5.5	0.5	90	338
22	0.6	45˚	1.5	0.6	95	343
23	0.6	60˚	2.5	0.2	100	348
24	0.6	75˚	3.5	0.3	80	353
25	0.6	90˚	4.5	0.4	85	333

3.4. 熵值法获取PEMFC综合性能表达式CP

熵值法是一种定量权重分配法，它根据各指标数值变化对整体的影响，计算指标的熵值，进而确定权重[19]。具体计算步骤如下：

步骤1：使用最大(最小)归一化法对正交实验所得数据进行标准化处理。

${x^{\prime }}_{ij}=\frac{x_{ij}-{\left(x_{ij}\right)}_{\min }^j}{{\left(x_{ij}\right)}_{\max }^j-{\left(x_{ij}\right)}_{\min }^j} \left(正向目�\right)$ (11)

${x^{\prime }}_{ij}=\frac{{\left(x_{ij}\right)}_{\max }^j-x_{ij}}{{\left(x_{ij}\right)}_{\max }^j-{\left(x_{ij}\right)}_{\min }^j} \left(�向目�\right)$ (12)

步骤2：确定每个指标的熵值E_j。

$E_j=-K\cdot {\displaystyle \sum }_{i=1}^n{p}_{ij}\ln \left(p_{ij}\right)$ (13)

$p_{ij}=\frac{{x^{\prime }}_{ij}}{{{\displaystyle \sum }}_{i=1}^n{x^{\prime }}_{ij}}$ (14)

$K=\frac{1}{\ln \left(n\right)}$ (15)

步骤3：计算每个指标的权重W_j。

$W_j=\frac{1-E_j}{{{\displaystyle \sum }}_{j=1}^m\left(1-E_j\right)}$ (16)

步骤4：计算每个样本的综合得分Z_j

$Z_i={\displaystyle \sum }_{j=1}^m{W}_j{x^{\prime }}_{ij}$ (17)

式中：x_ij–第i个样本的第j个指标值，$i=1,2,3,\cdots ,n$ ；$j=1,2,3,\cdots ,m$ ；E_j–每个指标的熵值；W_j–每个指标的权重；Z_j–每个样本的综合得分；n–样本数；如果$p_{ij}=0$ ，则定义$\underset{p_{ij}\to 0}{\lim }{\left(1+\frac{1}{n}\right)}^n{p}_{ij}\ln \left(p_{ij}\right)=0$ 。

根据EWM方法以及COMSOL模拟计算所获取到的300组模拟数据，计算出了电流密度、压降、氧气不均匀度和净功率四个性能指标对应的权重，如图6所示。因此，综合性能表达式CP的数学表达式为$\text{CP}=0.5155I+0.3655J+0.0751P+0.0439N$ 。

Figure 6. Contribution of significant variables to the four performance indicators

图6. 显著变量对四个性能指标的贡献度

3.5. ANN模型训练与测试

人工神经网络(ANN)模型可以根据得到的输入和输出信息对非线性函数进行处理，并很好地映射输入数据和输出数据。ANN的优点是利用经验数据具有良好的泛化能力，计算速度非常快。因此，本研究使用ANN来解决回归问题。ANN包含三部分：输入层、隐藏层和输出层。通过隐藏层神经元的多种组合，得到从输入到输出的映射关系：

$y_j=f_{\text{act}}\left({\displaystyle {\sum }_{i=1}^n{x}_i{\omega }_{ij}}-{\theta }_j\right)$ (18)

其中x_i表示输入变量；y_j表示输出变量；ω_ij为输入层与隐藏层间的连接权值；θ_j表示阈值；f_act表示非线性函数的激活函数。图7显示了ANN的原理图。

Figure 7. ANN model structure diagram

图7. ANN模型结构图

该模型的七个输入参数是ɛ、θ、ζ_c、H、RH_a、T和p，输出参数是综合性能指标CP。采用试错法选择隐藏层数，模型中的隐藏层数为2，每个隐藏层的节点数分别为20和10。ANN模型使用ReLU作为激活函数，ANN模型由Python中的scikit-learn库实现。在建立ANN模型后，按照5:1的比例把300组数据随机分为训练集和测试集，其中训练集包含250组数据，测试集包含50组数据。采用训练集进行训练，通过测试集评估预测性能。

决定系数(R²)用于测量ANN模型的预测精度，其表示为：

${\text{R}}^2=1-\frac{{\displaystyle {\sum }_{i=1}^n{\left(y-y_{\text{pre}}\right)}^2}}{{\displaystyle {\sum }_{i=1}^n{\left(y-y_{\text{ave}}\right)}^2}}$ (19)

其中n表示样本数量，y，y_pre表示第i个样本的真实值和预测值，y_ave表示所有样本真实值的平均值。R²的值越接近1表示预测能力越准确。如图8所示，模型经过训练后，得到了很高的拟合度，在训练过程中，R²达到了0.99833，在测试过程中，R²达到了0.99741，由此表明ANN模型的拟合结果良好。

Figure 8. Correlation between ANN model predictions and actual values

图8. ANN模型预测值与实际值之间的相关性

3.6. 多目标优化结果

根据训练好的ANN模型，对扩散层孔隙率($0.2,0.25,\cdots ，0.55,0.6$ )，挡板侧壁角度 (${30}^{\circ },{31}^{\circ },\cdots ,{89}^{\circ },{90}^{\circ }$ )，阴极化学计量比($1.5,2.0,\cdots ,5.0,5.5$ )，挡板高度 ($0.2\text{ mm},0.25\text{ mm},\cdots ,0.55\text{ mm},0.6\text{ mm}$ )，阳极相对湿度($80\%,81\%,\cdots ,99\%,100\%$ )，操作温度 ($333\text{ K},332\text{ K},\cdots ,352\text{ K},353\text{ K}$ )共9 × 61 × 9 × 9 × 21 × 21 = 19,610,829组数据进行预测，为了能获得最优的PEMFC性能，本文对预测CP值进行最优选择。最终，当扩散层孔隙率、挡板侧壁角度、阴极化学计量比、挡板高度、阳极相对湿度、操作温度、运行电压分别为0.57、51˚、5.0、0.39 mm、99%、333 K、0.4 V时，CP达到最大值0.9622，PEMFC达到最优性能。图8显示了优化流道模型与基础流道模型(0.2、30˚、1.5、0.2 mm、80%、333 K)的性能比较。如图9(a)和图9(b)所示，在0.4 V电压下，优化后流道模型相比较基础流道而言，功率密度和净功率分别提升了53.71%和53.70%，然而阴极压降却仅仅只上升了7.93 Pa。并且从图9(c)中可看出，优化后流道模型在阴极GDL/CL界面处的氧气不均匀度明显小于基础流道模型，经计算，基础流道模型的氧气不均匀度1.59，优化后流道模型的氧气不均匀度降至0.19，氧气分布的更加均匀，使得氧气更容易扩散至催化层表面参与反应，进而提升了PEMFC性能。

Figure 9. Comparison of the performance of the optimized model with the base model: (a) Polarization curves; (b) Net power and cathode pressure drop; (c) Oxygen mole fraction

图9. 优化模型与基础模型的性能比较：(a) 极化曲线；(b) 净功率和阴极压降；(c) 氧气摩尔分数

此外，在Intel(R) Xeon(R) Gold 6240 CPU @ 2.60GHz计算机上通过Python运行此ANN模型的总时间为15分钟46秒。相比之下，通过CFD仿真获得一组结果大约需要25分钟。如果像以前的文献一样优化算法计算目标函数，而不是使用代理模型，则该算法可能需要近14天的时间才能迭代800代。因此，所提出的综合性能表达式CP与训练好的ANN模型相结合，可以显著减少优化时间，同时确保不影响优化的准确性。

4. 结论

本文提出一种快速、系统的多目标优化方法，利用该方法对具有梯形挡板的PEMFC流道模型进行了多目标优化分析，全面考虑了其几何参数和操作参数，同时优化4个PEMFC性能指标电流密度、净功率、压降和氧气不均匀度，最大限度地提升了考虑多目标的PEMFC的综合性能。得出以下结论：

1) 通过将ANN模型与正交实验、方差分析法和熵值法结合的多目标优化方法仅用15分钟46秒便得出了优化结果，与传统的优化算法相比，可以显著缩短优化时间。

2) 针对PEMFC的几何参数和操作参数通过正交实验和方差分析，确定扩散层孔隙率、阴极化学计量比、挡板侧壁角度、挡板高度、阳极相对湿度和操作温度这6个参数对PEMFC性能的影响显著。因此将它们作为后续优化的决策变量，将优化空间的维度从10降至6。

3) 由熵值法确定PEMFC综合性能表达式$\text{CP}=0.5155I+0.3655J+0.0751P+0.0439N$ ，根据所获得的数据训练并验证ANN模型，训练和验证的R²分别高达0.99833和0.99741，由此看出ANN模型体现了较好的预测精度。

4) 用验证过的ANN模型在给定输入参数范围条件下对综合性能表达式CP进行预测，并对预测CP值进行最优选择，当扩散层孔隙率、挡板侧壁角度、阴极化学计量比、挡板高度、阳极相对湿度、操作温度、运行电压分别为0.57、51˚、5.0、0.39 mm、99%、333 K、0.4 V时，CP值达到最大值0.9622，PEMFC达到最佳性能。

5) 通过将优化后流道模型与基础模型的性能对比可以得出：在0.4 V电压下，优化后流道模型的功率密度相比较基础流道而言，功率密度和净功率分别提升了53.71%和53.70%，然而阴极压降却仅仅只上升了7.93 Pa，氧气不均匀度降至0.19。

NOTES

^*通讯作者。