CN100483276C - 一种基于噪声估计的自适应状态反馈预测控制方法 - Google Patents

Abstract

本发明公开了一种基于噪声估计的自适应状态反馈预测控制方法，包括以下步骤：A、利用Singular Pencil Model模型表示实际被控过程；B、采用递推的扩展卡尔曼滤波参数估计方法估计所述利用Singular PencilModel模型表示的实际被控过程的状态和模型参数；C、利用估计出来的状态和模型参数计算当前的控制作用，并将计算出来的控制作用加到实际被控过程中；D、在下一个执行周期，重复执行步骤B和C，实现对实际被控过程的自适应状态反馈预测控制。利用本发明，充分考虑了被控过程中含有的噪声信息和状态信息，大大提高了控制方法抑制噪声的能力。另外，本发明提供的控制方法计算量小，便于工业应用。

Description

一种基于噪声估计的自适应状态反馈预测控制方法

技术领域

本发明涉及自动化控制技术领域，尤其涉及一种基于噪声估计的自适应状态反馈预测控制方法。

背景技术

在工业自动化控制领域，应用预测控制技术对被控过程进行控制具有广泛的应用。

现有的应用预测控制技术多是基于离线辨识的模型，或是根据物理和化学反应机理推导出来的模型，例如模型算法控制(Model AlgorithmControl)、动态矩阵控制(Dynamic Matrix Control)和广义预测控制(GeneralPredictive Control)等。另外，现有的应用预测控制技术也有将在线辨识和上述预测控制相结合的自适应预测控制方法，例如自适应的广义预测控制，模型参考自适应预测控制方法等。

在上述各种现有的应用预测控制技术中，由于噪声信息难以确定，尤其是工业上被控过程中含有不可测的噪声信息，导致现有的控制方法无法充分考虑被控过程中含有的噪声信息。即使个别现有的控制方法考虑了被控过程中含有的噪声信息，考虑的也大都是可测噪声信息，而在工业上被控过程中含有的噪声信息多是不可测噪声信息，所以现有的控制方法并没有充分考虑被控过程中含有的不可测噪声信息。

另外，由于大部分应用预测控制技术都没有引入状态反馈，所以现有的控制方法都没有充分考虑被控过程中含有的状态信息。即使个别现有的方法虽引入了状态反馈，例如状态反馈预测控制，但是受到状态是否可测的限制，如果状态是可测的，则可以直接利用测出的状态信息进行状态反馈预测控制；如果状态是不可测的，则还要单独设计状态观测器，利用设计的状态观测器对状态信息进行观测，然后再利用测出的状态信息进行状态反馈预测控制，实现过程非常复杂，实现成本也很高。

因此，基于上述两个原因，即现有的控制方法没有充分考虑被控过程中含有的噪声信息和状态信息，导致了现有的控制方法抑制噪声的能力很差，这也是目前应用预测控制技术急需解决的重要问题。

发明内容

(一)要解决的技术问题

针对上述现有技术存在的问题，本发明的主要目的在于提供一种基于噪声估计的自适应状态反馈预测控制方法，以充分考虑被控过程中含有的噪声信息和状态信息，提高控制方法抑制噪声的能力。

(二)技术方案

为达到上述目的，本发明的技术方案是这样实现的：

一种基于噪声估计的自适应状态反馈预测控制方法，该方法包括：

a、利用Singular Pencil Model模型表示实际被控过程；

假设x_k是状态变量，u_k是输入变量，y_k是输出变量，e_k为噪声变量，E_*，A_*，B_*，C_*，E₀为矩阵，I为单位阵，z为前向移位算子，则建立公式

表示实际被控过程；

假设

$w_k={[y_k^T,u_k^T]}^T,$ ${\tilde{G}}_{*}\left(w_k\right)r=[-A_{*}|B_{*}]w_k,$ ${\tilde{G}}_0\left(w_k\right)r=[-A_0|B_0]w_k+{\mathrm{Jw}}_k,$ $C_{*}{e}_k={\tilde{C}}_{*}\left(e_k\right){\mathrm{\η}}_k,$ 其中矩阵

和

根据r中参数的特定位置构建，矩阵

根据η中参数的特定位置构建，则将用于表示实际被控过程的公式

变形为公式 $\left\{\begin{array}{c}{x}_{k+1}=E_{*}{x}_k+{\tilde{G}}_{*}\left(w_k\right)r+C_{*}{e}_k\\ 0=E_0{x}_k+{\tilde{G}}_0\left(w_k\right)r-{\mathrm{Jw}}_k+e_k\end{array}\right.;$

假设 ${\tilde{F}}_k=\left[\begin{array}{ccc}{E}_{*}& {\tilde{G}}_{*}\left(w_k\right)& {\tilde{C}}_{*}\left(e_k\right)\\ 0& I& 0\\ 0& 0& I\end{array}\right],$ $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T,$ ${\tilde{H}}_k=\left[\begin{array}{ccc}{E}_0& {\tilde{G}}_0\left(w_k\right)& 0\end{array}\right],$ 则将公式 $\left\{\begin{array}{c}{x}_{k+1}=E_{*}{x}_k+{\tilde{G}}_{*}\left(w_k\right)r+C_{*}{e}_k\\ 0=E_0{x}_k+{\tilde{G}}_0\left(w_k\right)r-{\mathrm{Jw}}_k+e_k\end{array}\right.$ 变形为公式 $\left\{\begin{array}{c}{s}_{k+1}={\tilde{F}}_k{s}_k\\ y_k={\tilde{H}}_k{s}_k+e_k\end{array}\right.;$

其中，I是单位矩阵；w_k是由输入输出数据组成的矩阵；r是待估计的***模型参数向量；η是待估计的噪声模型参数向量；

b、采用递推的扩展卡尔曼滤波参数估计方法估计所述利用SingularPencil Model模型表示的实际被控过程的状态和模型参数；

所述递推的扩展卡尔曼滤波参数估计算法公式为

$\left\{\begin{array}{c}{K}_k={\mathrm{\Ξ}}_k{P}_k{\tilde{H}}_k^T{({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)}^{-1}\\ P_{k+1}={\mathrm{\Ξ}}_k{P}_k{\mathrm{\Ξ}}_k^T-K_k({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)K_k^T\\ {\hat{s}}_{k+1}={\tilde{F}}_k{\hat{s}}_k+K_k(y_k-{\tilde{H}}_k{\hat{s}}_k)\end{array}\right.,$

式中， ${\hat{R}}_{k+1}={\hat{R}}_k+\frac{1}{k+1}({\hat{e}}_k{\hat{e}}_k^T-{\hat{R}}_k),$ ${\mathrm{\Ξ}}_k=\frac{\∂\left({\tilde{F}}_k{s}_k\right)}{\∂s_k}{|}_{s_k={\hat{s}}_k},$

$\frac{\∂\left({\tilde{F}}_k{s}_k\right)}{\∂s_k}=\left[\begin{array}{ccc}{E}_{*}-C_{*}{E}_0& {\tilde{G}}_{*}\left(w_k\right)-C_{*}{\tilde{G}}_0\left(w_k\right)& {\tilde{C}}_{*}\left(e_k\right)\\ 0& I& 0\\ 0& 0& I\end{array}\right]$

公式 $P_{k+1}={\mathrm{\Ξ}}_k{P}_k{\mathrm{\Ξ}}_k^T-K_k({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)K_k^T$ 用于修正K_k，将修正后的K_k代入公式 ${\hat{s}}_{k+1}={\tilde{F}}_k{\hat{s}}_k+K_k(y_k-{\tilde{H}}_k{\hat{s}}_k),$ 计算得到包含所述实际被控过程状态和模型参数的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T;$

c、利用估计出来的状态和模型参数计算当前的控制作用，并将计算出来的控制作用加到实际被控过程中；其中所述利用估计出来的状态和模型参数计算当前的控制作用，具体包括：

c1、利用包含实际被控过程状态和模型参数的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T$ 构建矩阵A、矩阵B、矩阵V和状态反馈向量x；

c2、利用构建的矩阵A、B和V构建状态反馈预测控制方法公式

$\left\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Ve}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+e\left(k\right)\end{array}\right.;$

c3、根据参考轨迹y^s(k+p)＝αy(k)+(1-α)ω、性能指标 $J_p=[\hat{y}(k+p)-y^s(k+p)]Q[\hat{y}(k+p)-y^s(k+p)]+\mathrm{\Δ}{u}^T\left(k\right)\mathrm{\λ\Δu}\left(k\right)$ 和预测公式 $\left\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Ve}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+e\left(k\right)\end{array}\right.,$ 计算增量形式的最优控制率；

Δu(k)＝(S^TQS+λ)^-1S^TQ(y^s(k+p)-y(k)-Δ_KX-Δ_SU-Δ_TY)，

其中，ω是输入设定值，α为滤波系数，

${\mathrm{\Δ}}_{\mathrm{KX}}=\left[\begin{array}{c}(C_1{A}^{p_1}-C_1{A}^{p_1}{\mathrm{VC}}_1)(x\left(k\right)-x(k-p_1))\\ (C_2{A}^{p_2}-C_2{A}^{p_2}{\mathrm{VC}}_2)(x\left(k\right)-x(k-p_2))\\ \·\\ \·\\ \·\\ (C_r{A}^{p_r}-C_r{A}^{p_r}{\mathrm{VC}}_r)(x\left(k\right)-x(k-p_r))\end{array}\right],$ ${\mathrm{\Δ}}_{\mathrm{SU}}=\left[\begin{array}{c}\underset{i=1}{\overset{p_1}{\mathrm{\Σ}}}(S_1\left(p_1\right)-S_1\left(i\right))\mathrm{\Δu}(k-i)\\ \underset{i=1}{\overset{p_2}{\mathrm{\Σ}}}(S_2\left(p_2\right)-S_2\left(i\right))\mathrm{\Δu}(k-i)\\ \·\\ \·\\ \·\\ \underset{i=1}{\overset{p_r}{\mathrm{\Σ}}}(S_r\left(p_r\right)-S_r\left(i\right))\mathrm{\Δu}(k-i)\end{array}\right],$

${\mathrm{\Δ}}_{\mathrm{TY}}=\left[\begin{array}{c}{C}_1{A}^{p_j-1}V(y\left(k\right)-y(k-p_1))\\ C_2{A}^{p_2-1}V(y\left(k\right)-y(k-p_2))\\ \·\\ \·\\ \·\\ C_r{A}^{p_r-1}V(y\left(k\right)-y(k-p_r))\end{array}\right],$ $S\left(p\right)=\left[\begin{array}{c}{S}_1\left(p_1\right)\\ S_2\left(p_2\right)\\ \·\\ \·\\ \·\\ S_r\left(p_r\right)\end{array}\right],$ $y^s=\left[\begin{array}{c}{y}_1^s(k+p_1)\\ y_2^s(k+p_2)\\ \·\\ \·\\ \·\\ y_r^s(k+p_r)\end{array}\right],$

$S_j\left(p_j\right)=\underset{i=1}{\overset{p_j}{\mathrm{\Σ}}}{C}_j{A}^{i-1}B,$ Q为输出预测误差加权阵，λ为控制量的加权系数，Δu(k)为控制作用u(k)的增量；

c4、将Δu(k)＝(S^TQS+λ)^-1S^TQ(y^s(k+p)-y(k)-Δ_KX-Δ_SU-Δ_TY)代入公式u(k)＝u(k-1)+Δu(k)，计算得到当前的控制作用u(k)；

d、在下一个执行周期，重复执行步骤b和c。

所述被控过程用二阶单入单出模型表示，r^T＝[a₂ a₁ b₂ b₁]，η^T＝[c₂ c₁]，x^T＝[x₂ x₁]，所述步骤c1包括：用r^T中的a_i构造矩阵 $A=\left[\begin{array}{cc}0& -a_2\\ 1& -a_1\end{array}\right],$ 用r^T中的b_i构造矩阵 $B=\left[\begin{array}{c}{b}_2\\ b_1\end{array}\right],$ 用η^T中的c_i构造矩阵 $V=\left[\begin{array}{c}{c}_2\\ c_1\end{array}\right],$ 用状态x^T作状态反馈。

步骤c中所述将计算出来的控制作用加到实际被控过程中包括：将计算出来的控制作用通过数模转化，由数字量转化成模拟量，再送到执行元件中。

(三)有益效果

从上述技术方案可以看出，本发明具有以下有益效果：

1、利用本发明，通过将扩展卡尔曼滤波参数估计方法和状态反馈预测控制方法相结合，可实现同时估计被控过程中含有的状态信息和模型参数，并且估计出来的状态信息可直接用于状态反馈，估计出来的模型参数中的噪声信息可直接用于控制器设计，相当于在控制作用中引入了噪声的校正项，充分考虑了被控过程中含有的噪声信息和状态信息，使该控制方法和其它控制方法相比具有较强的抑制噪声能力，大大提高了控制方法抑制噪声的能力。

2、利用本发明，由于采用的是单值状态反馈预测控制算法，所以本发明提供的控制方法计算量小，便于工业应用。

3、利用本发明，由于采用了SPM(Singular Pencil Model)模型将状态信息和参数信息进行分离，因此可以非常方便地对状态信息和参数信息同时进行估计。

4、利用本发明，由于可以直接利用在线估计出来的状态信息做状态反馈，所以无需再设计状态观测器，使控制器的设计更简单，简化了实现过程，降低了实现成本。

5、利用本发明，由于在每个执行周期内都在线估计被控过程的状态信息和参数信息，动态更新控制器的参数，所以控制器能适应被控过程的变化。

附图说明

图1为本发明提供的自适应状态反馈预测控制总体技术方案的实现流程图；

图2为本发明提供的基于噪声估计的自适应状态反馈预测控制结构示意图；

图3为现有技术中广义预测控制(GPC)方法和本发明提供的控制方法对阶跃扰动的抑制能力仿真对比示意图；

图4为本发明实施例中燃料电池流量***控制Simulink模块示意图；

图5为本发明实施例中燃料电池流量***方法流程图；

图6为本发明实施例中燃料电池流量***的实际控制效果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明的核心思想是：通过将扩展卡尔曼滤波参数估计方法和状态反馈预测控制方法相结合，可实现同时估计被控过程中含有的状态信息和模型参数，并且估计出来的状态信息可直接用于状态反馈，估计出来的模型参数中的噪声信息可直接用于控制器设计，相当于在控制作用中引入了噪声的校正项，充分考虑了被控过程中含有的噪声信息和状态信息，使该控制方法和其它控制方法相比具有较强的抑制噪声能力，大大提高了控制方法抑制噪声的能力。

上述扩展卡尔曼滤波参数估计方法和状态反馈预测控制方法都是目前应用预测控制技术领域广泛应用的控制方法，下面对其原理分别进行详细阐述。

卡尔曼滤波器是由美国学者卡尔曼(Rudolph Kalman)开发设计，最初是专为航天器导航而研发的，后逐渐应用到许多领域，包括各种核电站设备、人口统计建模、制造业、地层放射性探测以及模糊逻辑和神经网络学习等。卡尔曼滤波器主要用来预估那些只能被***本身间接或不精确观测的***状态。普通卡尔曼滤波器用于估计线性***的状态，对于非线性***的状态，则应使用扩展卡尔曼滤波器进行估计。事实上，在工程中几乎所有的进程都是非线性的。一些非线性***非常接近线性***，而另一些则差别甚大。因而，扩展卡尔曼滤波器更有应用价值。扩展卡尔曼滤波器只是简单地将线性卡尔曼滤波原理推广到非线性***。由于本发明方法中涉及噪声信息的估计，而噪声信息需要用输入输出信息来估计，这是一个非线性问题，因而应用扩展卡尔曼滤波进行参数估计。

状态反馈预测控制方法是将状态空间模型用于预测控制中的一种算法，由于引入了状态信息，而状态中又包含了不可测干扰的影响，相当于引入了未知干扰的前馈，大大提高了控制方法的抗干扰效果。单值的状态反馈预测控制算法还具有计算量小的优点，且和多值的算法控制效果基本相当，便于工业应用。但由于需要状态信息，往往由于状态的不可测而要另外构造状态观测器。本发明方法在参数估计出来的同时就可获得状态信息，使状态反馈预测控制方法更有工业应用价值。

本发明在将扩展卡尔曼滤波参数估计方法和状态反馈预测控制方法结合在一起的过程中，需要在状态反馈预测控制方法中考虑噪声信息，推导出带有噪声信息的状态反馈预测控制方法公式，并将在线估计出来的状态信息和参数信息用于控制方法中，形成自适应的状态反馈预测控制方法。

如图1所示，图1为本发明提供的自适应状态反馈预测控制总体技术方案的实现流程图，该方法包括以下步骤：

步骤101：利用Singular Pencil Model模型表示实际被控过程；

步骤102：采用递推的扩展卡尔曼滤波参数估计方法估计所述利用Singular Pencil Model模型表示的实际被控过程的状态和模型参数；

步骤103：利用估计出来的状态和模型参数计算当前的控制作用，并将计算出来的控制作用加到实际被控过程中；

步骤104：在下一个执行周期，重复执行步骤102和步骤103，实现对实际被控过程的自适应状态反馈预测控制。

上述步骤101中所述Singular Pencil Model模型是一种特殊的模型结构，可以用来表示线性动态***，其模型结构如公式(1)所示：

$\left[\begin{array}{cc}E-\mathrm{FD}& G\end{array}\right]\left[\begin{array}{c}x\\ w\end{array}\right]=0---\left(1\right)$

其中，是状态变量，也叫内部变量或辅助变量。

是外部变量，包括***的输入、输出变量。E、F是(n+p)×n矩阵，G是(n+p)×(p+m)矩阵。D是线性算子。对于连续***，D可以是Laplace算子s，对于离散***，D可以是前向移位算子z。

以状态空间模型表示的多入多出***可以用Singular Pencil Model模型表示为公式(2)：

$\left[\begin{array}{ccc}A-\mathrm{Iz}& 0& B\\ C& -I& D\end{array}\right]\left[\begin{array}{c}{x}_k\\ y_k\\ u_k\end{array}\right]=0---\left(2\right)$

其中，u_k是输入变量，y_k是输出变量，(A，B，C，D)是状态空间矩阵。滑动自回归平均模型 $\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{A}_i{z}^{n-i}{y}_k=\underset{i=0}{\overset{m}{\mathrm{\Σ}}}{B}_i{z}^{m-i}{u}_k$ 也可以用Singular PencilModel模型表示为公式(3)：

其中，A_i和B_i是矩阵，z是前向移位算子。

由于Singular Pencil Model模型具有如下的特征，因而得到了广泛应用。

1)、统一：Singular Pencil Model模型可包含线性控制***中最常用的状态空间模型和输入输出模型。

2)、简单：可通过对公式(1)中矩阵E，F，G做一系列简单的初等行和列变换完成线性变换。

3)、通用：和状态空间描述不同的是，Singular Pencil Model模型还可以描述非正***，如PID控制器、严格正***的逆***等等。

利用Singular Pencil Model模型，应用扩展卡尔曼滤波方法可同时估计出被控过程的状态信息和参数信息。

上述步骤101中所述利用Singular Pencil Model模型表示实际被控过程包括：假设x_k是状态变量，u_k是输入变量，y_k是输出变量，e_k为噪声变量，E_*，A_*，B_*，C_*，E₀为矩阵，I为单位阵，z为前向移位算子，则建立公式(4)表示实际被控过程，所述公式(4)如下：

上述在建立公式(4)表示实际被控过程之后进一步包括：假设 $w_k={[y_k^T,u_k^T]}^T,$ ${\tilde{G}}_{*}\left(w_k\right)r=[-A_{*}|B_{*}]w_k,$ ${\tilde{G}}_0\left(w_k\right)r=[-A_0|B_0]w_k+{\mathrm{Jw}}_k,$ $C_{*}{e}_k={\tilde{C}}_{*}\left(e_k\right){\mathrm{\η}}_k,$ 其中矩阵

和

根据r中参数的特定位置构建，矩阵

根据η中参数的特定位置构建，则将用于表示实际被控过程的公式(4)变形为公式(5)，所述公式(5)如下：

$\left\{\begin{array}{c}{x}_{k+1}=E_{*}{x}_k+{\tilde{G}}_{*}\left(w_k\right)r+C_{*}{e}_k\\ 0=E_0{x}_k+{\tilde{G}}_0\left(w_k\right)r-{\mathrm{Jw}}_k+e_k\end{array}\right.---\left(5\right)$

假设 ${\tilde{F}}_k=\left[\begin{array}{ccc}{E}_{*}& {\tilde{G}}_{*}\left(w_k\right)& {\tilde{C}}_{*}\left(e_k\right)\\ 0& I& 0\\ 0& 0& I\end{array}\right],$ $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T,$ ${\tilde{H}}_k=\left[\begin{array}{ccc}{E}_0& {\tilde{G}}_0\left(w_k\right)& 0\end{array}\right],$ 则将公式(5)变形为公式(6)，所述公式(6)如下：

$\left\{\begin{array}{c}{s}_{k+1}={\tilde{F}}_k{s}_k\\ y_k={\tilde{H}}_k{s}_k+e_k\end{array}\right.---\left(6\right)$

由于e_k不能直接测量，可用公式(7)来估计e_k，所述公式(7)如下：

${\hat{e}}_k=y_k-{\hat{H}}_k{\hat{s}}_k---\left(7\right)$

将上述公式(7)代入公式(6)可以估计出包含状态信息和参数信息的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T.$ 这样，应用扩展卡尔曼滤波器可以估计s_k，也就是说，在噪声不可测的情况下，状态信息和参数信息可以同时估计出来。

上述步骤102中所述采用递推的扩展卡尔曼滤波参数估计方法估计所述利用Singular Pencil Model模型表示的实际被控过程的状态信息和模型参数信息包括：

所述递推的扩展卡尔曼滤波参数估计算法公式为公式(8)，所述公式(8)如下：

$\left\{\begin{array}{c}{K}_k={\mathrm{\Ξ}}_k{P}_k{\tilde{H}}_k^T{({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)}^{-1}\\ P_{k+1}={\mathrm{\Ξ}}_k{P}_k{\mathrm{\Ξ}}_k^T-K_k({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)K_k^T\\ {\hat{s}}_{k+1}={\tilde{F}}_k{\hat{s}}_k+K_k(y_k-{\tilde{H}}_k{\hat{s}}_k)\end{array}\right.---\left(8\right)$

式中， ${\hat{R}}_{k+1}={\hat{R}}_k+\frac{1}{k+1}({\hat{e}}_k{\hat{e}}_k^T-{\hat{R}}_k),$ ${\mathrm{\Ξ}}_k=\frac{\∂\left({\tilde{F}}_k{s}_k\right)}{\∂s_k}{|}_{s_k={\hat{s}}_k},$

$\frac{\∂\left({\tilde{F}}_k{s}_k\right)}{\∂s_k}=\left[\begin{array}{ccc}{E}_{*}-C_{*}{E}_0& {\tilde{G}}_{*}\left(w_k\right)-C_{*}{\tilde{G}}_0\left(w_k\right)& {\tilde{C}}_{*}\left(e_k\right)\\ 0& I& 0\\ 0& 0& I\end{array}\right]$

公式 $P_{k+1}={\mathrm{\Ξ}}_k{P}_k{\mathrm{\Ξ}}_k^T-K_k({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)K_k^T$ 用于修正K_k，将修正后的K_k代入公式 ${\hat{s}}_{k+1}={\tilde{F}}_k{\hat{s}}_k+K_k(y_k-{\tilde{H}}_k{\hat{s}}_k),$ 计算得到包含所述实际被控过程状态和模型参数的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T.$

另外，步骤102中所述估计实际被控过程的状态信息和模型参数信息中的模型参数信息包括***矩阵参数和噪声参数。

上述步骤103中所述利用估计出来的状态和模型参数计算当前的控制作用包括：

步骤1031：利用包含实际被控过程状态和模型参数的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T$ 构建A、B、V和状态反馈x；

这里，以二阶单入单出模型为例，即假设所述被控过程用二阶单入单出模型表示，在二阶单入单出模型中，r^T＝[a₂ a₁ b₂ b₁]，η^T＝[c₂ c₁]，x^T＝[x₂ x₁]，则构建A、B、V和状态反馈x包括：

用r^T中的a_i构造矩阵 $A=\left[\begin{array}{cc}0& -a_2\\ 1& -a_1\end{array}\right],$ 用r^T中的b_i构造矩阵 $B=\left[\begin{array}{c}{b}_2\\ b_1\end{array}\right],$ 用η^T中的c_i构造矩阵 $V=\left[\begin{array}{c}{c}_2\\ c_1\end{array}\right],$ 用状态x^T作状态反馈。

步骤1032：利用构建的A、B、V构建状态反馈预测控制方法公式(9)，所述公式(9)如下：

$\left\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Ve}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+e\left(k\right)\end{array}\right.---\left(9\right)$

步骤1033：根据参考轨迹y^s(k+p)＝αy(k)+(1-α)ω、性能指标 $J_p=[\hat{y}(k+p)-y^s(k+p)]Q[\hat{y}(k+p)-y^s(k+p)]+\mathrm{\Δ}{u}^T\left(k\right)\mathrm{\λ\Δu}\left(k\right)$ 和预测公式(9)计算增量形式的控制率公式(10)，所述公式(10)如下：

Δu(k)＝(S^TQS+λ)^-1S^TQ(y^s(k+p)-y(k)-Δ_KX-Δ_SU-Δ_TY)   (10)

其中，ω是输入设定值，α为滤波系数， ${\mathrm{\Δ}}_{\mathrm{KX}}=\left[\begin{array}{c}(C_1{A}^{p_1}-C_1{A}^{p_1}{\mathrm{VC}}_1)(x\left(k\right)-x(k-p_1))\\ (C_2{A}^{p_2}-C_2{A}^{p_2}{\mathrm{VC}}_2)(x\left(k\right)-x(k-p_2))\\ \·\\ \·\\ \·\\ (C_r{A}^{p_r}-C_r{A}^{p_r}{\mathrm{VC}}_r)(x\left(k\right)-x(k-p_r))\end{array}\right],$ ${\mathrm{\Δ}}_{\mathrm{SU}}=\left[\begin{array}{c}\underset{i=1}{\overset{p_1}{\mathrm{\Σ}}}(S_1\left(p_1\right)-S_1\left(i\right))\mathrm{\Δu}(k-i)\\ \underset{i=1}{\overset{p_2}{\mathrm{\Σ}}}(S_2\left(p_2\right)-S_2\left(i\right))\mathrm{\Δu}(k-i)\\ \·\\ \·\\ \·\\ \underset{i=1}{\overset{p_r}{\mathrm{\Σ}}}(S_r\left(p_r\right)-S_r\left(i\right))\mathrm{\Δu}(k-i)\end{array}\right],$ ${\mathrm{\Δ}}_{\mathrm{TY}}=\left[\begin{array}{c}{C}_1{A}^{p_j-1}V(y\left(k\right)-y(k-p_1))\\ C_2{A}^{p_2-1}V(y\left(k\right)-y(k-p_2))\\ \·\\ \·\\ \·\\ C_r{A}^{p_r-1}V(y\left(k\right)-y(k-p_r))\end{array}\right],$ $S\left(p\right)=\left[\begin{array}{c}{S}_1\left(p_1\right)\\ S_2\left(p_2\right)\\ \·\\ \·\\ \·\\ S_r\left(p_r\right)\end{array}\right],$ $y^s=\left[\begin{array}{c}{y}_1^s(k+p_1)\\ y_2^s(k+p_2)\\ \·\\ \·\\ \·\\ y_r^s(k+p_r)\end{array}\right],$ $S_j\left(p_j\right)=\underset{i=1}{\overset{p_j}{\mathrm{\Σ}}}{C}_j{A}^{i-1}B,$ Q为输出预测误差加权阵，λ为控制量的加权系数，Δu(k)为控制作用u(k)的增量；

在本步骤中，由于采用了单值状态反馈预测控制算法，该算法预测值要在一步内调整到给定值，势必造成很大的控制作用，可采用两种手段减少控制作用，一是采用公式y^s(k+p)＝αy(k)+(1-α)ω所述的参考轨迹，一是在优化性能指标中考虑控制作用的影响，即采用公式 $J_p=[\hat{y}(k+p)-y^s(k+p)]Q[\hat{y}(k+p)-y^s(k+p)]+\mathrm{\Δ}{u}^T\left(k\right)\mathrm{\λ\Δu}\left(k\right)$ 所述的性能指标；其中，ω是输入设定值，α为滤波系数，α值越大，***的柔性越好，鲁棒性越强，但控制的快速性却变差，α对闭环***的动态特性和鲁棒性起着重要作用。这里，所述的鲁棒性用于表示***的抗干扰能力，鲁棒性越强，***的抗干扰能力越强。

步骤1034：将Δu(k)＝(S^TQS+λ)^-1S^TQ(y^s(k+p)-y(k)-Δ_KX-Δ_SU-Δ_TY)代入公式u(k)＝u(k-1)+Δu(k)，计算得到当前的控制作用u(k)。

上述步骤103中所述将计算出来的控制作用加到实际被控过程中包括：将计算出来的控制作用通过数模转化，由数字量转化成模拟量，如电压信号，再送到执行元件中。

在将计算出来的控制作用加到实际被控过程后，可以通过控制阀的开度大小调节实际被控过程，使被控量跟随给定值的变化。

另外，由于本发明对于每一个执行周期，都需要采用递推的扩展卡尔曼滤波参数估计方法估计所述利用Singular Pencil Model模型表示的实际被控过程的状态和模型参数，并利用估计出来的状态和模型参数计算当前的控制作用，然后将计算出来的控制作用加到实际被控过程中。所以，在下一个执行周期，需要重复执行步骤102和步骤103，这样就可以实现对实际被控过程的自适应状态反馈预测控制。

基于图1所述的本发明提供的自适应状态反馈预测控制总体技术方案的方法流程图，图2示出了本发明提供的基于噪声估计的自适应状态反馈预测控制结构示意图。

在图2中，各部分模块的功能可通过下述步骤体现：

步骤201：用Singular Pencil Model模型描述被控过程。

步骤202：用扩展卡尔曼滤波器同时估计状态和参数

步骤203：用估计出来的状态和参数构造矩阵A、B和V，用输出y(k)和

估计噪声e(k)，将其代入预测模型 $\left\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Ve}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+e\left(k\right)\end{array}\right.$ 中预测

步骤204：用实际过程输出y(k)和预测模型输出

之差e(k)校正模型预测值

得到y_o(k+j)。

步骤205：给定值ω经过参考模型后得到参考轨迹y^s(k+j)。

步骤206：将历史输入u(k-i)，估计出来的状态

预测的模型输出

和实际的过程输出y(k)等代入性能指标中进行滚动优化。

步骤207：最后根据参考轨迹y^s(k+j)和预测模型的校正输出y_o(k+j)的差值，计算最优控制率作用于被控过程。

步骤208：在下一个周期，重复上述步骤202至步骤207，实现对实际被控过程的自适应状态反馈预测控制。

基于图1所述的方法流程图和图2所述的结构示意图，图3示出了现有技术中GPC方法和本发明提供的控制方法对阶跃扰动的抑制能力仿真对比示意图。

在图3中，假设阶跃扰动的幅值为5，采样时间T_s为0.1秒，两种方法分别是GPC和基于噪声估计的自适应状态反馈预测控制。为衡量闭环***的控制性能，采用误差绝对值之和SAE(Sum of Absolute Value of Error)作为衡量标准，用公式表示为：

$\mathrm{SAE}=\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}|y^s\left(k\right)-y\left(k\right)|\·T_s$

经计算，两种控制方法的SAE分别为53.5010、15.3621，由此可以看出基于噪声估计的自适应状态反馈预测控制方法抑制扰动的能力很强。

为使本发明的目的、技术方案和优点更加清楚明白，下面以质子交换膜燃料电池流量***控制为具体实施例，并参照附图，对本发明进一步详细说明。

燃料电池是一种能将化学能转换成电能的装置。质子交换膜燃料电池由于具有温度低、效率高和无污染的特点而受到广泛关注。它以氢气和氧气，或氢气和空气为原料，反应副产物是水。在氢氧燃料电池中，氢和氧在各自的电极反应。氢气在阳极催化剂的作用下电离成电子和氢离子，发生还原反应，负电子经外电路到达导电的阴极，同时余下的氢离子通过质子渗透膜被送到阴极。在阴极，氢离子与氧气发生氧化反应并从阴极吸收电子生成水。其化学反应方程式为：

阴极：1/2O₂+2H⁺+2e^-→H₂O

阳极：H₂→2H⁺+2e^-

总反应：1/2O₂+H₂→H₂O

在本实施例中，以中国科学院自动化研究所1KW质子交换膜燃料电池流量***为被控对象，流量***的控制目标是通过自动和快速地调节阴极和阳极进气流量使电堆的压力保持在合适的数值，同时使空气和氢气的进气流量快速地跟踪***负载的变化。

如图4所示，图4为本发明实施例中燃料电池流量***控制Simulink模块示意图。整套***包括的硬件主要有燃料电池发电装置，网络化控制器(NetCon)***，负载电动摩托车和一台上位机。这里，Simulink是Matlab软件中的一个工具箱，此外，Matlab软件中还包括实时工具箱Real-Time Workshop。

Netcon***是一个基于32位ARM核微处理器的网络化控制器。其上可运行uclinux操作***。整个***包括可视化的控制组态软件和可视化的监视组态软件，运行于上位机，上位机与控制器之间通过以太网进行通讯。控制组态软件基于Matlab、Simulink和Real-Time Workshop，能够实现和Simulink的无缝结合。

控制算法可在上位机的Matlab/Simulink环境编写，利用Simulink提供的Real-Time Workshop工具和交叉编译器可将Simulink模块图生成可以在Netcon控制器上直接运行的代码，然后通过以太网下载到高级精简指令计算机处理器(Advanced RISC Machines，ARM)上执行。核心算法通过C-Sfunction实现，具体的实现步骤如图5所示，包括以下步骤：

步骤501：利用Singular Pencil Model模型表示实际被控过程；

具体过程如下：置初值，包括将P₀、协方差阵R₀、预测时域P，模型阶次n、参考轨迹滤波***α、输出预测误差加权阵Q、控制量的加权系数λ和确定待估计的参数等置初值。

步骤502：利用公式(8)递推扩展卡尔曼滤波器在线估计

，包括状态信息、***矩阵参数和噪声参数。

步骤503：利用

构造预测模型参数

步骤504：利用

估计出来的状态反馈和公式(10)计算Δu(k)。

步骤505：计算u(k)＝u(k-1)+Δu(k)，加到被控过程中。

步骤506：在下一个执行周期，重复执行步骤502至步骤505，实现对实际被控过程的自适应状态反馈预测控制。

基于图4所示的燃料电池流量***控制Simulink模块示意图和图5所示的燃料电池流量***控制方法流程图，图6示出了本发明实施例中燃料电池流量***的实际控制效果示意图。

在图6中，假设预测时域取6，模型阶次n为2，参考轨迹滤波***α为0.5，输出预测误差加权阵Q为1，控制量的加权系数λ为500，采样时间取0.2秒。其控制策略是对阴极压力采用基于噪声估计的自适应状态反馈预测控制，阳极跟随阴极压力的变化采用简单的PI控制，PID参数为比例0.2，积分为0.001，阴极压力的给定值为25KPa。整个控制策略用Matlab的Simulink模块图实现。

从图6所示的实际控制效果可以看出，本发明提供的这种基于噪声估计的自适应状态反馈预测控制方法非常实用。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1、一种基于噪声估计的自适应状态反馈预测控制方法，其特征在于，该方法包括：

a、利用Singular Pencil Model模型表示实际被控过程；

假设x_k是状态变量，u_k是输入变量，y_k是输出变量，e_k为噪声变量，E_*，A_*，B_*，C_*，E₀为矩阵，I为单位阵，z为前向移位算子，则建立公式

表示实际被控过程；

假设

$w_k={[y_k^T,u_k^T]}^T,$ ${\tilde{G}}_{*}\left(w_k\right)r=\left[{-A}_{*}\right|B_{*}]w_k,$ ${\tilde{G}}_0\left(w_k\right)r=\left[{-A}_0\right|B_0]w_k+{\mathrm{Jw}}_k,$ $C_{*}{e}_k={\tilde{C}}_{*}\left(e_k\right){\mathrm{\η}}_k,$ 其中矩阵

和

根据r中参数的特定位置构建，矩阵

根据η中参数的特定位置构建，则将用于表示实际被控过程的公式

变形为公式 $\left\{\begin{array}{c}{x}_{k+1}=E_{*}{x}_k+{\tilde{G}}_{*}\left(w_k\right)r+C_{*}{e}_k\\ 0=E_0{x}_k+{\tilde{G}}_0\left(w_k\right)r-{\mathrm{Jw}}_k+e_k\end{array}\right.;$

假设 ${\tilde{F}}_k=\left[\begin{array}{ccc}{E}_{*}& {\tilde{G}}_{*}\left(w_k\right)& {\tilde{C}}_{*}\left(e_k\right)\\ 0& I& 0\\ 0& 0& I\end{array}\right],$ $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T,$ ${\tilde{H}}_k=\left[\begin{array}{ccc}{E}_0& {\tilde{G}}_0\left(w_k\right)& 0\end{array}\right],$ 则将公式 $\left\{\begin{array}{c}{x}_{k+1}=E_{*}{x}_k+{\tilde{G}}_{*}\left(w_k\right)r+C_{*}{e}_k\\ 0=E_0{x}_k+{\tilde{G}}_0\left(w_k\right)r-{\mathrm{Jw}}_k+e_k\end{array}\right.;$ 变形为公式 $\left\{\begin{array}{c}{s}_{k+1}={\tilde{F}}_k{s}_k\\ y_k={\tilde{H}}_k{s}_k+e_k\end{array}\right.;$

其中，I是单位矩阵；w_k是由输入输出数据组成的矩阵；r是待估计的***模型参数向量；η是待估计的噪声模型参数向量；

b、采用递推的扩展卡尔曼滤波参数估计方法估计所述利用SingularPencil Model模型表示的实际被控过程的状态和模型参数；

所述递推的扩展卡尔曼滤波参数估计算法公式为

$\left\{\begin{array}{c}{K}_k={\mathrm{\Ξ}}_k{P}_k{\tilde{H}}_k^T{({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)}^{-1}\\ P_{k+1}={\mathrm{\Ξ}}_k{P}_k{\mathrm{\Ξ}}_k^T-K_k({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)K_k^T\\ {\hat{s}}_{k+1}={\tilde{F}}_k{\hat{s}}_k+K_k(y_k-{\tilde{H}}_k{\hat{s}}_k)\end{array}\right.,$

式中， ${\hat{R}}_{k+1}={\hat{R}}_k+\frac{1}{k+1}({\hat{e}}_k{\hat{e}}_k^T-{\hat{R}}_k),{\mathrm{\Ξ}}_k=\frac{\∂\left({\tilde{F}}_k{s}_k\right)}{\∂s_k}{|}_{s_k={\hat{s}}_k},$

$\frac{\∂\left({\tilde{F}}_k{s}_k\right)}{\∂s_k}=\left[\begin{array}{ccc}{E}_{*}-C_{*}{E}_0& {\tilde{G}}_{*}\left(w_k\right)-C_{*}{\tilde{G}}_0\left(w_k\right)& {\tilde{C}}_{*}\left(e_k\right)\\ 0& I& 0\\ 0& 0& I\end{array}\right]$

公式 $P_{k+1}={\mathrm{\Ξ}}_k{P}_k{\mathrm{\Ξ}}_k^T-K_k({\tilde{H}}_k{P}_k{\tilde{H}}_k^T+R_k)K_k^T$ 用于修正K_k，将修正后的K_k代入公式 ${\hat{s}}_{k+1}={\tilde{F}}_k{\hat{s}}_k+K_k(y_k-{\tilde{H}}_k{\hat{s}}_k),$ 计算得到包含所述实际被控过程状态和模型参数的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T;$

c、利用估计出来的状态和模型参数计算当前的控制作用，并将计算出来的控制作用加到实际被控过程中；其中所述利用估计出来的状态和模型参数计算当前的控制作用，具体包括：

c1、利用包含实际被控过程状态和模型参数的参数 $s_k={\left[\begin{array}{ccc}{x}_k^T& r_k^T& {\mathrm{\η}}_k^T\end{array}\right]}^T$ 构建矩阵A、矩阵B、矩阵V和状态反馈向量x；

c2、利用构建的矩阵A、B和V构建状态反馈预测控制方法公式

$\left\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Ve}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+e\left(k\right)\end{array}\right.;$

c3、根据参考轨迹y^s(k+p)＝αy(k)+(1-α(ω、性能指标 $J_p=[\hat{y}(k+p)-y^s(k+p)]Q[\hat{y}(k+p)-y^s(k+p)]+{\mathrm{\Δu}}^T\left(k\right)\mathrm{\λ\Δu}\left(k\right)$ 和预测公式 $\left\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bu}\left(k\right)+\mathrm{Ve}\left(k\right)\\ y\left(k\right)=\mathrm{Cx}\left(k\right)+e\left(k\right)\end{array}\right.,$ 计算增量形式的最优控制率；

Δu(k)＝(S^T QS+λ)^-1S^T Q(y^s(k+p)-y(k)-Δ_KX-Δ_SU-Δ_TY)，

其中，ω是输入设定值，α为滤波系数， ${\mathrm{\Δ}}_{\mathrm{KX}}=\left[\begin{array}{c}(C_1{A}^{p_1}-C_1{A}^{p_1}{\mathrm{VC}}_1)(x\left(k\right)-x(k-p_1))\\ (C_2{A}^{p_2}-C_2{A}^{p_2}{\mathrm{VC}}_2)(x\left(k\right)-x(k-p_2))\\ \·\\ \·\\ \·\\ (C_r{A}^{p_r}-C_r{A}^{p_r}{\mathrm{VC}}_r)(x\left(k\right)-x(k-p_r))\end{array}\right],$ ${\mathrm{\Δ}}_{\mathrm{SU}}=\left[\begin{array}{c}\underset{i=1}{\overset{p_1}{\mathrm{\Σ}}}(S_1\left(p_1\right)-S_1\left(i\right))\mathrm{\Δu}(k-i)\\ \underset{i=1}{\overset{P_2}{\mathrm{\Σ}}}(S_2\left(p_2\right)-S_2\left(i\right))\mathrm{\Δu}(k-i)\\ \·\\ \·\\ \·\\ \underset{i=1}{\overset{p_r}{\mathrm{\Σ}}}(S_r\left(p_r\right)-S_r\left(i\right))\mathrm{\Δu}(k-i)\end{array}\right],$ ${\mathrm{\Δ}}_{\mathrm{TY}}=\left[\begin{array}{c}{C}_1{A}^{p_j-1}V(y\left(k\right)-y(k-p_1))\\ C_2{A}^{p_2-1}V(y\left(k\right)-y(k-p_2))\\ \·\\ \·\\ \·\\ (C_r{A}^{p_r-1}V(y\left(k\right)-y(k-p_r))\end{array}\right],$ $S\left(p\right)=\left[\begin{array}{c}{S}_1\left(p_1\right)\\ S_2\left(p_2\right)\\ \·\\ \·\\ \·\\ S_r\left(p_r\right)\end{array}\right],$ $y^s=\left[\begin{array}{c}{y}_1^s(k+p_1)\\ y_2^s(k+p_2)\\ \·\\ \·\\ \·\\ y_r^s(k+p_r)\end{array}\right],$ $S_j\left(p_j\right)=\underset{i=1}{\overset{p_j}{\mathrm{\Σ}}}{C}_j{A}^{i-1}B,$ Q为输出预测误差加权阵，λ为控制量的加权系数，Δu(k)为控制作用u(k)的增量；

c4、将Δu(k)＝(S^T QS+λ)^-1S^T Q(y^s(k+p)-y(k)-Δ_KX-Δ_SU-Δ_TY)代入公式u(k)＝u(k-1)+Δu(k)，计算得到当前的控制作用u(k)；

d、在下一个执行周期，重复执行步骤b和c。

2、根据权利要求1所述的方法，其特征在于，所述被控过程用二阶单入单出模型表示，r^T＝[a₂ a₁ b₂ b₁]，η^T＝[c₂ c₁]，x^T＝[x₂ x₁]，所述步骤c1包括：

用r^T中的a_i构造矩阵 $A=\left[\begin{array}{cc}0& {-a}_2\\ 1& {-a}_1\end{array}\right],$ 用r^T中的b_i构造矩阵 $B=\left[\begin{array}{c}{b}_2\\ b_1\end{array}\right],$ 用η^T中的c_i构造矩阵 $V=\left[\begin{array}{c}{c}_2\\ c_1\end{array}\right],$ 用状态x^T作状态反馈。

3、根据权利要求1所述的方法，其特征在于，步骤c中所述将计算出来的控制作用加到实际被控过程中包括：

将计算出来的控制作用通过数模转化，由数字量转化成模拟量，再送到执行元件中。

Images