CN107632522A - 一种质子交换膜燃料电池非线性状态空间模型辨识方法 - Google Patents
一种质子交换膜燃料电池非线性状态空间模型辨识方法 Download PDFInfo
- Publication number
- CN107632522A CN107632522A CN201710772967.5A CN201710772967A CN107632522A CN 107632522 A CN107632522 A CN 107632522A CN 201710772967 A CN201710772967 A CN 201710772967A CN 107632522 A CN107632522 A CN 107632522A
- Authority
- CN
- China
- Prior art keywords
- mrow
- mtd
- msub
- mtr
- msubsup
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Fuel Cell (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
Abstract
本发明公开了一种质子交换膜燃料电池非线性状态空间模型辨识方法。首先选择氢气流量和负载电流作为输入变量,选择电压作为输出变量,采集大量的数据。其次,通过采集的数据构造Hankel矩阵,求解对偶矩阵。然后,利用矩阵方程的结果构造矩阵投影。最后,对斜投影利用奇异值分解得到系统的状态序列估计。重复使用上述方法求解系统矩阵,并利用低秩近似技术获得系统非线性特性估计。以实测输入数据为横坐标,辨识非线性特性作为纵坐标利用MATLAB曲线拟合工具箱对系统非线性特性进行多项式拟合。该方法可以对质子交换膜燃料电池的非线性特性进行十分准确的描述,不仅对实际工程人员对质子交换膜燃料电池建模以及后续控制系统设计提供解决方案,而且对此类系统建模均有较好的参考价值。
Description
技术领域
本发明属于工业控制领域,具体是一种质子交换膜燃料电池非线性状态空间模型辨识方法。
背景技术
能源短缺、环境污染形势严峻的当下,发展新型清洁、可持续能源是世界任何国家的必走之路。质子交换膜燃料电池是一种新型发电装置,仅需要通入氢气和空气即可产生电能,拥有高效、清洁、噪声小等突出优点,在新能源领域扮演着越来越重要的角色。燃料电池堆的设计研发、以及发电控制系统的设计、仿真,需要精确的数学模型。
现有常用的燃料电池建模方法主要集中在机理建模,通过能量守恒定律、电化学反应方程等物理化学原理对电堆的静态/动态、电特性/温度进行建模。机理建模往往能取得不错的效果,但是依然存在以下几个主要问题:
1)需要对电堆有较为清楚的了解,需要大量的先验知识。2)而且建模过程繁琐,需要确定大量的模型参数(模型参数往往因为电堆型号的不同而产生较大的出入)。3)所建模型往往因为结构问题导致无法直接用于控制仿真。
为了克服以上问题,许多学者对如何减小建模复杂度进行了相关研究。徐夏吟.PEMFC温度子空间辨识模型[J].工业控制计算机,2015(9):77-78.仅利用燃料电池的温度数据辨识得到了单输入单输出状态空间模型,但是模型误差较大,且该方法对非线性特性的描述能力欠佳。Buchholz M,Eswein M,Krebs V.Modelling PEM fuel cell stacks forFDI using linear subspace identification[J].2008:341-346.该文利用子空间辨识的方法得到电压模型,但是对于电压这一典型的非线性特性环节的描述精度依然有待提高。
发明内容
本发明公开了一种质子交换膜燃料电池非线性状态空间模型辨识方法。
实现本发明目的的技术解决方案为:一种质子交换膜燃料电池非线性状态空间模型辨识方法,包括以下步骤:
步骤一、采集实验数据
根据控制要求选择合适的输入输出变量,并在测控平台上在无控制作用下采集实验
数据,实验数据以1000-2000组为宜;
步骤二、数据预处理
对采集的实测数据进行数据完整性分析、数据滤波,剔除一部分变化率超过实际
PEMFC系统工作实际的测试数据;
步骤三、构造输入输出Hankel矩阵
构造方式参考离散状态空间方程迭代得到的广义输入输出方程,实测组数据输入与输出矩阵的构造方式相同,满足关系式2i+j-1=N,i为Hankel矩阵行数,j为Hankel矩阵列数,N为数据组数。
u为输入变量,R为实数;
步骤四、斜投影的计算
求解输出数据行空间在输入数据行空间上的正交投影
其中参数矩阵A,B,从以下对偶系统求出
其中kp,S从下述矩阵求得
其中p,q=1,...,j
步骤五、求取广义能观矩阵和系统状态矩阵
采取传统线性子空间辨识的奇异值分解(SVD)
式中W1∈Ri×i,W2∈Ri×i为单位加权矩阵,W1满秩且W2满足rank(Wp)=rank(Wp*W2),由矩阵理论得到模型的广义观测矩阵Γi,Γi-1
Γi-1=Γi(1:l(i-1),:)
进一步确定模型状态序列为
符号“+”为Moore-Penrose伪逆;
步骤六、求取模型矩阵A,B,C,D
模型矩阵A,C可由下式得到
ΘAC利用下式的对偶系统求取
式中kBD,kBD可同样按照步骤四中kp,S得到,γ为可调参数;
模型矩阵B,D可由下式得到
ΘBD可由下式得到
ABD,BBD由本步骤第二步求解对偶系统得到。
步骤七、求取非线性特性
对步骤六最后一个公式右边进行奇异值分解计算得到非线性状态空间模型非线性特性部分的估计。
步骤八、利用MATLAB曲线拟合工具箱求取非线性多项式。
多项式阶次以四阶为宜。
本发明与现有技术相比,其显著优点:1)无需电堆模型先验知识,辨识速度快;2)通用性好;3)非线性特性拟合精度高;4)便于后续的控制器设计和仿真。
附图说明
图1是本发明一种质子交换膜燃料电池非线性状态空间模型辨识方法的流程图。
图2是50W空冷型电堆实测氢气流量信号。
图3是50W空冷型电堆实测负载电流信号。
图4是在采用相同输入信号的情况下,采用本发明辨识方法得到的PEMFC非线性状态空间模型输出电压与实测电堆输出电压的对比。
具体实施方式
本发明一种质子交换膜燃料电池非线性状态空间模型辨识方法,无需复杂的电堆先验知识和参数,提高了模型的描述精度。首先选择氢气流量和负载电流作为输入变量,选择电压作为输出变量,采集大量的数据。其次,通过采集的数据构造Hankel矩阵,求解对偶矩阵。然后,利用矩阵方程的结果构造矩阵投影。最后,对斜投影利用奇异值分解得到系统的状态序列估计。重复使用上述方法求解系统矩阵,并利用低秩近似技术获得系统非线性特性估计。以实测输入数据为横坐标,辨识非线性特性作为纵坐标利用MATLAB曲线拟合工具箱对系统非线性特性进行多项式拟合。该方法可以对质子交换膜燃料电池的非线性特性进行十分准确的描述,不仅对实际工程人员对质子交换膜燃料电池建模以及后续控制系统设计提供解决方案,而且对此类系统建模均有较好的参考价值。
下面结合附图1、2、3、4对本发明作进一步说明。
本发明一种质子交换膜燃料电池非线性状态空间模型辨识方法过程如下:
步骤一、采集实测数据
输入输出变量的选择。作为新型发电装置,保证质子交换膜燃料电池的输出电压稳定自然是首要任务。然而对燃料电池的电特性有影响的因素有很多,例如氢气/氧气流量、氢气/氧气压力、电堆/环境温度、负载电流等。氢气流量直接影响电堆氧化还原反应的效率,且流量变化已于调节和测量,只需要控制氢瓶的流量阀即可。因此,将氢气流量作为一个输入。本文的建模对象为空冷型自增湿电堆,氧气流量通过风扇调节和改变,无法调节和测量。因此不将氧气流量作为输入变量。由于实测实验条件假设电堆温度通过风扇调节维持恒定,而环境温度变化微乎其微,因此不考虑两者作为输入变量。负载电流对电压的影响是最明显、最直接。因此将负载电流作为系统输入变量。
步骤二、数据预处理
对实测进行数据完整性分析,数据滤波的准备工作。剔除一部分变化率不符合PEMFC系统工作实际的测试数据。
步骤三、构造输入/输出Hankel矩阵。
利用质子交换膜燃料电池测试平台采集N=1200组输入输出数据根据下式,利用前600组数据构造输入Hankel矩阵:
i和j满足关系式2i+j-1=N。以同样的方式可以构造输出Hankel矩阵,其中Yp,Yf∈Ril×j,Yp +∈R(i+1)l×j,Yp -∈R(i-1)l×j并定义:
步骤四、计算斜投影。
从如下对偶系统求解
式中,A∈Rj×li,B∈R2im×li分别是包含一组Lagrange乘子的矩阵,T,S定义如下:
矩阵kp∈Rj×j,kf∈Rj×j,其元素分别为:
p,q=1,...,j
斜投影Oi计算如下:
斜投影Oi+1计算类似,如下:
式中可通过对偶系统求解。
步骤五、求取系统状态估计和广义观测矩阵估计。
对斜投影Oi进行SVD分解
式中W1∈Ri×i,W2∈Ri×i为单位加权矩阵,W1满秩且W2满足rank(Wp)=rank(Wp*W2)。由矩阵理论可以得到系统的广义观测矩阵Γi,Γi-1
Γi-1=Γi(1:l(i-1),:)
进一步可确定系统状态序列为
符号“+”为Moore-Penrose伪逆。
步骤六、计算系统矩阵。
质子交换膜燃料电池非线性状态空间模型的参数矩阵A,B,C,D和广义非线性特性可通过以下最小二乘问题得到
其中
上式所示的系统矩阵A,C可由如下所示的对偶问题求出
式中,
系统矩阵A,C的估计可以通过分解ΘAC得到,由矩阵方程可以得到下式
步骤七、求取广义非线性函数。
对上式右边进行奇异值分解计算可以得到ΘBD和非线性特性的估计
因此,质子交换膜燃料电池非线性状态空间模型参数矩阵的计算公式为
A=ΘAC(1:2,1:2),C=ΘAC(3:3,1:2)
B=ΘBD(1:2,1:2),D=ΘBD(3:3,1:2)
至此,已经全部完成质子交换膜燃料电池非线性状态空间模型参数矩阵A,B,C,D和非线性特性的辨识。
步骤七、求取非线性多项式。
根据辨识得到的非线性函数通过MATLAB曲线拟合工具箱得到非线性环节的多项式。实际辨识多项式阶次为四阶。利用多项式和辨识系统矩阵可以搭建出系统模型。利用剩下的600组实测数据对模型进行验证。从附图4可以看出,辨识得到的模型可以对质子交换膜燃料电池的非线性进行很好的描述。
Claims (6)
1.一种质子交换膜燃料电池非线性状态空间模型辨识方法,其特征在于步骤如下:
步骤一、采集实验数据
根据控制要求选择合适的输入输出变量,并在测控平台上在无控制作用下采集实验数据;
步骤二、数据预处理
对采集的实测数据进行数据完整性分析、数据滤波,剔除一部分变化率超过实际PEMFC系统工作实际的测试数据;
步骤三、构造输入输出Hankel矩阵
构造方式参考离散状态空间方程迭代得到的广义输入输出方程,实测数据输入与输出矩阵的构造方式相同,满足关系式2i+j-1=N,i为Hankel矩阵行数,j为Hankel矩阵列数,N为数据组数;
<mrow>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mn>0</mn>
<mo>|</mo>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mi>p</mi>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mi>f</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mn>0</mn>
<mo>|</mo>
<mi>i</mi>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>U</mi>
<mi>p</mi>
<mo>+</mo>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mi>U</mi>
<mi>f</mi>
<mo>-</mo>
</msubsup>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
步骤四、斜投影的计算
求解输出数据行空间在输入数据行空间上的正交投影
<mrow>
<msub>
<mi>O</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<msup>
<mi>A</mi>
<mi>T</mi>
</msup>
<msub>
<mi>k</mi>
<mi>p</mi>
</msub>
<mo>+</mo>
<msup>
<mi>B</mi>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>:</mo>
<mi>m</mi>
<mi>i</mi>
<mo>,</mo>
<mo>:</mo>
<mo>)</mo>
</mrow>
<msup>
<mi>S</mi>
<mi>T</mi>
</msup>
<mrow>
<mo>(</mo>
<mo>:</mo>
<mo>,</mo>
<mn>1</mn>
<mo>:</mo>
<mi>m</mi>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mover>
<mi>L</mi>
<mo>^</mo>
</mover>
<mi>y</mi>
</msub>
<msub>
<mi>Y</mi>
<mi>p</mi>
</msub>
</mrow>
其中参数矩阵A,B,从以下对偶系统求出
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>Y</mi>
<mi>P</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mi>Y</mi>
<mi>p</mi>
<mi>T</mi>
</msubsup>
</mtd>
<mtd>
<mrow>
<msub>
<mi>k</mi>
<mi>p</mi>
</msub>
<mo>+</mo>
<msub>
<mi>k</mi>
<mi>f</mi>
</msub>
<mo>+</mo>
<msup>
<mi>&gamma;</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
<mtd>
<mi>S</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msup>
<mi>S</mi>
<mi>T</mi>
</msup>
</mtd>
<mtd>
<mi>T</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>L</mi>
<mi>y</mi>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>A</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>B</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mi>Y</mi>
<mi>f</mi>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中kp,S从下述矩阵求得
<mrow>
<msub>
<mi>k</mi>
<mi>p</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>h</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>i</mi>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mi>k</mi>
<mi>k</mi>
</msup>
<mrow>
<mo>(</mo>
<msub>
<mi>u</mi>
<mrow>
<mi>h</mi>
<mo>+</mo>
<mi>p</mi>
<mo>-</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>,</mo>
<msub>
<mi>u</mi>
<mrow>
<mi>h</mi>
<mo>+</mo>
<mi>q</mi>
<mo>-</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
其中p,q=1,...,j
步骤五、求取广义能观矩阵和系统状态矩阵
采取传统线性子空间辨识的奇异值分解(SVD)
<mrow>
<msub>
<mi>W</mi>
<mn>1</mn>
</msub>
<msub>
<mi>O</mi>
<mi>i</mi>
</msub>
<msub>
<mi>W</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>U</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<msub>
<mi>S</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>S</mi>
<mn>2</mn>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "(" close = ")">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>V</mi>
<mn>1</mn>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mi>V</mi>
<mn>2</mn>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>&ap;</mo>
<msub>
<mi>U</mi>
<mn>1</mn>
</msub>
<msub>
<mi>S</mi>
<mn>1</mn>
</msub>
<msubsup>
<mi>V</mi>
<mn>1</mn>
<mi>T</mi>
</msubsup>
</mrow>
式中W1∈Ri×i,W2∈Ri×i为单位加权矩阵,W1满秩且W2满足rank(Wp)=rank(Wp*W2),由矩阵理论得到模型的广义观测矩阵Γi,Γi-1
<mrow>
<msub>
<mi>&Gamma;</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<msub>
<mi>L</mi>
<mn>1</mn>
</msub>
<msubsup>
<mi>S</mi>
<mn>1</mn>
<mrow>
<mn>1</mn>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msubsup>
</mrow>
Γi-1=Γi(1:l(i-1),:)
进一步确定模型状态序列为
<mrow>
<msub>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>i</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>&Gamma;</mi>
<mi>i</mi>
<mo>+</mo>
</msubsup>
<msub>
<mi>O</mi>
<mi>i</mi>
</msub>
</mrow>
<mrow>
<msub>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<msubsup>
<mi>&Gamma;</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>+</mo>
</msubsup>
<msub>
<mi>O</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
</mrow>
符号“+”为Moore-Penrose伪逆;
步骤六、求取模型矩阵A,B,C,D
模型矩阵A,C由下式得到
<mrow>
<msub>
<mi>&Theta;</mi>
<mrow>
<mi>A</mi>
<mi>C</mi>
</mrow>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>A</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>C</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
ΘAC利用下式的对偶系统求取
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mover>
<mi>X</mi>
<mo>~</mo>
</mover>
<mi>i</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mover>
<mi>X</mi>
<mo>~</mo>
</mover>
<mi>i</mi>
<mi>T</mi>
</msubsup>
</mtd>
<mtd>
<mrow>
<msub>
<mi>k</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mo>+</mo>
<msubsup>
<mi>&gamma;</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mi>I</mi>
</mrow>
</mtd>
<mtd>
<msub>
<mi>S</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msubsup>
<mi>S</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
<mi>T</mi>
</msubsup>
</mtd>
<mtd>
<msub>
<mi>T</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>&Theta;</mi>
<mrow>
<mi>A</mi>
<mi>C</mi>
</mrow>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>B</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mi>&chi;</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
式中kBD,kBD同样按照步骤四中kp,S得到;
模型矩阵B,D由下式得到
<mrow>
<msub>
<mi>&Theta;</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>B</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>D</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
ΘBD由下式得到
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&Theta;</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mo>:</mo>
<mo>,</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<msub>
<mover>
<mi>f</mi>
<mo>&OverBar;</mo>
</mover>
<mi>k</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>u</mi>
<mn>0</mn>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>,</mo>
<msub>
<mover>
<mi>f</mi>
<mo>&OverBar;</mo>
</mover>
<mi>k</mi>
</msub>
<mo>(</mo>
<mrow>
<msub>
<mi>u</mi>
<mn>1</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
<mo>)</mo>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<msub>
<mover>
<mi>f</mi>
<mo>&OverBar;</mo>
</mover>
<mi>k</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>u</mi>
<mrow>
<mi>N</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mo>&rsqb;</mo>
<mo>=</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>A</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
<mi>T</mi>
</msubsup>
<msup>
<mi>K</mi>
<mi>k</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>:</mo>
<mi>i</mi>
<mo>+</mo>
<mi>j</mi>
<mo>,</mo>
<mo>:</mo>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msubsup>
<mi>B</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
<mi>T</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>,</mo>
<mo>:</mo>
<mo>)</mo>
</mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>t</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<msubsup>
<mi>K</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
<mi>k</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>,</mo>
<mo>:</mo>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
ABD,BBD由本步骤中的求解对偶系统得到;
步骤七、求取非线性特性
对步骤六最后一个公式右边进行奇异值分解,计算得到非线性状态空间模型非线性特性部分的估计;
步骤八、利用MATLAB曲线拟合工具箱求取非线性多项式。
2.根据权利要求1所述的辨识方法,其特征在于:步骤一所述采集实验数据组数优选1000-2000组。
3.根据权利要求1所述的辨识方法,其特征在于:步骤二所述变化率为
4.根据权利要求1所述的辨识方法,其特征在于:步骤三所述Up,Uf由下式求出,
<mrow>
<msub>
<mi>U</mi>
<mrow>
<mn>0</mn>
<mo>|</mo>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mn>0</mn>
<mo>|</mo>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mi>p</mi>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mi>f</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mn>0</mn>
<mo>|</mo>
<mi>i</mi>
</mrow>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>U</mi>
<mrow>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>|</mo>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mover>
<mo>=</mo>
<mi>&Delta;</mi>
</mover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msubsup>
<mi>U</mi>
<mi>p</mi>
<mo>+</mo>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msubsup>
<mi>U</mi>
<mi>f</mi>
<mo>-</mo>
</msubsup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>.</mo>
</mrow>
5.根据权利要求1所述的辨识方法,其特征在于:步骤四所述T,Sq,kf分别由以下三式求出
<mrow>
<msub>
<mi>S</mi>
<mi>q</mi>
</msub>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>t</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msup>
<mi>K</mi>
<mn>1</mn>
</msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msup>
<mi>K</mi>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<mrow>
<msup>
<mi>K</mi>
<mi>m</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
<mrow>
<msub>
<mi>k</mi>
<mi>f</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>h</mi>
<mo>=</mo>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mrow>
<mn>2</mn>
<mi>i</mi>
</mrow>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mi>k</mi>
<mi>k</mi>
</msup>
<mrow>
<mo>(</mo>
<msub>
<mi>u</mi>
<mrow>
<mi>h</mi>
<mo>+</mo>
<mi>p</mi>
<mo>-</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>,</mo>
<msub>
<mi>u</mi>
<mrow>
<mi>h</mi>
<mo>+</mo>
<mi>q</mi>
<mo>-</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
6.根据权利要求1所述的辨识方法,其特征在于:步骤六所述
<mrow>
<msub>
<mi>k</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msubsup>
<mi>K</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
<mi>K</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>p</mi>
<mo>,</mo>
<mi>q</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
<msub>
<mi>A</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mo>&Element;</mo>
<msup>
<mi>R</mi>
<mrow>
<mi>j</mi>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>,</mo>
<msub>
<mi>B</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mo>&Element;</mo>
<msup>
<mi>R</mi>
<mrow>
<mi>m</mi>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>,</mo>
</mrow>
<mrow>
<msub>
<mi>T</mi>
<mrow>
<mi>B</mi>
<mi>D</mi>
</mrow>
</msub>
<mo>=</mo>
<mi>d</mi>
<mi>i</mi>
<mi>a</mi>
<mi>g</mi>
<mo>{</mo>
<msubsup>
<mn>1</mn>
<mi>N</mi>
<mi>T</mi>
</msubsup>
<msup>
<mi>K</mi>
<mi>k</mi>
</msup>
<msub>
<mn>1</mn>
<mi>N</mi>
</msub>
<mo>}</mo>
<mo>,</mo>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<mi>m</mi>
</mrow>
SBD=[Si+1 … Si+j]T.。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710772967.5A CN107632522B (zh) | 2017-08-31 | 2017-08-31 | 一种质子交换膜燃料电池非线性状态空间模型辨识方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710772967.5A CN107632522B (zh) | 2017-08-31 | 2017-08-31 | 一种质子交换膜燃料电池非线性状态空间模型辨识方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107632522A true CN107632522A (zh) | 2018-01-26 |
CN107632522B CN107632522B (zh) | 2020-06-19 |
Family
ID=61101003
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710772967.5A Active CN107632522B (zh) | 2017-08-31 | 2017-08-31 | 一种质子交换膜燃料电池非线性状态空间模型辨识方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107632522B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109143074A (zh) * | 2018-06-28 | 2019-01-04 | 中国科学院光电研究院 | 一种动力电池模型参数辨识方法及系统 |
CN109921072A (zh) * | 2019-03-20 | 2019-06-21 | 南京理工大学 | 一种质子交换膜燃料电池输出功率的预测控制方法 |
CN113110068A (zh) * | 2021-05-22 | 2021-07-13 | 北京理工大学 | 一种子空间系统辨识方法及系统 |
CN114137829A (zh) * | 2021-09-28 | 2022-03-04 | 南京理工大学 | 基于almbo优化算法的质子交换膜燃料电池子空间辨识方法 |
Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102103376A (zh) * | 2010-12-16 | 2011-06-22 | 浙江大学 | 基于子空间辨识的pid回路控制性能评估方法 |
CN102968056A (zh) * | 2012-12-07 | 2013-03-13 | 上海电机学院 | 质子交换膜燃料电池的建模系统及其智能预测控制方法 |
CN103970964A (zh) * | 2014-05-23 | 2014-08-06 | 哈尔滨工业大学 | 一种挠性卫星模态参数在轨辨识方法 |
CN104993480A (zh) * | 2015-07-22 | 2015-10-21 | 福州大学 | 基于递推随机子空间的电力系统低频振荡在线辨识方法 |
CN105425587A (zh) * | 2015-11-16 | 2016-03-23 | 北京理工大学 | 迟滞非线性电机辨识与控制方法 |
CN105446347A (zh) * | 2015-11-30 | 2016-03-30 | 上海卫星工程研究所 | 针对卫星太阳电池阵的在轨模态辨识系统及方法 |
DE112014006345T5 (de) * | 2014-02-07 | 2016-10-20 | Mitsubishi Electric Corporation | Systemidentifikationsvorrichtung |
CN106505587A (zh) * | 2016-10-19 | 2017-03-15 | 福州大学 | 基于广义形态滤波与改进mp算法的低频振荡模态辨识方法 |
CN106546847A (zh) * | 2016-10-20 | 2017-03-29 | 西南交通大学 | 基于prce的低频振荡模式在线辨识方法 |
CN106654319A (zh) * | 2016-12-27 | 2017-05-10 | 东南大学 | 一种基于变异粒子群和差分进化混合算法的pemfc系统温度建模方法 |
-
2017
- 2017-08-31 CN CN201710772967.5A patent/CN107632522B/zh active Active
Patent Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102103376A (zh) * | 2010-12-16 | 2011-06-22 | 浙江大学 | 基于子空间辨识的pid回路控制性能评估方法 |
CN102968056A (zh) * | 2012-12-07 | 2013-03-13 | 上海电机学院 | 质子交换膜燃料电池的建模系统及其智能预测控制方法 |
DE112014006345T5 (de) * | 2014-02-07 | 2016-10-20 | Mitsubishi Electric Corporation | Systemidentifikationsvorrichtung |
CN103970964A (zh) * | 2014-05-23 | 2014-08-06 | 哈尔滨工业大学 | 一种挠性卫星模态参数在轨辨识方法 |
CN104993480A (zh) * | 2015-07-22 | 2015-10-21 | 福州大学 | 基于递推随机子空间的电力系统低频振荡在线辨识方法 |
CN105425587A (zh) * | 2015-11-16 | 2016-03-23 | 北京理工大学 | 迟滞非线性电机辨识与控制方法 |
CN105446347A (zh) * | 2015-11-30 | 2016-03-30 | 上海卫星工程研究所 | 针对卫星太阳电池阵的在轨模态辨识系统及方法 |
CN106505587A (zh) * | 2016-10-19 | 2017-03-15 | 福州大学 | 基于广义形态滤波与改进mp算法的低频振荡模态辨识方法 |
CN106546847A (zh) * | 2016-10-20 | 2017-03-29 | 西南交通大学 | 基于prce的低频振荡模式在线辨识方法 |
CN106654319A (zh) * | 2016-12-27 | 2017-05-10 | 东南大学 | 一种基于变异粒子群和差分进化混合算法的pemfc系统温度建模方法 |
Non-Patent Citations (2)
Title |
---|
LENG BOYANG等: "Research on Electrical Characteristics of State Space Modeling and Parameter Identifying of PEMFC", 《PROCEEDINGS OF THE 34TH CHINESE CONTROL CONFERENCE》 * |
戚志东等: "PEMFC动态建模与模糊分数阶PIλDµ控制", 《控制与决策》 * |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109143074A (zh) * | 2018-06-28 | 2019-01-04 | 中国科学院光电研究院 | 一种动力电池模型参数辨识方法及系统 |
CN109921072A (zh) * | 2019-03-20 | 2019-06-21 | 南京理工大学 | 一种质子交换膜燃料电池输出功率的预测控制方法 |
CN109921072B (zh) * | 2019-03-20 | 2022-04-01 | 南京理工大学 | 一种质子交换膜燃料电池输出功率的预测控制方法 |
CN113110068A (zh) * | 2021-05-22 | 2021-07-13 | 北京理工大学 | 一种子空间系统辨识方法及系统 |
CN113110068B (zh) * | 2021-05-22 | 2023-03-10 | 北京理工大学 | 一种子空间系统辨识方法及系统 |
CN114137829A (zh) * | 2021-09-28 | 2022-03-04 | 南京理工大学 | 基于almbo优化算法的质子交换膜燃料电池子空间辨识方法 |
CN114137829B (zh) * | 2021-09-28 | 2022-12-27 | 南京理工大学 | 基于almbo优化算法的质子交换膜燃料电池子空间辨识方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107632522B (zh) | 2020-06-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Yuan et al. | Developed coyote optimization algorithm and its application to optimal parameters estimation of PEMFC model | |
CN107632522A (zh) | 一种质子交换膜燃料电池非线性状态空间模型辨识方法 | |
Chen et al. | Mechanism analysis of starvation in PEMFC based on external characteristics | |
Lee et al. | Empirical modeling of polymer electrolyte membrane fuel cell performance using artificial neural networks | |
CN113097542B (zh) | 一种基于Amesim的燃料电池空气系统建模仿真方法 | |
CN108038340A (zh) | 一种质子交换膜燃料电池分数阶状态空间模型辨识方法 | |
CN112038671B (zh) | 一种固体氧化物燃料电池温度分布估计方法与系统 | |
Chavan et al. | System identification black box approach for modeling performance of PEM fuel cell | |
CN102968056A (zh) | 质子交换膜燃料电池的建模系统及其智能预测控制方法 | |
CN103336867B (zh) | 质子交换膜燃料电池模型优化处理方法 | |
CN106295082A (zh) | 一种平板式固体氧化物燃料电池的数值模拟方法 | |
Ural et al. | Mathematical models of PEM fuel cells | |
CN116231000A (zh) | 一种燃料电池/电堆仿真模型的构建方法 | |
CN114551937B (zh) | 一种燃料电池的性能检测系统及方法 | |
CN106557627A (zh) | 基于状态空间维纳模型的递推参数估计方法 | |
CN114357806B (zh) | 基于物质流接口的燃料电池电堆的双模式仿真方法和设备 | |
CN102693215A (zh) | 幂函数型水位流量关系的一种拟合方法 | |
Deng et al. | Data-driven reconstruction of interpretable model for air supply system of proton exchange membrane fuel cell | |
Tang et al. | Electrical power prediction of proton exchange membrane fuel cell by using support vector regression | |
Donateo | Semi-Empirical Models for Stack and Balance of Plant in Closed-Cathode Fuel Cell Systems for Aviation | |
Koeppel et al. | Use of a Reduced Order Model (ROM) to simulate SOFC performance in system models | |
Zhong et al. | A hybrid multi-variable experimental model for a PEMFC | |
Guo et al. | Marginalized particle filtering for online parameter estimation of PEMFC applied to hydrogen UAVs | |
Bao et al. | A computationally efficient steady-state electrode-level and 1D+ 1D cell-level fuel cell model | |
Ansari | Reduced-order modeling of PEM fuel cell based on POD and PODI: an efficient approach toward combining highest accuracy with real-time performance |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |