质子交换膜燃料电池过氧比LQR最优控制策略

摘要

建立了质子交换膜燃料电池空气供给系统六阶模型。将神经网络前馈控制分别结合PID控制及线性二次型调节器（LQR）控制，来比较两种复合控制策略在燃料电池变载情况下过氧比的快速响应情况。在LQR控制策略中，首先通过单一平衡点进行模型线性化以获取状态反馈增益及降维状态观测器增益。仿真结果表明，基于单平衡点线性化的LQR结合神经网络前馈的复合控制策略在全工况范围内明显优于神经网络前馈加PID的控制策略。针对偏离平衡点较远工况处控制效果有所下降的情况，进一步采用多平衡点线性化的方法，优化设计LQR以动态调整相应增益，结果表明多平衡点LQR控制方法在过氧比控制中展现出最佳的快速响应性和稳定性。

关键词

质子交换膜燃料电池; 线性二次型调节器; 多点平衡; 过氧比; 最优控制

质子交换膜燃料电池（PEMFC）是一种能量转化装置，具有高功率密度、高能量效率的优点，同时以其零污染、低噪音的优点，被广泛应用于车载动力能源装置^［1］。在车载条件下，由于质子交换膜燃料电池工作条件的频繁变化，即拉载电流的频繁变化，对燃料的供给尤其是空气侧的氧气供给提出了较高的要求。一方面，如果拉载电流要求较高，而阴极侧氧气流量过低，会导致三相反应界面反应物不足，造成阴极氧饥饿的现象，严重的氧饥饿会进一步的造成关键材料的腐蚀从而影响燃料电池的耐久性^［2］；另一方面，如果拉载电流要求降低，而阴极侧氧气流量过高，会造成聚合物膜快速干燥的现象，增加了离子在其中传递的阻力，除此之外，过高的氧气流量也意味着空压机的功耗过大，造成不必要的功率损耗，使整个系统的能量转换效率降低。因此，有必要对变载条件下的OER（过氧比）进行精准的控制，通过设计合适的控制策略，使得OER随变载电流变化而快速响应，从而提高燃料电池的效率和耐久性。

近年来，已有很多学者针对车用质子交换膜燃料电池过氧比的控制策略进行了研究。Pukrushpan等^［3］最先建立面向控制的燃料电池系统九阶模型，基于该模型提出了前馈加比例积分⁃微分（PID）控制及基于单工作平衡点线性化后的状态反馈控制方法。Niknezhadi等^［4］沿用了Pukrushpan^［3］提出的九阶模型，设计了线性二次型调节器（LQR）对过氧比进行控制，并成功在单电池台架上得到验证，实验结果较好，然而这种控制方法仅在工作平衡点附近，或电池功率较小，非线性程度不高的台架上较为有效，在大功率电堆系统的全局范围内表现不佳。Ou等^［5］设计了一个模糊前馈PID控制器（FFPID），来调整氧气过量比，避免燃料电池负载突然变化的情况下出现氧气不足。Li等^［6］设计了一种超扭曲滑膜控制器来控制过氧比，并提出了一种参数识别和控制器参数调整的方法。此外，模型预测控制（MPC）也因其滚动更新寻找全局最优的特性被广泛运用在燃料电池系统控制领域。Hahn等^［7］开发了一种基于模型预测控制的空气供给系统自适应控制策略，基于对不用工况下的仿真研究，证明了模型预测控制器的自适应和效率优化行为，然而由于模型预测控制在计算量和计算时间上的复杂性，很难在车载台架中得以实现。陈凤祥等^［8］基于多工作点的线性化模型结合多内膜控制技术将非线性模型中的过氧比控制在2.0~2.4，使系统的净功率输出最大，仿真结果表明其比单内膜具有更佳的全局控制性能。

本文基于建立的空气供给系统六阶控制模型，首先选取一个工作点对非线性系统状态空间方程线性化，采用LQR控制器结合神经网络前馈的控制方法对过氧比进行控制，将控制结果与神经网络前馈加PID控制结果进行比较；其次提出基于T-S模糊多模型方法的多平衡点LQR控制策略，选取扰动电流分别为100、140、180 A时的工作点对空气供给系统模型线性化，针对三个线性模型分别设计LQR控制器，通过扰动电流的隶属度函数决定控制器的输出；最后通过比较过氧比在多平衡点LQR控制策略下与单平衡点LQR控制策略下的响应性能验证本文提出的多平衡点LQR控制的优越性。

1 空气供给系统数学模型

基于参考文献［3］的工作，为了兼顾计算速度和模型准确度，本文建立了空气供给系统的六阶控制模型，其状态空间方程如下：

$${\dot{x}}_{\mathrm{1}}=\frac{x_{{\mathrm{O}}_{\mathrm{2}},\mathrm{c}\mathrm{a},\mathrm{i}\mathrm{n}}}{\mathrm{1}+{\Theta }_{\mathrm{a}\mathrm{t}\mathrm{m}}}{k}_{\mathrm{s}\mathrm{m},\mathrm{o}\mathrm{u}\mathrm{t}}(x_{\mathrm{5}}-P_{\mathrm{c}\mathrm{a}})-\frac{M_{{\mathrm{O}}_{\mathrm{2}}}{P}_{{\mathrm{O}}_{\mathrm{2}}}}{M_a{P}_a+M_v{P}_{\mathrm{s}\mathrm{a}\mathrm{t}}}{k}_{\mathrm{s}\mathrm{m},\mathrm{o}\mathrm{u}\mathrm{t}}(P_{\mathrm{c}\mathrm{a}}-x_{\mathrm{3}})-M_{{\mathrm{O}}_{\mathrm{2}}}\frac{n{I}_{\mathrm{s}\mathrm{t}}}{\mathrm{4}F}$$

$${\dot{x}}_{\mathrm{2}}=\frac{\mathrm{1}-x_{{\mathrm{O}}_{\mathrm{2}}}}{\mathrm{1}+{\Theta }_{\mathrm{a}\mathrm{t}\mathrm{m}}}{k}_{\mathrm{s}\mathrm{m},\mathrm{o}\mathrm{u}\mathrm{t}}(x_{\mathrm{5}}-P_{\mathrm{c}\mathrm{a}})-\frac{M_a{P}_a-M_{{\mathrm{O}}_{\mathrm{2}}}{P}_{{\mathrm{O}}_{\mathrm{2}}}}{M_a{P}_a+M_v{P}_{\mathrm{s}\mathrm{a}\mathrm{t}}}{k}_{\mathrm{c}\mathrm{a},\mathrm{o}\mathrm{u}\mathrm{t}}(P_{\mathrm{c}\mathrm{a}}-x_{\mathrm{3}})$$

$${\dot{x}}_{\mathrm{3}}=\frac{R_a{T}_{\mathrm{r}\mathrm{m}}}{V_{\mathrm{r}\mathrm{m}}}(k_{\mathrm{c}\mathrm{a},\mathrm{o}\mathrm{u}\mathrm{t}}(P_{\mathrm{c}\mathrm{a}}-x_{\mathrm{3}})-W_{\mathrm{r}\mathrm{m},\mathrm{o}\mathrm{u}\mathrm{t}})$$

$${\dot{x}}_{\mathrm{4}}=W_{\mathrm{c}\mathrm{p}}-k_{\mathrm{c}\mathrm{a},\mathrm{i}\mathrm{n}}(x_{\mathrm{5}}-P_{\mathrm{c}\mathrm{a}})$$

$${\dot{x}}_{\mathrm{5}}=\frac{\gamma R_a}{V_{\mathrm{s}\mathrm{m}}}(W_{\mathrm{c}\mathrm{p}}(T_{\mathrm{a}\mathrm{t}\mathrm{m}}+\frac{T_{\mathrm{a}\mathrm{t}\mathrm{m}}}{{\eta }_{\mathrm{c}\mathrm{p}}}({(\frac{x_{\mathrm{5}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}})}^{\frac{\gamma -\mathrm{1}}{\gamma }}-\mathrm{1}))-k_{\mathrm{c}\mathrm{a},\mathrm{i}\mathrm{n}}(x_{\mathrm{5}}-P_{\mathrm{c}\mathrm{a}})\frac{x_{\mathrm{5}}{V}_{\mathrm{s}\mathrm{m}}{M}_a}{R{x}_{\mathrm{4}}})$$

$${\dot{x}}_{\mathrm{6}}=-\frac{{\eta }_{\mathrm{c}\mathrm{m}}{k}_v{k}_t}{J_{\mathrm{c}\mathrm{p}}{R}_{\mathrm{c}\mathrm{m}}}{\omega }_{\mathrm{c}\mathrm{p}}-\frac{C_p{T}_{\mathrm{a}\mathrm{t}\mathrm{m}}}{{\eta }_{\mathrm{c}\mathrm{p}}}\frac{{(\frac{x_{\mathrm{5}}}{P_{\mathrm{a}\mathrm{t}\mathrm{m}}})}^{\frac{\gamma -\mathrm{1}}{\gamma }}}{J_{\mathrm{c}\mathrm{p}}{\omega }_{\mathrm{c}\mathrm{p}}}{W}_{\mathrm{c}\mathrm{p}}+\frac{{\eta }_{\mathrm{c}\mathrm{m}}{k}_t}{J_{\mathrm{c}\mathrm{p}}{R}_{\mathrm{c}\mathrm{m}}}{U}_{\mathrm{c}\mathrm{m}}$$

$$z=\frac{\frac{x_{{\mathrm{O}}_{\mathrm{2}},\mathrm{c}\mathrm{a},\mathrm{i}\mathrm{n}}}{\mathrm{1}+\Theta \left(T_{\mathrm{s}\mathrm{t}}\right)}{k}_{\mathrm{c}\mathrm{a},\mathrm{i}\mathrm{n}}(x_{\mathrm{5}}-P_{\mathrm{c}\mathrm{a}})}{M_{{\mathrm{O}}_{\mathrm{2}}}\frac{n{I}_{\mathrm{s}\mathrm{t}}}{\mathrm{4}F}}$$

（1）

式中：${\left[x_{\mathrm{1}}{x}_{\mathrm{2}}{x}_{\mathrm{3}}{x}_{\mathrm{4}}{x}_{\mathrm{5}}{x}_{\mathrm{6}}\right]}^{\mathrm{T}}={\left[m_{{\mathrm{O}}_{\mathrm{2}}}{m}_{{\mathrm{N}}_{\mathrm{2}}}{P}_{\mathrm{r}\mathrm{m}}{m}_{\mathrm{s}\mathrm{m}}{P}_{\mathrm{s}\mathrm{m}}{\omega }_{\mathrm{c}\mathrm{p}}\right]}^{\mathrm{T}}$；$z={\lambda }_{{\mathrm{O}}_{\mathrm{2}}}$为氧气过量比； $x_{{\mathrm{O}}_{\mathrm{2}},\mathrm{c}\mathrm{a},\mathrm{i}\mathrm{n}}$为阴极入口处氧气质量分数；$P_{\mathrm{c}\mathrm{a}}$为阴极侧压力，Pa；${\Theta }_{\mathrm{a}\mathrm{t}\mathrm{m}}$为标准大气压下比湿度；$W_{\mathrm{c}\mathrm{p}}$为空压机出口空气质量流量，$\mathrm{k}\mathrm{g}\cdot s^{-\mathrm{1}}$；${\omega }_{\mathrm{c}\mathrm{p}}$为空压机转速，$\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}$；$I_{\mathrm{s}\mathrm{t}}$为电流，A；$U_{\mathrm{c}\mathrm{m}}$为空压机电机控制电压，V；其余均为常数。由于燃料电池在运行过程中始终有液态水产生，因此本模型设定为阴极出口的水蒸气分压力始终处于饱和状态。本文基于该模型进行过氧比控制器的设计和仿真分析，将上述六阶空气供给系统非线性模型线性化，将空压机出口空气质量流量$W_{\mathrm{c}\mathrm{p}}$，供气管路压力$P_{\mathrm{s}\mathrm{m}}$，下游管路压力$P_{\mathrm{r}\mathrm{m}}$作为可测量输出，建立降维状态观测器，在此基础上设计LQR控制器实现状态反馈控制，控制输入为空压机电压$U_{\mathrm{c}\mathrm{m}}$，控制目标为将过氧比控制在${\lambda }_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{2}$。

2 单平衡点LQR控制器设计

2.1　模型线性化

六阶状态空间方程的完整表达式如下：

$$\left\{\begin{array}{c}\dot{x}=f(x,u,w)\\ y=g(x)\\ z=h(x,w)\end{array}\right.$$

（2）

式中： ${\left[y_{\mathrm{1}}{y}_{\mathrm{2}}{y}_{\mathrm{3}}\right]}^{\mathrm{T}}={\left[W_{\mathrm{c}\mathrm{p}}{P}_{\mathrm{s}\mathrm{m}}{P}_{\mathrm{r}\mathrm{m}}\right]}^{\mathrm{T}}$；$w=I_{\mathrm{s}\mathrm{t}}$。通过在稳定工作状态点进行泰勒一阶展开从而将模型进行线性化，获得系统的线性定常系数矩阵如下：

$$\begin{array}{c}\frac{\mathrm{d}({\mathit{x}}_{\mathrm{0}}+\mathrm{\delta }\mathit{x})}{\mathrm{d}t}={\dot{\mathit{x}}}_{\mathrm{0}}+\mathrm{\delta }\dot{\mathit{x}}=\\ f({\mathit{x}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}})+\frac{\partial f}{\partial \mathit{x}}{|}_{(x_{\mathrm{0}},u_{\mathrm{0}})}\mathrm{\delta }\mathit{x}+\frac{\partial f}{\partial \mathit{u}}{|}_{\left(x_{\mathrm{0}},u_{\mathrm{0}}\right)}\mathrm{\delta }\mathit{u}\end{array}$$

（3）

式中：${\mathit{x}}_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}}$分别是稳定工作状态点的各状态值和此时系统的输入值； $\mathrm{\delta }\mathit{x}\mathrm{、}\mathrm{\delta }\mathit{u}$为系统标称轨迹在平衡点邻域内的微小增量，在燃料电池系统中，还需要确定平衡点的电流扰动，即$w_{\mathrm{0}}$。通过离线仿真得到系统在扰动电流$w_{\mathrm{0}}=\mathrm{140}\mathrm{A}$时的各状态值，${\mathit{x}}_{\mathrm{0}}=[\mathrm{0.001}\mathrm{914},\mathrm{0.012}\mathrm{6},\mathrm{1.444}\times {\mathrm{10}}^{\mathrm{5}},\mathrm{0.033}\mathrm{63},\mathrm{1.743}\times {\mathrm{10}}^{\mathrm{5}}{,\mathrm{6}\mathrm{684}]}^{\mathrm{T}}$，此时$u_{\mathrm{0}}=\mathrm{122.11}\mathrm{V}$，从而得到系统在平衡点处的六阶线性状态空间方程：

$$\left\{\begin{array}{c}\mathrm{\delta }\dot{\mathit{x}}=A\mathrm{\delta }\mathit{x}+B\mathrm{\delta }\mathit{u}+D_{\mathrm{1}}\mathrm{\delta }\mathit{w}\\ \mathrm{\delta }\mathit{y}=C_{\mathrm{1}}\mathrm{\delta }\mathit{x}\\ \mathrm{\delta }\mathit{z}=C_{\mathrm{2}}\mathrm{\delta }\mathit{x}+D_{\mathrm{2}}\mathrm{\delta }\mathit{w}\end{array}\right.$$

（4）

2.2　LQR控制设计

LQR最优控制的本质是一种满足控制要求的状态反馈控制，该控制易于实现，且具有较强的工程性。文献［4］表明，LQR控制的特点在于线性最优控制的结果，可以用于非线性程度低的小型燃料电池系统，其计算和实现比一般的非线性控制方法更加容易。

此外，LQR最优控制除满足性能指标最小外，还具有优秀的频响特性和稳定裕度，即实现最优的极点配置，同时将其与误差积分控制相结合，可以增强系统的鲁棒性，抗扰动能力增强。

LQR最优控制的设计目标在于以最小的输入代价，控制过氧比在电流扰动变化的情况下，快速稳定的维持在${\lambda }_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{2}$，即寻找最优控制输入，使得如下代价方程取最小值：

$$J={\displaystyle \underset{\mathrm{0}}{\overset{\infty }{\int }}}\delta {\mathit{z}}^T{Q}_Z\delta \mathit{z}+{\mathit{q}}^{\mathrm{T}}{Q}_{\mathrm{I}}\mathit{q}+\mathrm{\delta }{\mathit{u}}^{\mathrm{T}}R\mathrm{\delta }\mathit{u}\mathrm{d}t$$

（5）

其中$\dot{\mathit{q}}={\mathit{W}}_{\mathrm{c}\mathrm{p}}^{*}-{\mathit{W}}_{\mathrm{c}\mathrm{p}}$是为了增加系统克服扰动能力而引入的系统扩张状态，其目的是在存在扰动的情况下，仍然希望系统输出能够很好的跟随参考输入。选择空压机出口流量作为扩张状态主要是由于其便于观测的特点，同时出口质量流量也是直接影响阴极过氧比的上游参数，通过对其误差的积分来消除静态误差，增强了系统的鲁棒性，也增加了整个系统的动态特性。LQR控制器状态反馈矩阵由如下ARE（黎卡提）方程解得：

$$\begin{array}{l}\mathit{P}\widehat{\mathit{A}}+{\widehat{\mathit{A}}}^{\mathrm{T}}\mathit{P}+{\widehat{\mathit{Q}}}_x-\mathit{P}\widehat{\mathit{B}}{\mathit{R}}^{-\mathrm{1}}{\widehat{\mathit{B}}}^{\mathrm{T}}\mathit{P}=\mathrm{0}\\ K={\mathit{R}}^{-\mathrm{1}}{\widehat{\mathit{B}}}^{\mathrm{T}}\mathit{P}\end{array}$$

（6）

式中： $\widehat{\mathit{A}}=[A,\mathrm{0};-C_{\mathrm{1}},\mathrm{0}]$，$\widehat{\mathit{B}}=[B;\mathrm{0}]$，${\widehat{\mathit{Q}}}_x=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(Q_x,Q_{\mathrm{I}})$，$Q_x={\mathit{C}}_{\mathrm{2}}^{\mathrm{T}}{Q}_Z{\mathit{C}}_{\mathrm{2}}$。最终，使得式（5）取得最小值的最优控制输入有如下形式：

$$\mathrm{\delta }\mathit{u}=-K{[\mathrm{\delta }\widehat{\mathit{x}}\mathit{q}]}^{\mathrm{T}}=-K_p\mathrm{\delta }\widehat{\mathit{x}}-K_{\mathrm{I}}\mathit{q}$$

（7）

其中$\delta \widehat{\mathit{x}}$是由降维状态观测器对状态估计得到的。由式（4）可知，阴极过氧比$\mathrm{\delta }\mathit{z}$除受系统状态$\mathrm{\delta }\mathit{x}$影响外，还受到电流扰动$\mathrm{\delta }\mathit{w}$的影响；然而在考虑代价函数时，不得不将电流扰动影响忽略，构成标准的仅包含状态变量和控制输入的线性二次型指标，故在得到LQR的状态反馈增益后，根据线性系统对输出线性相加的原理^［9］，对控制输入进行前馈补偿，该补偿采用系统在$I_{\mathrm{s}\mathrm{t}}=\mathrm{140}\mathrm{A}$、${\lambda }_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{2}$下的线性化后的前馈控制器设计结果。进一步可以通过神经网络对系统在多个平衡点处的数据进行训练和拟合，找到一个适用于全局的从电流$I_{\mathrm{s}\mathrm{t}}$到控制输入$U_{\mathrm{c}\mathrm{m}}$的映射。

2.3　降维状态观测器设计

为了实现状态反馈控制，必须得到每一时刻系统内部所有状态的值，虽然在仿真过程中可以直接得到这些数值，然而在实际工程实现中，除输出变量${\left[W_{\mathrm{c}\mathrm{p}}{P}_{\mathrm{s}\mathrm{m}}{P}_{\mathrm{r}\mathrm{m}}\right]}^{\mathrm{T}}$可由传感器直接观测得到外，其余3个状态变量均需通过输入和输出信息构造状态观测器得到。针对已经在$I_{\mathrm{s}\mathrm{t}}=\mathrm{140}\mathrm{A}$平衡点处线性化后的线性时不变燃料电池空气供给系统，构造降维状态观测器，由$\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}{C}_{\mathrm{1}}=\mathrm{3}$，任取$\mathit{R}$使$\mathit{P}=\left[\begin{array}{c}{C}_{\mathrm{1}}\\ R\end{array}\right]$为非奇异，对状态向量作线性非奇异变换$\left[\begin{array}{c}{\dot{x}}_{\mathrm{1}}\\ {\dot{x}}_{\mathrm{2}}\end{array}\right]=\mathit{P}\mathit{x}$，可证${\dot{\mathit{x}}}_{\mathrm{1}}=\mathit{y}$，于是便可将降维状态观测器的设计简化为对于3维状态向量${\dot{\mathit{x}}}_{\mathrm{2}}$的全维状态观测器设计：

$$\begin{array}{l}\dot{\mathit{z}}=({\dot{\mathit{A}}}_{\mathrm{22}}-\dot{L}{\dot{\mathit{A}}}_{\mathrm{12}})\mathit{z}+[({\dot{\mathit{A}}}_{\mathrm{22}}-\dot{L}{\dot{\mathit{A}}}_{\mathrm{12}})\dot{\mathit{L}}+\\ ({\dot{\mathit{A}}}_{\mathrm{21}}-\dot{L}{\dot{\mathit{A}}}_{\mathrm{11}})]\mathit{y}+({\dot{\mathit{B}}}_{\mathrm{2}}-\dot{L}{\dot{\mathit{B}}}_{\mathrm{1}})\mathit{u}\\ {\widehat{\dot{\mathit{x}}}}_{\mathrm{2}}=\mathit{z}+\dot{L}\mathit{y}\end{array}$$

（8）

式中： ${\dot{A}}_{\mathrm{11}},{\dot{A}}_{\mathrm{12}},{\dot{A}}_{\mathrm{21}},{\dot{A}}_{\mathrm{22}},{\dot{B}}_{\mathrm{1}},{\dot{B}}_{\mathrm{2}}$是基于变换$\dot{\mathit{x}}=\mathit{P}\mathit{x}$得到的状态方程$\mathrm{\delta }\dot{\dot{\mathit{x}}}=\dot{A}\mathrm{\delta }\dot{\mathit{x}}+\dot{B}\mathrm{\delta }\mathit{u}$的系数矩阵并根据相应维数划分的分块矩阵，最后再通过观测到的变换状态$\dot{\mathit{x}}$可以重构状态$\mathit{x}$：

$$\widehat{\mathit{x}}=Q_{\mathrm{1}}\mathit{y}+Q_{\mathrm{2}}\left(\mathit{z}+\dot{\mathit{L}}\mathit{y}\right)$$

（9）

其中$\mathit{Q}={\mathit{P}}^{-\mathrm{1}}=\left[Q_{\mathrm{1}}|Q_{\mathrm{2}}\right]$。基于单平衡点线性化设计的LQR控制器与燃料电池系统的控制过程如图1所示，其中${\mathit{W}}_{\mathrm{c}\mathrm{p}}^{*}$可根据当前扰动电流计算得到。

图1  单平衡点LQR控制框图

Fig.1  Single balance point LQR control block diagram

3 多平衡点LQR控制器

通过对单平衡点线性化后的模型进行LQR控制器设计，在平衡点附近小范围内出现工况变化时，可以有效的控制过氧比的动态响应，然而在偏离平衡点较远处，过氧比可能会由于模型失配的原因产生静态误差，并且其动态响应指标也会比平衡点位置稍差。虽然可以采用神经网络训练多平衡点处的数据来消除过氧比与期望值之间的静态误差，但采用的离线数据过多，在实际工程实现上不够简便，并且其动态性能仍然会受到影响。

为了让LQR控制器在全局工况下都有良好的控制效果，引入多模型控制理论的思想，利用非线性过程在多个工作点上的线性化模型，将系统模型依照状态参数分成一系列子空间，在每个子空间内找到一个平衡点对该子空间内的非线性模型线性化，同时基于该线性子模型设计独立的LQR控制器。文献［4］表明，在一个小工况范围内，数学模型的失配度较小，采用LQR控制器可以较好的满足小范围内过氧比的快速响应性能，这也印证了通过多模型方法，独立设计各自LQR控制器的可行性。关键点在于设计一个有效的控制器合成方案，将所得到的所有子空间的线性控制器组合起来，最终实现全工况范围内的有效控制。

T-S模糊多模型方法是Takagi-Sugeno^［10］提出的一种多模型组织方法，本文采用这样的方法将系统划分成多个相互重叠的模糊集合，各模糊集合上建立局部线性模型，通过一组“如果-则”模糊规则来描述非线性系统，整个非线性系统最终就是各个子空间数学模型与相应隶属度加权的线性组合。T-S模糊模型可描述如下：

$$\begin{array}{l}{R}^i:\mathrm{I}\mathrm{F}{x}_{\mathrm{1}}\mathrm{i}\mathrm{s}{M}_{\mathrm{1}}\mathrm{a}\mathrm{n}\mathrm{d}...\mathrm{a}\mathrm{n}\mathrm{d}{x}_k\mathrm{i}\mathrm{s}{M}_k\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}\dot{\mathit{x}}=\\ {\mathit{A}}_i\mathit{x}+{\mathit{B}}_i\mathit{u}i=\mathrm{1,2},...,l\end{array}$$

（10）

式中： $M_i(i=\mathrm{1,2},\dots ,k)$是模糊集合，$l$为模糊规则个数，${\mathit{A}}_i\mathrm{、}{\mathit{B}}_i$是第i个子系统的系数矩阵，最终可得整个系统的状态方程：

$$\dot{\mathit{x}}={\displaystyle \sum _{i=\mathrm{1}}^l}{h}_i(\mathit{x})[{\mathit{A}}_i\mathit{x}+{\mathit{B}}_i\mathit{u}]$$

（11）

其中$h_i(x)={\mu }^i(\mathit{x})/{\displaystyle \sum _{j=\mathrm{1}}^l}{\mu }^j(\mathit{x});{\displaystyle \sum _{i=\mathrm{1}}^l}{h}_i(\mathit{x})=\mathrm{1}$，${\mu }^i(\mathit{x})$表示x属于$M_i$的隶属度。本文以干扰电流为指标，当干扰电流分别为100、140、180 A时，对系统进行线性化，并以所得线性化模型计算LQR控制器的状态反馈增益，得到3个局部模型及相应的线性控制器，其T-S模糊多模型可描述如下：

$$\begin{array}{l}{R}^{\mathrm{1}}:\mathrm{I}\mathrm{F}{I}_{\mathrm{s}\mathrm{t}}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{100}\mathrm{A}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}\mathit{u}=-K_{\mathrm{1}}\mathit{x}\\ R^{\mathrm{2}}:\mathrm{I}\mathrm{F}{I}_{\mathrm{s}\mathrm{t}}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{140}\mathrm{A}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}\mathit{u}=-K_{\mathrm{2}}\mathit{x}\\ R^{\mathrm{3}}:\mathrm{I}\mathrm{F}{I}_{\mathrm{s}\mathrm{t}}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{180}\mathrm{A}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{n}\mathit{u}=-K_{\mathrm{3}}\mathit{x}\end{array}$$

（12）

其中$K_{\mathrm{1}},K_{\mathrm{2}},K_{\mathrm{3}}$分别为3个局部模型的状态反馈增益。在Matlab/Simulink平台搭建模糊多模型控制系统结构如图2所示，其中任一时刻的实际控制输入是3个LQR控制器输出的加权和，即$\mathit{u}=-{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{3}}}{h}_i{K}_i\mathit{x}$，规则式（12）中$I_{\mathrm{s}\mathrm{t}}$对应的隶属度函数如图3所示。

图2  Simulink搭建模糊多模型控制系统结构图

Fig.2  Fuzzy multi-model control system structure by Simulink

图3  扰动电流的隶属度函数

Fig.3  Membership function corresponding to disturbance current

4 数值仿真

数值仿真在Matlab/Simulink平台进行，仿真中将前馈分别结合单平衡点LQR、多平衡点LQR及PID控制器进行对比。在单平衡点LQR方案中选择扰动电流$I_{\mathrm{s}\mathrm{t}}=\mathrm{140}\mathrm{A}$、过氧比${\lambda }_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{2}$时进行线性化并设计相应的LQR控制器，在多平衡点LQR方案中，将非线性系统依据扰动电流的数值分为3个模糊集合，并采用上一节所展示的隶属度函数对3个LQR控制器的输出权重进行在线计算，最后加权得到实时的控制输入。系统所施加的扰动电流从80 A阶跃上升到180 A，然后再阶跃下降到100 A，如图4所示。

图4  燃料电池电堆输出电流

Fig.4  Output current of PEMFC

使用神经网络前馈分别与单平衡点LQR和PID控制相结合的两种复合控制策略的效果如图5所示，其中单平衡点LQR策略选择扰动电流为140 A时的系统工作点对系统线性化，并基于线性模型设计LQR控制器。由图5可知，在扰动电流在15 s处从120 A阶跃到140 A时，过氧比都存在下冲的现象，原因在于此刻拉载电流突然上升，导致氧气消耗量突然增加；为维持过氧比不变，控制输入传递给空压机更大的电压输入，使其转速增加；由于系统的迟滞和容腔效应，短时间供气没有及时供上，导致过氧比出现下冲，下冲后迅速回升，可见LQR策略回升极快，有少量超调随后过氧比在0.5 s时间内达到稳定值2，而PID没有超调，呈现过阻尼现象，需要3 s才能达到平衡状态。在变载情况下，平衡点附近处，前馈加单平衡点LQR控制过氧比的效果要优于前馈加PID控制，但是在远离平衡点位置LQR策略会由于模型失配产生静态误差，但其动态响应效果仍然优于传统PID控制策略。

图5  单平衡点LQR控制效果图

Fig.5  Single balance point LQR control effect

采用基于T-S模糊多模型方法的多平衡点LQR策略与单平衡点LQR策略控制过氧比的效果对比如图6所示，多平衡点策略分别选取100、140、180 A时的工况作为平衡点，所以在15 s和35 s处，扰动电流都阶跃变化到140 A水平，此时多平衡点策略与单平衡点策略的曲线重合。在非平衡点处，如10 s和30 s处，此时扰动电流分别为120 A和 160 A，如图6中局部放大图所示。

图6  多平衡点LQR控制效果图

Fig.6  Multi-balance point LQR control effect

由图6可见，在10 s处，单平衡点LQR控制下过氧比的稳态误差为0.014，多平衡点LQR策略控制下过氧比稳态误差为0.007，稳态误差降低了50%，此时过氧比的超调量0.017，而单平衡点策略下过氧比的超调量为0.037，超调量降低了46%；在30 s处，单平衡点LQR控制下过氧比的稳态误差为0.013，多平衡点LQR策略控制下过氧比稳态误差为0.006，稳态误差降低了46%，此时过氧比的超调量0.006，而单平衡点策略下过氧比的超调量为0.012，超调量降低了50%。

仿真结果验证多平衡点LQR在非平衡点处对于过氧比的控制有更好的动态响应，同时明显降低了静态误差，能够将过氧比更好地稳定在${\lambda }_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{2}$。多平衡点LQR控制策略在全局范围内都展现了更佳的控制效果。

5 结论

本文建立了六阶空气供给系统控制模型。基于T-S模糊多模型方法，将系统依据扰动电流分为3个模糊子空间，对每个子空间的LQR控制器的输出加权求和，最终得到实时的系统输入，实现了对空气供给系统的最优过氧比控制。

将其与前馈加PID和前馈加单平衡点LQR控制策略对比，仿真结果表明：神经网络前馈加PID控制策略下，过氧比的平均稳定时长为3.597 s，过渡过程时间过长，不适合在复杂变工况下作为过氧比的控制策略；单平衡点LQR控制策略下，过氧比的平均稳态误差为0.026，平均稳定时长为1.1 s，而多平衡点LQR控制策略下，过氧比平均稳态误差为0.007，平均稳定时长为1.3 s，所以基于T-S模糊多模型理论的多平衡点LQR控制器不仅能在平衡点处快速有效的将过氧比稳定在期望值，而且在偏离平衡点位置，仍然展现了稳定且快速的控制能力，具有较短的整定时长，稳态误差在1%之内，因此可以支持车用燃料电池在复杂的工况下提供准确且快速的空气供应。