基于高阶全驱的质子交换膜燃料电池阴极进气系统控制方法

Cathode intake system control method of proton exchange membrane fuel cells based on high-order fully actuated system

摘要

针对质子交换膜燃料电池进气系统中存在的多输入多输出、强非线性强耦合、内外部扰动等问题，提出了一种新的控制策略实现对进气流量和阴极压力精确协同控制。本文将传统状态空间的四阶燃料电池欠驱系统转化为高阶全驱系统，同时基于高阶全驱系统设计扩张状态观测器来估计未知扰动。根据高阶全驱系统的特性推导出控制律，实现对进气流量和阴极压力的协同控制。最后，通过与其他控制策略的仿真对比分析，并且在实验台架上进行验证，结果表明可以实现对进气流量和阴极压力的精确协同控制。

Abstract

A new control strategy was proposed for the air supply system of proton exchange membrane fuel cell， which solves the multiple input and multiple output， strong nonlinear coupling and internal and external disturbance control problem of the air supply system， then realizes the precision accurate tracking control of the air flow rate and cathode pressure. The high-order fully actuated system is derived from fourth-order underactuated system of the traditional state space model， and the extended state observer is designed to estimate the unknown disturbance based on high-order fully actuated system. According to the characteristics of high-order fully actuated system， the control law is designed to realize the cooperative control of the air flow and cathode pressure. Finally， through the simulation result comparison with other control strategies， and verified on the experimental bench， the results show that the air flow rate and cathode pressure can be precise controlled.

Graphical abstract

关键词

自动控制技术 / 质子交换膜燃料电池 / 阴极进气系统 / 高阶全驱方法 / 扩张状态观测器

Key words

automatic control technology / proton exchange membrane fuel cell / air supply control / high-order fully actuated system approaches / extended state observer

0 引　言

近年来，随着人们对碳排放和全球气候变化的担忧加剧，几乎所有国家都在减少碳排放。为了应对这些挑战，越来越多的全球汽车制造商将目光转向了质子交换膜燃料电池（Proton exchange mem-brane fuel cell， PEMFC）汽车^［1］。质子交换膜燃料电池是一种将氢和氧的化学能转化为电能的电化学装置，已广泛应用于各个领域^［2］。在汽车应用中，PEMFC因其低噪声、零污染、高效率和高功率密度而被认为是传统发动机最有前途的替代品^［3，4］。

对于用于汽车的大功率PEMFC，供气系统的瞬态特征对实现PEMFC的期望功率和延长其使用寿命起着至关重要的作用^［5，6］。在众多的特性中，进入阴极的空气流量和阴极压力是两个关键变量^［7，8］。空气流速过低可能导致缺氧，从而使电化学反应中氧气不足而损坏质子交换膜。相反，空气流速过高使寄生功率增加到PEMFC总功率的近30%^［9］。此外，空气流速的变化会影响阴极压力，过大的阴极压力，可以增加通过扩散层到达催化层的空气量，从而提高催化速率。然而，阴极和阳极之间的瞬时压差过大可能会影响质子交换膜，导致PEMFC产生不可逆损坏。

为解决这些问题，目前国内外研究人员提出了许多控制策略。Pukrushpan等^［10］建立了面向控制的PEMFC进气系统的九阶模型，这是一个至今仍被广泛使用的经典模型，他们提出了线性最优控制来调节氧气流量，从而避免缺氧；Talj等^［11］提出了一种用于PEMFC的鲁棒策略迭代控制方案，来解决进气系统的数学模型中存在不确定性等问题；Ma等^［12］提出了一种非线性三步法，控制方案分为稳态控制、前馈控制和误差反馈控制3部分；Liu等^［13］提出了一种高阶滑模控制，来解决控制对干扰和不确定性不敏感的问题，以避免缺氧现象；Wang等^［14］提出了基于区间的Type-2模糊逻辑系统的自适应控制， 以描述未建模的系统动力学。同时反馈线性化和无模型自适应控制等技术也在PEMFC进气系统控制中取得了显著效果。

上述控制方案均基于Pukrushpan提出的PEMFC进气系统经典状态空间模型进行设计。为降低控制器设计的复杂性，对该模型进行了多种简化。然而，迄今为止，最简单的PEMFC送风系统模型仍然是三阶的。PEMFC的辅助设备如空压机和供气歧管仍具有较强的非线性，且供气系统的控制输入仅为一维。而且非线性动力学需要用Lyapunov稳定性设计方法来处理。在某些复杂情况下，甚至不能实现系统的局部稳定。近年来，Duan^［15］针对非线性系统提出了一种新颖的高阶全驱（High-order fully actuated，HOFA）系统控制方法，他指出HOFA系统是对动态系统的一种描述。与状态空间模型类似，HOFA系统是一个面向控制的模型。全驱动模型能够消除开环系统的所有动态特性，包括线性或非线性特征。这样就很容易建立新的闭环动态特性。与传统的非线性方法相比，HOFA系统方法提供了更大的设计自由度。

本文提出基于HOFA系统理论的PEMFC进气系统控制方案，给出了将燃料电池欠驱系统转换为HOFA系统模型的具体步骤。使用扩张状态观测器来估计模型中的未知扰动和未建模动态。通过一系列仿真和实验验证了所提出的基于HOFA方法的有效性。

1 系统建模与分析

本节主要叙述PEMFC阴极进气系统的建模过程，分别建立了空气压缩机模型、供气歧管模型、阴极模型和节气门模型。

1.1　控制问题

质子交换膜燃料电池系统如图1所示，主要包括送风系统、供氢系统、温度管理系统和湿度系统。空气进入阴极之前在压缩机中进行压缩，压缩空气依次通过中间冷却器和加湿器，然后进入燃料电池堆，与质子交换膜上的氢离子发生电化学反应。剩余的空气通过水-空气分离装置，其中水蒸气被排放到水箱中，而干燥的空气则通过节流阀排出。在此过程中，压缩机转速和节气门角度作为控制变量，共同调节进入阴极的气流和阴极处的压力。

1.2　PEMFC阴极动态建模

本文通过控制空压机转速来调节进入阴极的气流。采用电机驱动转矩方程来表征空压机动态特性：

J_{\mathrm{c}\mathrm{p}}\frac{\mathrm{d}{\omega }_{\mathrm{c}\mathrm{p}}}{\mathrm{d}t}={\tau }_{\mathrm{c}\mathrm{m}}-{\tau }_{\mathrm{c}\mathrm{p}}（1）

式中：J_{\mathrm{c}\mathrm{p}}为电机的转动惯量；{\omega }_{\mathrm{c}\mathrm{p}}为压缩机的角速度，与压缩机转速N_{\mathrm{c}\mathrm{p}}的关系为{\omega }_{\mathrm{c}\mathrm{p}}=\mathrm{2}\mathrm{\pi }{N}_{\mathrm{c}\mathrm{p}}/\mathrm{60}；{\tau }_{cm}为电机的驱动转矩；{\tau }_{\mathrm{c}\mathrm{p}}为电机的负载转矩。

{\tau }_{\mathrm{c}\mathrm{m}}=\frac{k_t{\eta }_{\mathrm{c}\mathrm{m}}}{R_{\mathrm{c}\mathrm{m}}}\left(V_{\mathrm{c}\mathrm{m}}-k_v{\omega }_{\mathrm{c}\mathrm{p}}\right)（2）

{\tau }_{\mathrm{c}\mathrm{p}}=\frac{C_p{T}_{\mathrm{a}\mathrm{t}\mathrm{m}}}{{\omega }_{\mathrm{c}\mathrm{p}}{\eta }_{\mathrm{c}\mathrm{p}}}\left[{\left(\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}\right)}^{\frac{\gamma -\mathrm{1}}{\gamma }}-\mathrm{1}\right]W_{\mathrm{c}\mathrm{p}}（3）

式中：{\eta }_{\mathrm{c}\mathrm{m}}、k_t、k_v和R_{\mathrm{c}\mathrm{m}}为电机常数；V_{\mathrm{c}\mathrm{m}}为空压机电压；C_p为恒压下的空气比热容；\gamma为空气的比热比；P_{\mathrm{a}\mathrm{t}\mathrm{m}}为大气压力；{\eta }_{\mathrm{c}\mathrm{p}}为空压机最大工作效率；T_{\mathrm{a}\mathrm{t}\mathrm{m}}为标准大气压力下温度；W_{\mathrm{c}\mathrm{p}}为空压机出口处空气的质量流量，由二阶多项式拟合得到。

\begin{array}{l}{W}_{\mathrm{c}\mathrm{p}}\left(N_{\mathrm{c}\mathrm{p}},\frac{{\mathrm{P}}_{\mathrm{s}\mathrm{m}}}{{\mathrm{P}}_{\mathrm{a}\mathrm{t}\mathrm{m}}}\right)=p_{\mathrm{1}}{N}_{\mathrm{c}\mathrm{p}}+p_{\mathrm{2}}{N}_{\mathrm{c}\mathrm{p}}\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}+p_{\mathrm{3}}\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}+\\ p_{\mathrm{4}}{N}_{\mathrm{c}\mathrm{p}}^{\mathrm{2}}+p_{\mathrm{5}}{\left(\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}\right)}^{\mathrm{2}}+p_{\mathrm{6}}{N}_{\mathrm{c}\mathrm{p}}^{\mathrm{2}}\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}+p_{\mathrm{7}}{N}_{\mathrm{c}\mathrm{p}}{\left(\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}\right)}^{\mathrm{2}}+p_{\mathrm{8}}\end{array}（4）

式中：p_i为拟合系数，i\in \{\mathrm{1},\dots ,\mathrm{8}\}，p_i的具体参数可从文献［16］可获。

根据理想气体质量守恒定律，供气歧管压力为：

\frac{\mathrm{d}{P}_{\mathrm{s}\mathrm{m}}}{\mathrm{d}t}=\frac{R{T}_{\mathrm{c}\mathrm{p}}}{V_{\mathrm{s}\mathrm{m}}{M}_{\mathrm{a}}}(W_{\mathrm{c}\mathrm{p}}-W_{\mathrm{s}\mathrm{m}})（5）

式中：P_{\mathrm{s}\mathrm{m}}为进气歧管压力；R为一般气体常数；T_{\mathrm{c}\mathrm{p}}为空压机的出口温度；V_{\mathrm{s}\mathrm{m}}为进气歧管容积；M_{\mathrm{a}}为空气的摩尔质量；W_{\mathrm{s}\mathrm{m}}为进气歧管出口处空气质量流量。

空气压缩机的出口温度T_{\mathrm{c}\mathrm{p}}为：

T_{\mathrm{c}\mathrm{p}}=T_{\mathrm{a}\mathrm{t}\mathrm{m}}+\frac{T_{\mathrm{a}\mathrm{t}\mathrm{m}}}{{\eta }_{\mathrm{c}\mathrm{p}}}\left[{\left(\frac{P_{\mathrm{s}\mathrm{m}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}}\right)}^{\frac{\gamma -\mathrm{1}}{\gamma }}-\mathrm{1}\right]（6）

通过喷嘴流动方程的线性化形式计算进气歧管出口处的空气质量流量W_{\mathrm{s}\mathrm{m}}：

W_{\mathrm{s}\mathrm{m}}=K_{\mathrm{s}\mathrm{m}}(P_{\mathrm{s}\mathrm{m}}-P_{\mathrm{c}\mathrm{a}})（7）

式中：K_{\mathrm{s}\mathrm{m}}为供应歧管常数。

使用理想气体方程来描述阴极压力P_{\mathrm{c}\mathrm{a}}动态：

\frac{\mathrm{d}{P}_{\mathrm{c}\mathrm{a}}}{\mathrm{d}t}=\frac{T_{\mathrm{f}\mathrm{c}}}{V_{\mathrm{c}\mathrm{a}}}\left(\frac{R}{M_{\mathrm{a}}}{W}_{\mathrm{s}\mathrm{m}}-\frac{R}{M_{\mathrm{a}}}{W}_{\mathrm{c}\mathrm{a},\mathrm{o}\mathrm{u}\mathrm{t}}-\frac{R}{M_{{\mathrm{O}}_{\mathrm{2}}}}{W}_{{\mathrm{O}}_{\mathrm{2}},\mathrm{r}\mathrm{e}\mathrm{c}}\right)（8）

式中：T_{\mathrm{f}\mathrm{c}}为燃料电池的工作温度；V_{\mathrm{c}\mathrm{a}}为阴极的体积；M_{\mathrm{a}}和M_{{\mathrm{O}}_{\mathrm{2}}}分别为空气和氧气的摩尔质量；W_{\mathrm{c}\mathrm{a},\mathrm{o}\mathrm{u}\mathrm{t}}为阴极出口处空气的质量流量；W_{{\mathrm{O}}_{\mathrm{2}},\mathrm{r}\mathrm{e}\mathrm{c}}为阴极反应中消耗氧气的质量流量。

通过法拉第定律计算电化学反应中消耗氧气的质量流量：

W_{{\mathrm{O}}_{\mathrm{2}},\mathrm{r}\mathrm{e}\mathrm{c}}=M_{{\mathrm{O}}_{\mathrm{2}}}\frac{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}\cdot I_{\mathrm{s}\mathrm{t}}}{\mathrm{4}F}（9）

式中：n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}为PEMFC中串联电池单元的个数；F为法拉第常数；I_{\mathrm{s}\mathrm{t}}为负载电流。

用非线性喷嘴方程计算阴极出口的空气质量流量：

W_{\mathrm{c}\mathrm{a},\mathrm{o}\mathrm{u}\mathrm{t}}=

\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\left(\frac{\theta \mathrm{\pi }}{\mathrm{180}}\right)\frac{C_{\mathrm{D},\mathrm{t}\mathrm{r}}{A}_{\mathrm{T},\mathrm{t}\mathrm{r}}{P}_{\mathrm{c}\mathrm{a}}}{\sqrt[]{R{T}_{\mathrm{f}\mathrm{c}}}}{(\gamma )}^{\frac{\mathrm{1}}{\mathrm{2}}}{\left(\frac{\mathrm{2}}{\gamma +\mathrm{1}}\right)}^{\frac{\gamma +\mathrm{1}}{\mathrm{2}(\gamma -\mathrm{1})}}（10）

式中：\theta为节气门的开合角度；C_{\mathrm{D},\mathrm{t}\mathrm{r}}节气门流量系数；A_{\mathrm{T},\mathrm{t}\mathrm{r}}为节气门的最大开启面积。

在燃料电池模型中，假定节气门受闭环控制，并且节气门的动态性能可以通过一阶系统来近似。

\frac{\mathrm{d}\theta }{\mathrm{d}t}=\frac{\mathrm{1}}{T_{\mathrm{t}\mathrm{r}}}\left({\theta }^{\mathrm{*}}-\theta \right)（11）

式中：T_{\mathrm{t}\mathrm{r}}为节气门响应时间常数；{\theta }^{\mathrm{*}}为节气门开启角度命令。

燃料电池系统模型参数从文献［16］可得。由式（1）~（11）可以得到状态空间方程为：

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}}{\mathrm{d}t}=-b_{\mathrm{1}}{x}_{\mathrm{1}}-\frac{b_{\mathrm{2}}}{x_{\mathrm{1}}}\left({\left(\frac{x_{\mathrm{2}}}{b_{\mathrm{3}}}\right)}^{b_{\mathrm{4}}}-\mathrm{1}\right)W_{\mathrm{c}\mathrm{p}}+b_{\mathrm{5}}{u}_{\mathrm{1}}\\ \frac{\mathrm{d}{x}_{\mathrm{2}}}{\mathrm{d}t}=b_{\mathrm{6}}\left[\mathrm{1}+b_{\mathrm{7}}\left({\left(\frac{x_{\mathrm{2}}}{b_{\mathrm{3}}}\right)}^{b_{\mathrm{4}}}-\mathrm{1}\right)\right]\left[W_{\mathrm{c}\mathrm{p}}-b_{\mathrm{8}}\left(x_{\mathrm{2}}-x_{\mathrm{3}}\right)\right]\\ \frac{\mathrm{d}{x}_{\mathrm{3}}}{\mathrm{d}t}=b_{\mathrm{9}}{x}_{\mathrm{2}}-b_{\mathrm{10}}\xi -\left(b_{\mathrm{11}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\left(\frac{x_{\mathrm{4}}\mathrm{\pi }}{\mathrm{180}}\right)+b_{\mathrm{9}}\right)x_{\mathrm{3}}\\ \frac{\mathrm{d}{x}_{\mathrm{4}}}{\mathrm{d}t}=\frac{\mathrm{1}}{b_{\mathrm{12}}}\left(-x_{\mathrm{4}}+u_{\mathrm{2}}\right)\end{array}\right.（12）

式中：[u_{\mathrm{1}} u_{\mathrm{2}}{]}^{\mathrm{T}}=[V_{\mathrm{c}\mathrm{m}} {\theta }^{\mathrm{*}}{]}^{\mathrm{T}}为控制量；[x_{\mathrm{1}} x_{\mathrm{2}} x_{\mathrm{3}} x_{\mathrm{4}}{]}^{\mathrm{T}}=[{\omega }_{\mathrm{c}\mathrm{p}} P_{\mathrm{s}\mathrm{m}} P_{\mathrm{c}\mathrm{a}}{ \theta ]}^{\mathrm{T}}为系统的状态；\xi =I_{\mathrm{s}\mathrm{t}}为可测量的扰动；参数b_i，i\in \{\mathrm{1},\dots ,\mathrm{12}\}如表1所示。

系统的输出方程可以表示为：

y_{\mathrm{1}}=b_{\mathrm{8}}\left(x_{\mathrm{2}}-x_{\mathrm{3}}\right)（13）

y_{\mathrm{2}}=x_{\mathrm{3}}（14）

2 高阶全驱控制器设计

2.1　高阶全驱系统

为简化计算过程，将式（12）改写为：

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}}{\mathrm{d}t}=f_{\mathrm{1}}\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)+b_{\mathrm{5}}{u}_{\mathrm{1}}\\ \frac{\mathrm{d}{x}_{\mathrm{2}}}{\mathrm{d}t}=f_{\mathrm{2}}\left(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}\right)\\ \frac{\mathrm{d}{x}_{\mathrm{3}}}{\mathrm{d}t}=f_{\mathrm{3}}\left(x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}}\right)\\ \frac{\mathrm{d}{x}_{\mathrm{4}}}{\mathrm{d}t}=f_{\mathrm{4}}\left(x_{\mathrm{4}}\right)+\frac{\mathrm{1}}{b_{\mathrm{12}}}{u}_{\mathrm{2}}\end{array}\right.（15）

分别对式（15）中的{\dot{x}}_{\mathrm{2}}和{\dot{x}}_{\mathrm{3}}求导，可得：

\frac{{\mathrm{d}}^{\mathrm{2}}{x}_{\mathrm{2}}}{\mathrm{d}{t}^{\mathrm{2}}}=\frac{\partial f_{\mathrm{2}}(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})}{\partial x_{\mathrm{1}}}\frac{\mathrm{d}{x}_{\mathrm{1}}}{\mathrm{d}t}+\frac{\partial f_{\mathrm{2}}(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})}{\partial x_{\mathrm{2}}}\frac{\mathrm{d}{x}_{\mathrm{2}}}{\mathrm{d}t}+

\frac{\partial f_{\mathrm{2}}(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})}{\partial x_{\mathrm{3}}}\frac{\mathrm{d}{x}_{\mathrm{3}}}{\mathrm{d}t}（16）

\frac{{\mathrm{d}}^{\mathrm{2}}{x}_{\mathrm{3}}}{\mathrm{d}{t}^{\mathrm{2}}}=\frac{\partial f_{\mathrm{3}}(x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}})}{\partial x_{\mathrm{2}}}\frac{\mathrm{d}{x}_{\mathrm{2}}}{\mathrm{d}t}+\frac{\partial f_{\mathrm{3}}(x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}})}{\partial x_{\mathrm{3}}}\frac{\mathrm{d}{x}_{\mathrm{3}}}{\mathrm{d}t}+\frac{\partial f_{\mathrm{3}}(x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}})}{\partial x_{\mathrm{4}}}\frac{\mathrm{d}{x}_{\mathrm{4}}}{\mathrm{d}t}（17）

结合公式：

\left\{\begin{array}{l}\frac{\mathrm{d}{x}_{\mathrm{1}}}{\mathrm{d}t}=f_{\mathrm{1}}\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)+b_{\mathrm{5}}{u}_{\mathrm{1}}\\ \frac{\mathrm{d}{x}_{\mathrm{4}}}{\mathrm{d}t}=f_{\mathrm{4}}\left(x_{\mathrm{4}}\right)+\frac{\mathrm{1}}{b_{\mathrm{12}}}{u}_{\mathrm{2}}\end{array}\right.（18）

将式（18）代入式（16）和（17）中，简化可得：

\frac{{\mathrm{d}}^{\mathrm{2}}{x}_{\mathrm{2}}}{\mathrm{d}{t}^{\mathrm{2}}}=Q_{\mathrm{2}}(x,\dot{x})+

[L_{\mathrm{21}}(x,\dot{x}) L_{\mathrm{22}}(x,\dot{x})][u_{\mathrm{1}} u_{\mathrm{2}}{]}^{\mathrm{T}}（19）

\frac{{\mathrm{d}}^{\mathrm{2}}{x}_{\mathrm{3}}}{\mathrm{d}{t}^{\mathrm{2}}}=Q_{\mathrm{2}}(x,\dot{x})+

[L_{\mathrm{21}}(x,\dot{x}) L_{\mathrm{22}}(x,\dot{x})][u_{\mathrm{1}} u_{\mathrm{2}}{]}^{\mathrm{T}}（20）

式中：Q_i(x,\dot{x},\xi ,\dot{\xi })和L_{ij}(x)，i\in \{\mathrm{2,3}\}j\in \{\mathrm{1,2}\}为：

\begin{array}{l}{Q}_{\mathrm{2}}=b_{\mathrm{6}}\left[b_{\mathrm{4}}{b}_{\mathrm{7}}{\left(\frac{x_{\mathrm{2}}}{b_{\mathrm{3}}}\right)}^{b_{\mathrm{4}}-\mathrm{1}}\frac{\mathrm{1}}{b_{\mathrm{3}}}\right]{\dot{x}}_{\mathrm{2}}\left[h(x_{\mathrm{1}},x_{\mathrm{2}})-b_{\mathrm{8}}(x_{\mathrm{2}}-x_{\mathrm{3}})\right]+b_{\mathrm{6}}\left[\mathrm{1}+b_{\mathrm{7}}\left({\left(\frac{x_{\mathrm{2}}}{b_{\mathrm{3}}}\right)}^{b_{\mathrm{4}}}-\mathrm{1}\right)\right]\times \\ \left[\frac{\partial h(x_{\mathrm{1}},x_{\mathrm{2}})}{\partial x_{\mathrm{1}}}\left(-b_{\mathrm{1}}{x}_{\mathrm{1}}-\frac{b_{\mathrm{2}}}{x_{\mathrm{1}}}\left({\left(\frac{x_{\mathrm{2}}}{b_{\mathrm{3}}}\right)}^{b_{\mathrm{4}}}-\mathrm{1}\right)h(x_{\mathrm{1}},x_{\mathrm{2}})\right)+\right.\\ \left.\frac{\partial h(x_{\mathrm{1}},x_{\mathrm{2}})}{\partial x_{\mathrm{2}}}{\dot{x}}_{\mathrm{2}}-b_{\mathrm{8}}{\dot{x}}_{\mathrm{2}}+b_{\mathrm{8}}{\dot{x}}_{\mathrm{3}}\right];\end{array}

Q_{\mathrm{3}}=b_{\mathrm{9}}{\dot{x}}_{\mathrm{2}}-b_{\mathrm{10}}\dot{\xi }-b_{\mathrm{11}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\left(\frac{x_{\mathrm{4}}\mathrm{\pi }}{\mathrm{180}}\right){\dot{x}}_{\mathrm{3}}-b_{\mathrm{9}}{\dot{x}}_{\mathrm{3}}-

\mathrm{2}{x}_{\mathrm{3}}{b}_{\mathrm{11}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{x_{\mathrm{4}}\mathrm{\pi }}{\mathrm{180}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{x_{\mathrm{4}}\mathrm{\pi }}{\mathrm{180}}\right)\frac{\mathrm{\pi }}{\mathrm{180}}\left(-\frac{\mathrm{1}}{b_{\mathrm{12}}}{x}_{\mathrm{4}}\right);L_{\mathrm{21}}=b_{\mathrm{5}}{b}_{\mathrm{6}}\frac{\partial h(x_{\mathrm{1}},x_{\mathrm{2}})}{\partial x_{\mathrm{1}}}\left[\mathrm{1}+b_{\mathrm{7}}\left({\left(\frac{x_{\mathrm{2}}}{b_{\mathrm{3}}}\right)}^{b_{\mathrm{4}}}-\mathrm{1}\right)\right];L_{\mathrm{22}}=\mathrm{0};L_{\mathrm{31}}=\mathrm{0};L_{\mathrm{32}}=-\mathrm{2}{x}_{\mathrm{3}}{b}_{\mathrm{11}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{x_{\mathrm{4}}\mathrm{\pi }}{\mathrm{180}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{x_{\mathrm{4}}\mathrm{\pi }}{\mathrm{180}}\right)\frac{\mathrm{\pi }}{\mathrm{180}}\frac{\mathrm{1}}{b_{\mathrm{12}}}。

将公式简化为：

\mathit{x}=\left[\begin{array}{c}{x}_{\mathrm{2}}\\ x_{\mathrm{3}}\end{array}\right];\mathit{Q}\left(x,\dot{x}\right)=\left[\begin{array}{c}{Q}_{\mathrm{2}}\left(x,\dot{x}\right)\\ Q_{\mathrm{3}}\left(x,\dot{x}\right)\end{array}\right]；\mathit{L}(x,\dot{x})=\left[\begin{array}{cc}{L}_{\mathrm{21}}(x,\dot{x})& L_{\mathrm{22}}(x,\dot{x})\\ L_{\mathrm{31}}(x,\dot{x})& L_{\mathrm{32}}(x,\dot{x})\end{array}\right];\mathit{u}=\left[\begin{array}{c}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\end{array}\right]。

式（19）和（20）可以转变为：

\ddot{\mathit{x}}=\mathit{Q}(\mathit{x},\dot{\mathit{x}})+\mathit{L}(\mathit{x},\dot{\mathit{x}})\mathit{u}（21）

2.2　扩张状态观测器

由于模型与实际燃料电池系统间存在稳态偏差，采用扩展状态观测器（Extended state observe， ESO）来观测模型中的未知扰动和未建模动态。用d_{\mathrm{2}}和d_{\mathrm{3}}表示扰动和未建模动态，因此，燃料电池系统为：

\ddot{\mathit{x}}=\mathit{Q}(x,\dot{x})+\mathit{L}(x,\dot{x})\mathit{u}+\mathit{d}（22）

式中：\mathit{d}={\left[\begin{array}{cc}{d}_{\mathrm{2}}& d_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}。

令{\Theta }_{\mathrm{1},i}=x_{\mathrm{1},i}，{\Theta }_{\mathrm{2},i}=x_i和{\Theta }_{\mathrm{3},i}=d_i，i\in \{\mathrm{2,3}\}。假设\mathit{d}可微，式（22）表述为：

\left\{\begin{array}{l}{\dot{\Theta }}_{\mathrm{1},i}={\Theta }_{\mathrm{2},i}\\ {\dot{\Theta }}_{\mathrm{2},i}=Q_i\left(x,\dot{x}\right)+{\Theta }_{\mathrm{3},i}+\left[\begin{array}{cc}{L}_{i\mathrm{1}}& L_{i\mathrm{2}}\end{array}\right]{\left[u_{\mathrm{1}} u_{\mathrm{2}}\right]}^{\mathrm{T}}\\ {\dot{\Theta }}_{\mathrm{3},i}={\dot{d}}_i\end{array}\right.（23）

式（23）的ESO形式可以构建为：

\left\{\begin{array}{l}{\dot{\widehat{\Theta }}}_{\mathrm{1},i}={\widehat{\Theta }}_{\mathrm{2},i}+l_{\mathrm{1},i}\left({\Theta }_{\mathrm{1},i}-{\widehat{\Theta }}_{\mathrm{1},i}\right),\\ {\dot{\widehat{\Theta }}}_{\mathrm{2},i}=Q_i(x,\dot{x})+{\widehat{\Theta }}_{\mathrm{3},i}+\left[\begin{array}{cc}{L}_{i\mathrm{1}}& L_{i\mathrm{2}}\end{array}\right][u_{\mathrm{1}} u_{\mathrm{2}}{]}^{\mathrm{T}},\\ +l_{\mathrm{2},i}\left({\Theta }_{\mathrm{1},i}-{\widehat{\Theta }}_{\mathrm{1},i}\right)\\ {\dot{\widehat{\Theta }}}_{\mathrm{3},i}=l_{\mathrm{3},i}\left({\Theta }_{\mathrm{1},i}-{\widehat{\Theta }}_{\mathrm{1},i}\right)\end{array}\right.（24）

式中：l_{j,i}为ESO的增益参数，j\in \{\mathrm{1,2},\mathrm{3}\}，i\in \{\mathrm{2,3}\}。

选择增益为：

[l_{\mathrm{1},i} l_{\mathrm{2},i} l_{\mathrm{3},i}]=[{\omega }_{\mathrm{0},i}{\alpha }_{\mathrm{1}} {\omega }_{\mathrm{0},i}^{\mathrm{2}}{\alpha }_{\mathrm{2}} {\omega }_{\mathrm{0},i}^{\mathrm{3}}{\alpha }_{\mathrm{3}}]（25）

式中：{\omega }_{\mathrm{0},i}>\mathrm{0}，i\in \{\mathrm{2,3}\}；选择{\alpha }_j使特征多项式s^{\mathrm{3}}+{\alpha }_{\mathrm{1}}{s}^{\mathrm{2}}+{\alpha }_{\mathrm{2}}s+{\alpha }_{\mathrm{3}}是赫尔维茨矩阵，j\in\{\mathrm{1,2},\mathrm{3}\}，则{\alpha }_{\mathrm{1}}=\mathrm{3}，{\alpha }_{\mathrm{2}}=\mathrm{3}，{\alpha }_{\mathrm{3}}=\mathrm{1}。

令{\tilde{\Theta }}_{j,i}\left(t\right)={\Theta }_{j,i}\left(t\right)-{\widehat{\Theta }}_{j,i}\left(t\right)，ESO的误差为：

\left\{\begin{array}{l}{\dot{\tilde{\Theta }}}_{\mathrm{1},i}={\tilde{\Theta }}_{\mathrm{2},i}-\mathrm{3}{\omega }_{\mathrm{0},i}{\tilde{\Theta }}_{\mathrm{1},i}\\ {\dot{\tilde{\Theta }}}_{\mathrm{2},i}={\tilde{\Theta }}_{\mathrm{3},i}-\mathrm{3}{\omega }_{\mathrm{0},i}^{\mathrm{2}}{\tilde{\Theta }}_{\mathrm{1},i}\\ {\dot{\tilde{\Theta }}}_{\mathrm{3},i}={\dot{\Theta }}_{\mathrm{3},i}-{\omega }_{\mathrm{0},i}^{\mathrm{3}}{\tilde{\Theta }}_{\mathrm{1},i}\end{array}\right.（26）

令{\epsilon }_{j,i}\left(t\right)={\tilde{\Theta }}_{j,i}\left(t\right)/{\omega }_{\mathrm{0},i}^{j-\mathrm{1}}，式（26）可以写为：

{\dot{\epsilon }}_i={\omega }_{\mathrm{0},i}\mathit{A}{\epsilon }_i+\mathit{B}\frac{\dot{d}}{{\omega }_{\mathrm{0},i}^{\mathrm{2}}}（27）

式中：\mathit{A}=\left[\begin{array}{ccc}-\mathrm{3}& \mathrm{1}& \mathrm{0}\\ -\mathrm{3}& \mathrm{0}& \mathrm{1}\\ -\mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right],\mathit{B}={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}。

定理 1　假设d_i是有界的，存在一个常数{\sigma }_i>\mathrm{0}，有限时间T_{\mathrm{1}}>\mathrm{0}，可以得出\left|{\tilde{\Theta }}_i\right|\le {\sigma }_i，i\in \{\mathrm{2,3}\}，\forall t\ge T_{\mathrm{1}}>\mathrm{0}和{\omega }_{\mathrm{0}}>\mathrm{0}。因此{\sigma }_i=\mathit{O}(\mathrm{1}/{\omega }_{\mathrm{0},i}^k)，k\in \{\mathrm{1,2},\mathrm{3}\}^［17］。

证明 1　求解式（27），可以得到：

\mathit{\epsilon }\left(t\right)={\mathit{e}}^{{\omega }_{\mathrm{0},i}\mathit{A}t}\epsilon \left(\mathrm{0}\right)+{\int }_{\mathrm{0}}^t{\mathit{e}}^{{\omega }_{\mathrm{0},i}\mathit{A}(t-\tau )}\mathit{B}\frac{\mathit{d}}{{\omega }_{\mathrm{0},i}^{\mathrm{2}}}\mathrm{d}\tau（28）

令：

\mathit{p}\left(t\right)={\int }_{\mathrm{0}}^t{\mathit{e}}^{{\mathrm{\omega }}_{\mathrm{0},i}\mathit{A}(t-\mathrm{\tau })}\mathit{B}\frac{\mathit{d}}{{\mathrm{\omega }}_{\mathrm{0},i}^{\mathrm{2}}}\mathrm{d}\mathrm{\tau }（29）

由于\mathit{d}是有界的，因此\left|d_i\right|\le \mathrm{\delta }，i\in \{\mathrm{2,3}\}，其中{\mathrm{\delta }}_i是一个正的常数，可得：

\left|p_j\left(t\right)\right|\le \frac{\mathrm{\delta }}{{\mathrm{\omega }}_{\mathrm{0},i}^{\mathrm{3}}}\left[\left|{\left({\mathit{A}}^{-\mathrm{1}}\mathit{B}\right)}_i\right|+\left|{\left({\mathit{A}}^{-\mathrm{1}}{\mathit{e}}^{{\mathrm{\omega }}_{\mathrm{0},i}\mathit{A}t}\mathit{B}\right)}_i\right|\right]（30）

式中：j=\mathrm{1},\mathrm{2},\mathrm{3}，i\in \{\mathrm{2,3}\}，由于{\mathit{A}}^{-\mathrm{1}}=\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& -\mathrm{1}\\ \mathrm{1}& \mathrm{0}& -\mathrm{3}\\ \mathrm{0}& \mathrm{1}& -\mathrm{3}\end{array}\right]，可得：

\left|{\left({\mathit{A}}^{-\mathrm{1}}\mathit{B}\right)}_j\right|=\left\{\begin{array}{c}\mathrm{1}{\left|\right.}_{j=\mathrm{1}}\\ \mathrm{3}{\left|\right.}_{j=\mathrm{2},\mathrm{3}}\end{array}\right.（31）

由于\mathit{A}是赫尔维茨矩阵，因此，存在一个有限时间T_{\mathrm{1}}>\mathrm{0}使得：

\left|{\left[{\mathit{e}}^{{\mathrm{\omega }}_{\mathrm{0},i}\mathit{A}t}\right]}_{mn}\right|\le \frac{\mathrm{1}}{{\mathrm{\omega }}_{o,i}^{\mathrm{3}}}（32）

对于所有的t\ge T_{\mathrm{1}}，m,n=\mathrm{1},\mathrm{2},\mathrm{3}，可得：

\left|{\left[{\mathit{e}}^{{\mathrm{\omega }}_{\mathrm{0},i}\mathit{A}t}\mathit{B}\right]}_m\right|\le \frac{\mathrm{1}}{{\mathrm{\omega }}_{\mathrm{0},i}^{\mathrm{3}}}（33）

对于所有的t\ge T_{\mathrm{1}}，m=\mathrm{1},\mathrm{2},\mathrm{3}，T_{\mathrm{1}}取决于{\mathrm{\omega }}_{\mathrm{0},i}A。

令{\mathit{A}}^{-\mathrm{1}}=\left[\begin{array}{ccc}{s}_{\mathrm{11}}& s_{\mathrm{12}}& s_{\mathrm{13}}\\ s_{\mathrm{21}}& s_{\mathrm{22}}& s_{\mathrm{23}}\\ s_{\mathrm{31}}& s_{\mathrm{32}}& s_{\mathrm{33}}\end{array}\right],{\mathit{e}}^{{\omega }_{\mathrm{0},i}\mathit{A}t}=\left[\begin{array}{ccc}{w}_{\mathrm{11}}& w_{\mathrm{12}}& w_{\mathrm{13}}\\ w_{\mathrm{21}}& w_{\mathrm{22}}& w_{\mathrm{23}}\\ w_{\mathrm{31}}& w_{\mathrm{32}}& w_{\mathrm{33}}\end{array}\right],

可得：

\left|{\left({\mathit{A}}^{-\mathrm{1}}{\mathit{e}}^{{\mathrm{\omega }}_{\mathrm{0},i}\mathit{A}t}\mathit{B}\right)}_m\right|=|s_{m\mathrm{1}}{w}_{\mathrm{13}}+s_{m\mathrm{2}}{w}_{\mathrm{23}}+s_{m\mathrm{3}}{w}_{\mathrm{33}}|\le

\left\{\begin{array}{c}\frac{\mathrm{1}}{{\mathrm{\omega }}_{o,i}^{\mathrm{3}}}{\left|\right.}_{m=\mathrm{1}}\\ \frac{\mathrm{4}}{{\mathrm{\omega }}_{o,i}^{\mathrm{3}}}{\left|\right.}_{m=\mathrm{2},\mathrm{3}}\end{array}\right.（34）

式中：t\ge T_{\mathrm{1}}。

从式（30）（31）（34）可得：

|p_m\left(t\right)|\le \frac{\mathrm{3}\mathrm{\delta }}{{\mathrm{\omega }}_{\mathrm{0},i}^{\mathrm{3}}}+\frac{\mathrm{4}\mathrm{\delta }}{{\mathrm{\omega }}_{\mathrm{0},i}^{\mathrm{6}}}（35）

令{\epsilon }_{\mathrm{s}\mathrm{u}\mathrm{m}}\left(\mathrm{0}\right)=\left|{\epsilon }_{\mathrm{1}}\left(\mathrm{0}\right)\right|+\left|{\epsilon }_{\mathrm{2}}\left(\mathrm{0}\right)\right|+\left|{\epsilon }_{\mathrm{3}}\left(\mathrm{0}\right)\right|，可得：

\left|{\left[{\mathit{e}}^{{\omega }_{\mathrm{0},i}\mathit{A}t}\epsilon \left(\mathrm{0}\right)\right]}_m\right|\le \frac{{\epsilon }_{\mathrm{s}\mathrm{u}\mathrm{m}}\left(\mathrm{0}\right)}{{\omega }_{o,i}^{\mathrm{3}}}（36）

从式（28）可得：

\left|{\epsilon }_m\left(t\right)\right|\le \left|{\left[e^{{\omega }_{\mathrm{0},i}At}\epsilon \left(\mathrm{0}\right)\right]}_m\right|+\left|p_m\left(t\right)\right|（37）

令{\tilde{\Theta }}_{\mathrm{s}\mathrm{u}\mathrm{m}}\left(\mathrm{0}\right)=\left|{\tilde{\Theta }}_{\mathrm{1}}\left(\mathrm{0}\right)\right|+\left|{\tilde{\Theta }}_{\mathrm{2}}\left(\mathrm{0}\right)\right|+\left|{\tilde{\Theta }}_{\mathrm{3}}\left(\mathrm{0}\right)\right|，由式（35）~（37）和{\epsilon }_{j,i}\left(t\right)={\tilde{\Theta }}_{j,i}\left(t\right)/{\omega }_{\mathrm{0},i}^{j-\mathrm{1}}可得：

\left|{\tilde{\Theta }}_m\left(t\right)\right|\le \left|\frac{{\tilde{\Theta }}_{\mathrm{s}\mathrm{u}\mathrm{m}}\left(\mathrm{0}\right)}{{\mathrm{\omega }}_{o,i}^{\mathrm{3}}}\right|+\frac{\mathrm{3}\mathrm{\delta }}{{\mathrm{\omega }}_{o,i}^{\mathrm{4}-i}}+\frac{\mathrm{4}\mathrm{\delta }}{{\mathrm{\omega }}_{o,i}^{\mathrm{7}-i}}={\mathrm{\sigma }}_m（38）

式中：t\ge T_{\mathrm{1}}，m=\mathrm{1},\mathrm{2},\mathrm{3}。

证明完毕。

2.3　控制器的设计

定义{\mathrm{y}}_i^{\mathrm{*}}，i\in \{\mathrm{1,2}\}为参考输出信号，输出误差为：

\mathit{z}=\left[\begin{array}{c}{z}_{\mathrm{2}}\\ z_{\mathrm{3}}\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{2}}-x_{\mathrm{2}}^{\mathrm{*}}\\ x_{\mathrm{3}}-x_{\mathrm{3}}^{\mathrm{*}}\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{2}}-\left(\frac{y_{\mathrm{1}}^{\mathrm{*}}}{b_{\mathrm{8}}}+y_{\mathrm{2}}^{\mathrm{*}}\right)\\ x_{\mathrm{3}}-y_{\mathrm{2}}^{\mathrm{*}}\end{array}\right]（39）

\mathit{Q}(x,\dot{x})和\mathit{L}(x,\dot{x})转化为关于\mathit{z}和t的函数^［17］：

\left\{\begin{array}{l}\mathit{Q}(x,\dot{x})=\mathit{Q}(\mathit{z},\dot{\mathit{z}},t)\\ \mathit{L}(x,\dot{x})=\mathit{L}(\mathit{z},\dot{\mathit{z}},t)\end{array}\right.（40）

因此，系统（22）转化为：

\ddot{\mathit{z}}=\mathit{Q}(\mathit{z},\dot{\mathit{z}},t)+\mathit{L}(\mathit{z},\dot{\mathit{z}},t)\mathit{u}+\mathit{d}-({\mathit{x}}^{\mathrm{*}}{)}^{(\mathrm{2})}（41）

当系统模型为式（41）的HOFA系统时，可采用以下控制律使系统稳定^［18］：

\left\{\begin{array}{c}\mathit{u}=-{\mathit{L}}^{-\mathrm{1}}(\mathit{z},\dot{\mathit{z}},t)({\mathit{A}}^{\mathrm{0}~\mathrm{1}}{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+{\mathit{u}}^{\mathrm{*}})\\ {\mathit{u}}^{\mathrm{*}}=\widehat{\mathit{d}}+\mathit{Q}(\mathit{z},\dot{\mathit{z}},t)-({\mathit{x}}^{\mathrm{*}}{)}^{(\mathrm{2})}\end{array}\right.（42）

得到下述形式的闭环系统：

{\mathit{z}}^{(\mathrm{2})}+{\mathit{A}}^{\mathrm{0}~\mathrm{1}}{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}=\tilde{\mathit{d}}（43）

将控制律代入燃料电池系统中，确保：①{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\to \mathrm{0}，t\to \mathrm{\infty }；②\Vert \mathit{d}-\widehat{\mathit{d}}\Vert是有界的。

2.4　稳定性分析

定理 2　假设\mathit{A}\in {\mathbb{R}}^{n\times n}满足^［19］：

\mathrm{R}\mathrm{e}{\lambda }_i\left(\mathit{A}\right)\le -\frac{\gamma }{\mathrm{2}},i=\mathrm{1,2},\dots ,n（44）

当\gamma >\mathrm{0}，则存在正定矩阵\mathit{P}\in {\mathbb{R}}^{n\times n}，满足：

{\mathit{\Phi }}^{\mathrm{T}}({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})\mathit{P}({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})+

\mathit{P}({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})\mathit{\Phi }({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})\le -\mu \mathit{P}（45）

其中，\mathit{P}\left({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}}\right)为：

\mathit{P}({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})=[{\mathit{P}}_{\mathrm{1}}{\mathit{P}}_{\mathrm{2}}\cdots {\mathit{P}}_n], {\mathit{P}}_i\in {\mathbb{R}}^{nr\times r}（46）

为更方便分析控制器的稳定性，引入以下符号：

{\mathit{P}}_L({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})\triangleq \mathit{P}({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}})\left[\begin{array}{c}\mathrm{0}\\ {\mathit{I}}_r\end{array}\right]={\mathit{P}}_n（47）

闭环系统（43）表示为如下状态空间形式：

{\dot{\mathit{z}}}^{(\mathrm{0}~\mathrm{1})}=\mathit{\Phi }({\mathit{A}}^{\mathrm{0}~\mathrm{1}}){\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+\left[\begin{array}{c}{\mathrm{0}}_{\mathrm{2}\times \mathrm{1}}\\ \tilde{\mathit{d}}\end{array}\right]（48）

证明 2　对于系统（48），可选择如下Lyapunov函数：V={\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\mathit{P}({\mathit{A}}^{\mathrm{0}~\mathrm{1}}){\mathit{z}}^{(\mathrm{0}~\mathrm{1})}，则：

\begin{array}{l}\dot{V}={\left({\dot{\mathit{z}}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\mathit{P}{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+{\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\mathit{P}{\dot{\mathit{z}}}^{(\mathrm{0}~\mathrm{1})}=\\ {\left(\mathit{\Phi }{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+\left[\begin{array}{c}{\mathrm{0}}_{\mathrm{2}\times \mathrm{1}}\\ \tilde{\mathit{d}}\end{array}\right]\right)}^{\mathrm{T}}\mathit{P}{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+{\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\\ \mathit{P}\left(\mathit{\Phi }{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+\left[\begin{array}{c}{\mathrm{0}}_{\mathrm{2}\times \mathrm{1}}\\ \tilde{\mathit{d}}\end{array}\right]\right)={\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\left({\mathit{\Phi }}^{\mathrm{T}}\mathit{P}+\mathit{P}\mathit{\Phi }\right){\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+\\ \mathrm{2}{\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\mathit{P}\left[\begin{array}{c}{\mathrm{0}}_{\mathrm{2}\times \mathrm{1}}\\ \tilde{\mathit{d}}\end{array}\right]\le -\mu {\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}\mathit{P}{\mathit{z}}^{(\mathrm{0}~\mathrm{1})}+\\ \mathrm{2}{\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}{\mathit{P}}_L\left({\mathit{A}}^{\mathrm{0}~\mathrm{1}}\right)\tilde{\mathit{d}}\le -\mu V+\mathrm{2}{\left({\mathit{z}}^{(\mathrm{0}~\mathrm{1})}\right)}^{\mathrm{T}}{\mathit{P}}_L\left({\mathit{A}}^{\mathrm{0}~\mathrm{1}}\right)\\ \tilde{\mathit{d}}\le -\mu V+\mathrm{2}{\sigma }_m‖{\mathit{P}}_L^{\mathrm{T}}({\mathit{A}}^{\mathrm{0}~\mathrm{1}}){\mathit{z}}^{\mathrm{0}~\mathrm{1}}‖\le -\mu V+\epsilon \mathrm{。}\end{array}

根据比较定理^［18］，可得：

V\le V(\mathrm{0})e^{-\mu t}+\frac{\epsilon }{\mu }\left(\mathrm{1}-e^{-\mu t}\right)（49）

因此：

V\le \left(V(\mathrm{0})-\frac{\epsilon }{\mu }\right)e^{-\mu t}+\frac{\epsilon }{\mu }\to \frac{\epsilon }{\mu }, t\to \mathrm{\infty }（50）

证明完毕。

2.5　参数化设计

命题 1对于任意选择的矩阵\mathit{F}\in {\mathbb{R}}^{nr\times nr}，矩阵{\mathit{A}}^{\mathrm{0}~n-\mathrm{1}}和非奇异矩阵\mathit{V}\in {\mathbb{R}}^{nr\times nr}满足^［15］：\mathit{\Phi }\left({\mathit{A}}^{\mathrm{0}~n-\mathrm{1}}\right)=\mathbf{V}\mathbf{F}{\mathit{V}}^{-\mathrm{1}}，其中：{\mathit{A}}^{\mathrm{0}~n-\mathrm{1}}=-\mathit{Z}{\mathit{F}}^n{\mathit{V}}^{-\mathrm{1}}\left(\mathit{Z},\mathit{F}\right)，\mathit{V}=\mathit{V}\left(\mathit{Z},\mathit{F}\right)=\left[\begin{array}{c}\mathit{Z}\\ \mathit{Z}\mathit{F}\\ ⋮\\ \mathit{Z}{\mathit{F}}^{\mathit{n}-\mathrm{1}}\end{array}\right]，其中：\mathit{Z}\in {\mathbb{R}}^{r\times nr}为任意参数矩阵，满足\mathrm{d}\mathrm{e}\mathrm{t}\mathit{V}\left(\mathit{Z},\mathit{F}\right)\ne \mathrm{0}。\mathit{F}为对角元素为负的对角矩阵。

\mathit{F}=\left[\begin{array}{cccc}-a& -b& \mathrm{0}& \mathrm{0}\\ b& -a& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -c& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -d\end{array}\right]（51）

式中：a、b、c和d是正标量。本文选择：

\mathit{F}=\left[\begin{array}{cccc}-\mathrm{10}& -\mathrm{10}& \mathrm{0}& \mathrm{0}\\ \mathrm{10}& -\mathrm{10}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{30}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{30}\end{array}\right],\mathit{Z}=\left[\begin{array}{c}\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}\end{array}\\ \begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\end{array}\right]\mathrm{。}因此，

\mathit{V}\left(\mathit{Z},\mathit{F}\right)=\left[\begin{array}{c}\mathit{Z}\\ \mathit{Z}\mathit{F}\end{array}\right]=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& -\mathrm{20}& -\mathrm{30}& -\mathrm{30}\\ -\mathrm{10}& -\mathrm{10}& \mathrm{0}& -\mathrm{30}\end{array}\right]，其中，满足\mathrm{d}\mathrm{e}\mathrm{t}\mathit{V}\left(\mathit{Z},\mathit{F}\right)\ne \mathrm{0}，可得：{\mathit{A}}^{\mathrm{0}~\mathrm{1}}=-\mathit{Z}{\mathit{F}}^{\mathrm{2}}{\mathit{V}}^{-\mathrm{1}}\left(\mathit{Z},\mathit{F}\right)=\left[\begin{array}{cccc}\mathrm{600}& -\mathrm{600}& \mathrm{50}& -\mathrm{20}\\ \mathrm{300}& \mathrm{0}& \mathrm{10}& \mathrm{30}\end{array}\right]。

2.6　输出参考值

进气流量的参考值W_{\mathrm{s}\mathrm{m}}^{\mathrm{*}}为与负载电流和最佳氧过剩比{\lambda }_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{*}}有关的函数。

W_{\mathrm{s}\mathrm{m}}^{\mathrm{*}}=\frac{\mathrm{1}}{x_{{\mathrm{O}}_{\mathrm{2}}}}\frac{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}{I}_{\mathrm{s}\mathrm{t}}}{\mathrm{4}F}{\lambda }_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{*}}（52）

式中：x_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{0.21}为空气中氧的质量分数；{\lambda }_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{*}}为^［20］：

{\lambda }_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{*}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}{I}_{\mathrm{s}\mathrm{t}}^{\mathrm{3}}-\mathrm{2.87}\times {\mathrm{10}}^{-\mathrm{5}}{I}_{\mathrm{s}\mathrm{t}}^{\mathrm{2}}+\mathrm{2.23}\times

{\mathrm{10}}^{-\mathrm{3}}{I}_{\mathrm{s}\mathrm{t}}+\mathrm{2.5}（53）

阴极压力的参考值P_{\mathrm{c}\mathrm{a}}^{\mathrm{*}}通常随阳极压力变化，并与负载电流呈正相关^［7］。

P_{\mathrm{c}\mathrm{a}}^{\mathrm{*}}=\mathrm{0.015}\mathrm{42}\times I_{\mathrm{s}\mathrm{t}}^{\mathrm{3}}-\mathrm{10.25}\times I_{\mathrm{s}\mathrm{t}}^{\mathrm{2}}+

\mathrm{2}\mathrm{327}\times I_{\mathrm{s}\mathrm{t}}-\mathrm{28}\mathrm{240}（54）

为了得到曲线光滑的{\dot{y}}_i^{\mathrm{*}}和{\ddot{y}}_i^{\mathrm{*}}，i\in \left\{\mathrm{1,2}\right\}，使用两个二阶低通滤波器处理信号。原始参考信号{\overline{y}}_i^{\mathrm{*}}通过滤波器产生。

\frac{y_i^{\mathrm{*}}}{{\overline{y}}_i^{\mathrm{*}}}=\frac{{\omega }_{n,i}^{\mathrm{2}}}{s^{\mathrm{2}}+\mathrm{2}{\zeta }_i{\omega }_{n,i}s+{\omega }_{n,i}^{\mathrm{2}}}（55）

燃料电池高阶全驱系统的控制策略总体结构如图2所示，将ESO与HOFA系统方法相结合，基于全驱系统的特性，消除燃料电池模型中的非线性部分，建立起新的闭环动态特性。

3 仿真与实验结果

3.1　仿真结果

通过一系列的仿真验证本文所提质子交换膜燃料电池HOFA控制方法的有效性。燃料电池的负载电流如图3所示，负载电流在120~280 A不断阶跃变化。本文将提出的控制策略（HOFA+ESO）与文献中基于扩展状态观测器的多输出反馈控制策略方法^［16］（MFB+ESO）和文献中的反馈线性化控制方法^［21］（FBL）进行比较，以验证HOFA+ESO方法的控制效果。在与文献［21］控制方法进行比较时，为体现仿真验证的严谨性，本文将该文献所提出的方法结合扩张状态观测器来补偿干扰和未建模动态。同时，由于该方法是单输出的控制，因此，本文在该方法的基础上改进为多输出控制。

3种不同控制策略下得到的电压和节气门开度如图4所示。在负载电流快速阶跃变化的情况下，基于HOFA+ESO的策略中，电压最大超调量为33%，最小超调量为4%，节气门最大超调量为48%，最小超调量为12%。电压达到稳定输出的响应速度为1 100 ms左右，节气门达到稳定输出的响应速度为400 ms左右。基于MFB+ESO的策略中，电压最大超调量为97%，最小超调量为15%。节气门最大超调量为87%，最小超调量为18%。电压达到稳定输出的响应速度为1 100 ms左右，节气门达到稳定输出的响应速度为600 ms左右。基于FBL+ESO的策略中，电压最大超调量为92%，最小超调量为13%。节气门最大超调量为79%，最小超调量为16%。电压达到稳定输出的响应速度为1 100 ms左右，节气门达到稳定输出的响应速度为600 ms左右。

图5为不同控制策略下的仿真结果。仿真结果表明在负载电流快速阶跃变化的情况下，基于HOFA+ESO的控制策略中，阴极压力的响应速度为1 000 ms， 进气流量的响应速度为600 ms，控制误差均在1%以内。基于MFB+ESO和FBL+ESO的策略中，阴极压力的响应速度为1 000 ms，进气流量的响应速度为800 ms，控制误差均在1%以内。结果表明：HOFA+ESO策略动态响应快，而且曲线平滑无超调。其原因为，与传统的非线性方法相比，HOFA系统方法通过提供更大的设计自由度，使其超调量和响应速度更快，通过设计不同的a、b、c、d值，使其能达到更好的控制效果^［22］。当a、b、c、d过大时，系统的超调量会变大；当a、b、c、d过小时，跟踪参考值的误差会变大，因此，选择a=b=10，c=d=30，使超调量和跟踪误差在合理的范围内。

3.2　实验结果

实验平台如图6所示，主要由空压机、中冷器、空气缓冲罐、节流阀、ECU、空压机控制器等组成。其主要目的是验证阴极进气系统的工作状态。该实验平台利用传感器实时监测关键参数，包括进气流量、阴极压力和温度。

本文实验验证进气歧管流量与阴极压力的控制效果。在实验验证中，阴极压力的参考值固定为104 kPa。由图7（b）可以看出，当进气歧管流量从40→45→50 g/s变化时，压缩机转速增大，导致进气歧管流量增大。由于进气流量与阴极压力之间的耦合关系，进气歧管流量的增加导致阴极压力升高。因此，增大节气门开角，使此时压力稳定。 同样，进气歧管流量从50→25 g/s变化时，空压机转速降低，节气门开度也减小。在图7（a）（b）中，进气歧管流量和阴极压力在参考值附近波动。进气歧管流量在±0.5 g/s范围内波动和阴极压力在±0.5 kPa范围内波动是由传感器的测量误差引起的。

4 结束语

本文针对质子交换膜燃料电池进气系统进气流量和阴极压力难以精确协同控制问题，提出了一种基于HOFA系统结合扩张状态观测器的控制方案。并通过仿真和实验验证了所提HOFA方法的有效性。仿真结果表明：在负载电流快速阶跃变化情况下，基于HOFA+ESO控制策略，电压达到稳定输出的响应速度为1 200 ms左右，节气门达到稳定输出的响应速度为400 ms左右。阴极压力的响应速度为1 000 ms， 进气流量的响应速度为600 ms，控制误差均在1%以内。该方法能够在负载电流变化情况下快速准确地跟踪到参考输出。但是，阴极压力对负载电流阶跃变化的响应时间还有不足，在后续的研究中，将对阴极的响应速度进行调节以实现更精确、快速的控制。

参考文献

[1] Ehsani M, Gao Y, Miller J M. Hybrid electric vehicles: architecture and motor drives[J]. Proceedings of the IEEE, 2007, 95(4): 719-728.

[2] Ma J, Liu X, Zou X, et al. Degradation prognosis for proton exchange membrane fuel cell based on hybrid transfer learning and intercell differences[J]. ISA Transactions, 2021, 113: 149-165.

[3] Wang Y, Chen K S, Mishler J, et al. A review of polymer electrolyte membrane fuel cells: technology, applications, and needs on fundamental research[J]. Applied Energy, 2011, 88(4): 981-1007.

[4] Zhang D, Baraldi P, Cadet C, et al. An ensemble of models for integrating dependent sources of information for the prognosis of the remaining useful life of proton exchange membrane fuel cells[J]. Mechanical Systems and Signal Processing, 2019, 124: 479-501.

[5] Wang G L, Wang Y, Shi J H, et al. Coordinating IMC-PID and adaptive SMC controllers for a PEMFC[J]. ISA Transactions, 2010, 49(1): 87-94.

[6] 高金武, 王义琳, 刘华洋, 等. 基于滑模观测器的质子交换膜燃料电池阴极进气系统解耦控制[J]. 吉林大学学报: 工学版, 2022, 52(9): 2156-2167.

[7] Gao Jin-wu, Wang Yi-lin, Liu Hua-yang, et al. Decoupling control for proton exchange membrane fuel cell air supply system based on sliding mode observer[J]. Journal of Jilin University (Engineering and Technology Edition), 2022, 52(9): 2156-2167.

[8] Aldurra A, Yurkovich S, Guezennec Y. Study of nonlinear control schemes for an automotive traction PEM fuel cell system[J]. International Journal of Hydrogen Energy, 2010, 35(20): 11291-11307.

[9] 胡云峰, 于彤, 杨惠策, 等. 低温环境下燃料电池启动优化控制方法[J]. 吉林大学学报: 工学版, 2022, 52(9): 2034-2043.

[10] Hu Yun-feng, Yu Tong, Yang Hui-ce, et al. Optimal control method of fuel cell start⁃up in low temperature environment[J]. Journal of Jilin University (Engineering and Technology Edition), 2022, 52(9): 2034-2043.

[11] Danzer M A, Wilhelm J, Aschemann H, et al. Model-based control of cathode pressure and oxygen excess ratio of a PEM fuel cell system[J]. Journal of Power Sources, 2008, 176(2): 515-522.

[12] Pukrushpan J T, Stefanopoulou A G, Peng H. Modeling and Control of Fuel Cell Systems and Fuel Processors,Modeling, Analysis, and Feedback Design[M]. London:Springer,2004.

[13] Talj R, Ortega R, Astolfi A. Passivity and robust PI control of the air supply system of a PEM fuel cell model[J]. Automatica, 2011, 47(12): 2554-2561.

[14] Ma Y, Zhang F, Gao J, et al. Oxygen excess ratio control of PEM fuel cells using observer-based nonlinear triple-step controller[J]. International Journal of Hydrogen Energy, 2020, 45(54): 29705-29717.

[15] Liu J, Gao Y, Su X, et al. Disturbance-observer-based control for air management of PEM fuel cell systems via sliding mode technique[J]. IEEE Transactions on Control Systems Technology, 2018, 27(3): 1129-1138.

[16] Wang Y, Wang Y, Xu J, et al. Observer-based discrete adaptive neural network control for automotive PEMFC air-feed subsystem[J]. IEEE Transactions on Vehicular Technology, 2021, 70(4): 3149-3163.

[17] Duan G R. High-order fully actuated system approaches: part I. models and basic procedure[J]. International Journal of Systems Science, 2021, 52(2): 422-435.

[18] Li M, Yin H, Ding T, et al. Air flow rate and pressure control approach for the air supply subsystems in PEMFCs[J]. ISA Transactions, 2022, 128: 624-634.

[19] Zheng Q, Dong L, Lee D, et al. Active disturbance rejection control for mems gyroscopes[J]. Transactions on Control Systems Technology, 2009, 17(6): 1432-1438.

[20] Duan G R. High-order fully actuated system approaches: adaptive control and high-order backstep[J]. International Journal of Systems Science, 2021, 52(5): 972-989.

[21] Duan G. High-order fully actuated system approaches: robust control and high-order backstepping[J]. International Journal of Systems Science, 2021, 52(5): 952-971.

[22] Zhang H K, Wang Y F, Wang D H, et al. Adaptive robust control of oxygen excess ratio for PEMFC system based on type-2 fuzzy logic system[J]. Information Sciences, 2020, 511: 1-17.

[23] Chen J, Liu Z, Wang F, et al. Optimal oxygen excess ratio control for PEM fuel cells[J]. IEEE Transactions on Control Systems Technology, 2017, 26(5): 1711-1721.

[24] 段广仁. 高阶系统方法——I. 全驱系统与参数化设计[J]. 自动化学报, 2020, 46(7): 1333-1345.

[25] Duan Guang-ren. High order systems approach—I. Fully actuated system and parametric design[J]. Acta Automatica Sinica, 2020, 46(7): 1333-1345.

基金资助

汽车底盘集成与仿生全国重点实验室开放项目(20210215)