ABSTRACT

Proton Exchange Membrane Fuel Cell (PEMFC) is a complex multi-variable, strongly coupled and non-linear system. Accurate modeling method is the basis of its research, and advanced control strategy is the key to improving the power generation performance. In recent years, many studies have shown that the dynamic processes of PEMFC power generation, such as gas diffusion, heat conduction and electrochemical reaction, have fractional order characteristics. In this paper, a fractional order state space model of PEMFC is established by combining fractional order theory with subspace identification method (SIM).

Keywords:PEMFC, Fractional Order Theory, Spatial Identification Method

PEMFC分数阶子空间建模

苏申，戈未平

南京理工大学自动化学院，江苏 南京

收稿日期：2019年8月7日；录用日期：2019年8月23日；发布日期：2019年8月30日

摘 要

质子交换膜燃料电池(Proton Exchange Membrane Fuel Cell, PEMFC)是一个复杂的多变量、强耦合、非线性系统，精确的建模方法是对其进行研究的基础，而先进的控制策略则是提高发电性能的关键。近年来，不少研究表明PEMFC发电过程中的气体扩散、热量传导以及电化学反应等动态过程存在分数阶特性。为此，本文将分数阶理论与子空间辨识方法(Subspace Identification Method, SIM)相结合，建立PEMFC的分数阶状态空间模型。

关键词 :PEMFC，分数阶理论，空间辨识方法

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

PEMFC是一个包括流动、传质、传热和电化学反应等多种物理化学现象的复杂机体，若要对其进行深入研究，建模是一种直观且快速的手段。根据建模方法划分，PEMFC模型可分为机理模型和经验模型。机理模型主要是针对电堆内部的气体、质子、水以及热量分布等进行建模，采用能量守恒方程、传质传热方程和电化学反应等方程来描述相应参数对PEMFC的输出特性的影响。Amphlett [1] 提出了PEMFC的电化学模型，描述了阴极/阳极气体分压、温度和电流密度对输出电压的影响。朱柳 [2] 根据能量守衡定律，考虑了电堆自身的热辐射及欧姆极化的影响，构建了PEMFC的热管理模型。Zhang [3] 研究了PEMFC电堆内部水的气、液两相分布问题，根据质量守恒定律，建立了PEMFC的水管理模型。此外，还有不少学者还考虑了反应物浓度、电堆温度、水、电流密度等因素在二维或者三维空间分布不均对PEMFC的影响，建立了更加精确和全面的PEMFC模型，但相应的模型也更加复杂 [4] [5] [6] [7] [8] 。机理建模需对PEMFC系统有大量的先验知识，理解相对复杂，但可以反映电池内部参数变化对输出特性的影响。机理模型建模过程复杂、模型参数众多且获取难度大，适合对电池性能进行分析，一般不适用于控制。

近年来已有研究表明，PEMFC是一个具有分数阶特性的系统，采用整数阶方法对其进行建模往往容易忽略系统的真实性。为此，Bian [9] 和Qi [10] 分别采用分数阶微分方程描述PEMFC的气体扩散和热量耗散过程，仿真表明分数阶模型能更加精确的刻画PEMFC的动态过程。张明 [11] 将分数阶电容理论引入到PEMFC等效电路模型当中，建立了PEMFC的分数阶等效电路模型。Table [12] 对PEMFC电堆进行电化学阻抗谱测试，实验表明结果表明采用分数阶等效电路模型能更为准确地描述PEMFC的阻抗特性。在针对PEMFC分数阶系统辨识方面，胡聪 [13] 将分数阶子空间辨识方法运用了PEMFC系统辨识当中，并建立了空冷型PEMFC的分数阶状态空间模型。文中作者将PEMFC的辨识过程分为两步，忽略了系统系数矩阵和分数阶次之间的相互影响，且对于辨识算法中其他参数的选取过于经验化，没有考虑到全局参数的最优估计问题。

综上所述，目前针对PEMFC系统的分数阶特性的建模研究较少，且主要集中在机理模型和等效电路模型，有关于PEMFC的分数阶系统辨识建模的研究更是少之又少。因此，本文将在前人研究的基础上，基于子空间辨识方法，将分数阶理论运用于PEMFC系统辨识当中，以获得更好的辨识效果。

2. PEMFC概述

EMFC单电池结构如图1所示，主要由双极板、密封圈和膜电极组件(Membrane Electrodes Assembly，MEA)构成，其中MEA为PEMFC的核心部分，包括质子交换膜、阴/阳极催化层和气体扩散层。

PEMFC的工作原理实际上是水电解的逆过程，如图2所示。氢气经过气体流场注入阳极室，然后通过气体扩散层到达阳极催化层，在阳极催化剂的作用下，氢气分解为质子和电子，质子穿过交换膜到达

图1. PEMFC单电池结构示意图

阴极，而电子则通过外电路到达阴极，在阴极催化剂作用下和氧气反应生成水，尾气通过排气装置排出。电极反应方程式如下：

阳极反应： ${\text{H}}_2=2{\text{H}}^{+}+2{\text{e}}^{-}$ (1)

阴极反应： ${\text{2H}}^{\text{+}}\text{+}2{\text{e}}^{-}\text{+}\frac{1}{2}{\text{O}}_2\to {\text{H}}_2\text{O}$ (2)

总反应方程式： ${\text{H}}_2+\frac{1}{2}{\text{O}}_2\to {\text{H}}_2\text{O}+热量$ (3)

图2. PEMFC工作原理示意图

当外电路连接负载时，电堆输出电压为：

$E=N\cdot \left(E_{nernst}-V_{act}-V_{ohm}-V{}_{con}\right)$ (4)

式中，N表示电堆中包含的单电池数目， $E_{nernst}$ 为热力学电动势，即单电池开路电压， $V_{act}$ 、 $V_{ohm}$ 和 $V_{con}$ 分别表示活化电压损耗、欧姆电压损耗以及浓差电压损耗。

$E_{nernst}$ 可以由下式表示：

$E_{nernst}=\frac{\Delta G}{2F}+\frac{\Delta S}{2F}\left(T-T_{ref}\right)+\frac{RT}{2F}\left(\ln P_{H_2}+\frac{1}{2}\ln P_{O_2}\right)$ (5)

式中， $\Delta G$ 为吉布斯自由能， $\Delta S$ 为气体摩尔熵，R和F分别表示气体常数和法拉第常数，T为电堆工作温度， $T_{ref}$ 为参考温度， $P_{H_2}$ 和 $P_{O_2}$ 分别表示氢气和氧气分压。

$V_{act}$ 表示电池电化学反应克服活化对外做功时带来的损耗，可以下式表示：

$V_{act}={\eta }_1+{\eta }_{2}T+{\eta }_{3}T\ln C_{O_2}+{\eta }_{4}T\ln i$ (6)

式中，i为电堆电流， ${\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4$ 表示和电堆相关的经验系数， $C_{O_2}$ 表示阴极催化层与反应气体交界面的氧气浓度，可由亨利定律求得：

$C_{O_2}=\frac{P_{O_2}}{5.08\times {10}^6\times \exp \left(-498/T\right)}$ (7)

$V_{ohm}$ 主要是由单电池等效内阻 $R_c$ 引起，满足欧姆定律，可表示为：

$V_{ohm}=i\cdot R_c$ (8)

$R_c$ 是与质子交换膜电阻率、电池材料、结构、工作温度以及湿度等因素相关。

$V_{con}$ 是在电化学反应中，由于质量传输导致反应物浓度(氢气、氧气、水以及离子等)发生变化而产生的，其表达式如下：

$V_{con}=-B\times \ln \left(1-\frac{J}{J_{\max }}\right)$ (9)

式中，B为与电池工况相关的参数，J为电池平均电流密度， $J_{\max }$ 为最大电流密度。

由式(4)可知，电堆输出电压除了与本身的材料及结构等参数相关，同时还与其工作状态，如电流、温度、湿度、反应物浓度及压力等因素相关。可见，PEMFC机理建模过程十分复杂，且相关参数获取难度大。因此，亟需寻找新的适用于PEMFC的建模方法。

3. PEMFC分数阶子空间建模

分数阶微积分操作算子 ${}_a\mathcal{D}{}_t^{\alpha }$ 定义如下：

${}_a\mathcal{D}{}_t^{\alpha }=\{\begin{array}{l}\frac{{\text{d}}^{\alpha }}{\text{d}t},R\left(\alpha \right)>0\\ 1,R\left(\alpha \right)=0\\ {\displaystyle {\int }_a^t{\left(\text{d}\tau \right)}^{-\alpha }},R\left(\alpha \right)<0\end{array}$ (10)

式中， $\alpha$ 为微积分阶次，可以为任意的实数甚至复数， $R\left(\alpha \right)$ 为 $\alpha$ 的实部， $t,a$ 为操作算子的上下线。

整数阶系统类似，线性定常分数阶系统可分为以下三种：分数阶传递函数、状态空间方程以及微分方程 [10] 。分数阶微分方程可描述如下：

${\mathcal{D}}^{{\alpha }_n}y\left(t\right)+a_{n-1}{\mathcal{D}}^{{\alpha }_{n-1}}y\left(t\right)+\cdots +a_0{\mathcal{D}}^{{\alpha }_0}y\left(t\right)=b_m{\mathcal{D}}^{{\beta }_m}u\left(t\right)+b_{m-1}{\mathcal{D}}^{{\beta }_{m-1}}u\left(t\right)+\cdots b_0{\mathcal{D}}^{{\beta }_0}u\left(t\right)$ (11)

式中，分数阶阶次满足 ${\alpha }_n>{\alpha }_{n-1}>\cdots >{\alpha }_0$ 且 ${\beta }_m>{\beta }_{m-1}>\cdots >{\beta }_0$ 。当上式中所有分数阶阶次均为某一常数的倍数时，则称该系统为同元阶次分数阶系统。若不满足该条件，则称系统为分布阶次系统。

假设函数 $f\left(t\right)$ 在 $t=0$ 时刻初始值为0，可以得到 $D^{\alpha }f\left(t\right)$ 的Laplace变换：

$\mathcal{L}\left[{\mathcal{D}}^{\alpha }f\left(t\right)\right]=s^{\alpha }F\left(s\right)$ (12)

观察上式可知，其形式与整数阶Laplace变换一致，当 $\alpha$ 为整数时，分数阶Laplace变换则退化为整数阶的Laplace变换。参考式(12)，对式(11)的进行Laplace变换形式，可以得到分数阶传递函数形式如下：

$G\left(s\right)=\frac{Y\left(s\right)}{U\left(s\right)}=\frac{b_m{s}^{{\beta }_m}+b_{m-1}{s}^{{\beta }_{m-1}}+\cdots +b_0}{s^{{\alpha }_n}+a_{n-1}{s}^{{\alpha }_{n-1}}+\cdots +a_0}$ (13)

对于多输入多输出系统，可以采用分数阶状态空间方程来描述：

$\{\begin{array}{l}{\mathcal{D}}^{\alpha }x\left(t\right)=Ax\left(t\right)+Bu\left(t\right)\hfill \\ y\left(t\right)=Cx\left(t\right)+Du\left(t\right)\hfill \end{array}$ (14)

式中， $0<\alpha <2$ ，这里只考虑同元阶次的分数阶系统。 $u\left(t\right)$ 和 $y\left(t\right)$ 分别为输入和输出信号， $A,B,C,D$ 为系统系数矩阵。

分数阶状态空间模型描述如下：

$\{\begin{array}{l}{\mathcal{D}}^{\alpha }x\left(t\right)=Ax\left(t\right)+Bu\left(t\right)+w\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Du\left(t\right)+v\left(t\right)\end{array}$ (15)

式中， $\alpha$ 为分数阶微分阶次，且满足 $0<\alpha <2$ ，这里只考虑同元阶次的分数阶系统。 $u\left(t\right)\in N^m$ 和 $y\left(t\right)\in N^l$ 分别为可测输入和输出，m和l分别为输入和输出维数。 $w\left(t\right)\in N^n$ 和 $v\left(t\right)\in N^l$ 称为系统过程噪声和测量噪声。 $A\in N^{n\times n}$ ， $B\in N^{m\times n}$ ， $C\in N^{l\times n}$ ， $D\in N^{l\times m}$ 为系统系数矩阵，n为系统阶次。对式(15)中输出方程左右两边求 $\alpha$ 阶导数，可得：

$\begin{array}{c}{\mathcal{D}}^{\alpha }y\left(t\right)=C\cdot {\mathcal{D}}^{\alpha }x\left(t\right)+D\cdot {\mathcal{D}}^{\alpha }u\left(t\right)+{\mathcal{D}}^{\alpha }v\left(t\right)\\ =CAx\left(t\right)+\left(CB+D\cdot {\mathcal{D}}^{\alpha }\right)u\left(t\right)+Cw\left(t\right)+{\mathcal{D}}^{\alpha }v\left(t\right)\end{array}$ (16)

对上式进行迭代求导，有：

$\begin{array}{c}{\mathcal{D}}^{2\alpha }y\left(t\right)=CA\cdot {\mathcal{D}}^{\alpha }x\left(t\right)+\left(CB\cdot {\mathcal{D}}^{\alpha }+D\cdot {\mathcal{D}}^{2\alpha }\right)u\left(t\right)+C\cdot {\mathcal{D}}^{\alpha }w\left(t\right)+{\mathcal{D}}^{2\alpha }v\left(t\right)\\ =C{A}^{2}x\left(t\right)+\left(CAB+CB\cdot {\mathcal{D}}^{\alpha }+D\cdot {\mathcal{D}}^{2\alpha }\right)u\left(t\right)+\left(CA+C\cdot {\mathcal{D}}^{\alpha }\right)w\left(t\right)+{\mathcal{D}}^{2\alpha }v\left(t\right)\end{array}$ (17)

$⋮$

$\begin{array}{l}{\mathcal{D}}^{\left(q-1\right)\alpha }y\left(t\right)=C{A}^{q-1}x\left(t\right)+\left(C{A}^{q-2}B+C{A}^{q-2}B\cdot {\mathcal{D}}^{\alpha }+\cdots +D\cdot {\mathcal{D}}^{\left(q-1\right)\alpha }\right)u\left(t\right)\\ +\left(C{A}^{q-2}+C{A}^{q-3}\cdot {\mathcal{D}}^{\alpha }+\cdots +C\cdot {\mathcal{D}}^{\left(q-1\right)\alpha }\right)w\left(t\right)+{\mathcal{D}}^{\left(q-1\right)\alpha }v\left(t\right)\end{array}$ (18)

将式(16)~(18)改写为矩阵形式，有

$Y\left(t\right)=\Gamma x\left(t\right)+HU\left(t\right)+GW\left(t\right)+V\left(t\right)$ (19)

式中，q为正整数，且需满足 $q\ge n$ 。 $\Gamma$ 为广义能观矩阵，H和G为下三角Toeplitz矩阵，定义如下：

$\Gamma =\left[\begin{array}{c}C\\ CA\\ ⋮\\ C{A}^{q-1}\end{array}\right],U\left(t\right)=\left[\begin{array}{c}u\left(t\right)\\ {\mathcal{D}}^{\alpha }u\left(t\right)\\ ⋮\\ {\mathcal{D}}^{\left(q-1\right)\alpha }u\left(t\right)\end{array}\right],Y\left(t\right)=\left[\begin{array}{c}y\left(t\right)\\ {\mathcal{D}}^{\alpha }y\left(t\right)\\ ⋮\\ {\mathcal{D}}^{\left(q-1\right)\alpha }y\left(t\right)\end{array}\right]$

$H=\left[\begin{array}{cccc}D& 0& \cdots & 0\\ CB& D& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C{A}^{q-2}B& C{A}^{q-3}B& \cdots & D\end{array}\right],G=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ C& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ C{A}^{q-2}& C{A}^{q-3}& \cdots & 0\end{array}\right]$

矩阵 $W\left(t\right)$ 和 $V\left(t\right)$ 的构造形式与 $U\left(t\right)$ 类似。

选取合适的采样时间 $T_s$ ，令 $t_k=k\cdot T_s,k=1,2\cdots N$ , $u\left(k\right)=u\left(t_k\right)$ , $y\left(k\right)=y\left(t_k\right)$ ，由式(19)可以得到：

$Y_N\left(k\right)=\Gamma X_N\left(k\right)+H{U}_N\left(k\right)+G{W}_N\left(k\right)+V_N\left(k\right)$ (20)

式中

$X_N\left(k\right)=\left[\begin{array}{cccc}x\left(k\right)& x\left(k+2\right)& \cdots & x\left(N\right)\end{array}\right]$

$U_N\left(k\right)=\left[\begin{array}{cccc}U\left(k\right)& U\left(k+1\right)& \cdots & U\left(N\right)\end{array}\right]$

$Y_N\left(k\right)=\left[\begin{array}{cccc}Y\left(k\right)& Y\left(k+1\right)& \cdots & Y\left(N\right)\end{array}\right]$

矩阵 $W_N\left(k\right)$ 和 $V_N\left(k\right)$ 的构造形式与 $U_N\left(k\right)$ 类似。

由式(19)可知，可测输入和输出信号必须满足 $\left(q-1\right)\alpha$ 阶分数阶可导。但在实际辨识中输入输出数据很难满足这样的条件，导致后续辨识过程无法进行。为解决此问题，本文引入Poisson滤波器对输入输出数据进行滤波处理。Poisson矩函数描述如下：

$G_k\left(s\right)={\left(\frac{\beta }{s+\lambda }\right)}^k$ (21)

式中，s为Laplace算子，k为滤波器阶次，可以理解为k个低通滤波器 $\beta /\left(s+\lambda \right)$ 的级联。 $g_k\left(t\right)$ 为 $G_k\left(s\right)$ 的时域形式，表达式如下：

$g_k\left(t\right)={\mathcal{L}}^{-1}\left(G_k\left(s\right)\right)={\beta }^k\frac{t^{k-1}}{\left(k-1\right)!}{\text{e}}^{-\lambda t}$ (22)

对运算符 $M\{\cdot \}$ 作如下定义：

$M\left\{f\left(t\right)\right\}=g_k\left(t\right)\ast f\left(t\right)$ (23)

式中， $\ast$ 为卷积运算符， $f\left(t\right)$ 为关于时间t的函数，根据卷积运算法则和分数阶理论基础有：

$M\left\{{\mathcal{D}}^{\alpha }f\left(t\right)\right\}={\mathcal{L}}^{-1}\left\{G_k\left(s\right)s^{\alpha }F\left(s\right)\right\}={\mathcal{L}}^{-1}\left\{s^{\alpha }{G}_k\left(s\right)F\left(s\right)\right\}={\mathcal{D}}^{\alpha }{g}_k\left(t\right)\ast f\left(t\right)$ (24)

上式可以解释为：通过使用运算符 $M\{\cdot \}$ 将原本对函数 $f\left(t\right)$ 的 $\alpha$ 阶微分转移到对Poisson矩函数 $g_k\left(t\right)$ 的 $\alpha$ 阶微分。如此，采用上述方法对输入输出信号进行滤波处理，可以把原本对输入输出数据的各阶分数阶微分转化为对的Poisson矩函数的各阶分数阶微分。对输入数据进行滤波处理可以得到：

$M\left\{U\left(k\right)\right\}=\left[\begin{array}{c}M\left\{u\left(k\right)\right\}\\ M\left\{{\mathcal{D}}^{\alpha }u\left(k\right)\right\}\\ ⋮\\ M\left\{{\mathcal{D}}^{\left(q-1\right)\alpha }u\left(k\right)\right\}\end{array}\right]=\left[\begin{array}{c}{g}_k\left(k\right)\ast u\left(k\right)\\ g_k\left(k\right)\ast {\mathcal{D}}^{\alpha }u\left(k\right)\\ ⋮\\ g_k\left(k\right)\ast {\mathcal{D}}^{\left(q-1\right)\alpha }u\left(k\right)\end{array}\right]=\left[\begin{array}{c}{g}_k\left(k\right)\ast u\left(k\right)\\ {\mathcal{D}}^{\alpha }{g}_k\left(k\right)\ast u\left(k\right)\\ ⋮\\ {\mathcal{D}}^{\left(q-1\right)\alpha }{g}_k\left(k\right)\ast u\left(k\right)\end{array}\right]$ (25)

由式(19)可知，令 $k\ge \left(q-1\right)\cdot \alpha$ ，则Poisson矩函数 $g_k\left(t\right)$ 是可以满足 $\left(q-1\right)\cdot \alpha$ 可导的。对式(20)进行整体滤波可得：

$M\left\{Y_N\left(k\right)\right\}=\Gamma M\left\{X_N\left(k\right)\right\}+HM\left\{U_N\left(k\right)\right\}+GM\left\{W_N\left(k\right)\right\}+M\left\{V_N\left(k\right)\right\}$ (26)

式中

$M\left\{U_N\left(k\right)\right\}=\left[\begin{array}{cccc}M\left\{U\left(k\right)\right\}& M\left\{U\left(k+1\right)\right\}& \cdots & M\left\{U\left(N\right)\right\}\end{array}\right]$

矩阵 $M\left\{Y_N\left(k\right)\right\}$ 和 $M\left\{X_N\left(k\right)\right\}$ ， $M\left\{W_N\left(k\right)\right\}$ 和 $M\left\{V_N\left(k\right)\right\}$ 的构造方式与 $M\left\{U_N\left(k\right)\right\}$ 类似。

分数阶子空间辨识算法流程图如图3所示：

图3. 分数阶子空间辨识算法流程图

4. 仿真及结果

在25˚C室温下，通过图4中所示PEMFC实验平台，每隔0.5 s采集共750组电堆正常工作状态下的数据，其中前500组用于模型训练，后250组用于模型验证。电堆电流变化周期为10 s，氢气流量变化周期为5 s。输入信号变化曲线如图5和图6所示：

图4. 空冷自增湿型PEMFC实验平台

图5. 氢气流量变化曲线图

图6. 电堆电流变化曲线图

在优化过程中令Hankel矩阵行数 $q=10$ 、系统阶次 $n=3$ 、离散时间 $T_s=0.5$ 。以式(24)为目标函数，采用Matlab自带遗传算法工具箱对辨识过程进行优化，得到最优参数值如下： $\alpha =0.947$ , $L=83$ , $k=10$ , $\lambda =14.6$ , $\beta =11.9$ 。图7和图8为辨识结果与实测数据的对比图：

从上图可以看出，辨识结果可以很好地跟随实测数据的变化。尤其是输出电压，基本和实测数据保持吻合，从而证明了分数阶子空间辨识算法的有效性。利用该辨识算法辨识得到PEMFC系统的分数阶状态空间模型如下：

图7. 分数阶子空间辨识输出电压结果对比图

图8. 分数阶子空间辨识输出功率结果对比图

$\alpha =0.947$

$A=\left[\begin{array}{ccc}0.9697& 0.2059& 0.0119\\ 0.0431& 0.5691& 0.5883\\ -0.0912& 0.2284& 0.6630\end{array}\right]$ , $B=\left[\begin{array}{cc}-0.0496& -0.0332\\ -1.0696& 0.1410\\ 0.6745& -0.0822\end{array}\right]$

$C=\left[\begin{array}{ccc}-0.1522& 0.8536& -0.9611\\ -1.4070& -1.5883& -0.5296\end{array}\right]$ , $D=\left[\begin{array}{cc}-2.1070& 0.0013\\ 14.8263& 0.0021\end{array}\right]$

辨识效果验证将从辨识模型精度和辨识算法运行时间两方面来进行。

在辨识模型精度方面，将整数阶子空间辨识结果与分数阶子空间辨识结果作对比，两种辨算法所得模型预测输出效果对比如图9和图10所示：

图9. 不同辨识算法模型预测输出电压与实测数据对比图

图10. 不同辨识算法模型预测输出功率与实测数据对比图

5. 总结

本文主要研究了PEMFC发电过程的分数阶状态空间模型辨识。首先，介绍分数阶微积分的基本定义，为本文的分数阶辨识算法提供预备知识；其次，研究并推导了分数阶子空间辨识的基本原理。采用Poisson滤波器对输入输出数据进行滤波，将对输入输出数据的分数阶微分转化为对Poisson矩函数的微分；引入短时记忆法，将分数阶微分的迭代求解过程转化为对短时记忆矩阵和增广输入输出矩阵的乘积运算，减少算法的计算量；利用遗传算法对辨识算法中存在的未知参数进行全局寻优，以获得更好的辨识效果。最后，结合本文研究对象，建立PEMFC系统的分数阶状态空间模型。仿真及实验结果表明，所提出的辨识算法能有效地提升建模精度。

文章引用

苏 申,戈未平. PEMFC分数阶子空间建模
Research on Fractional Order Subspace Modeling[J]. 建模与仿真, 2019, 08(03): 125-135. https://doi.org/10.12677/MOS.2019.83015

参考文献