[发明专利]一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法

摘要

本发明公开了一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法。该方法首先依据CSTR连续系统模型，获得了状态分量包含待辨识参数的状态空间表达式；接着，借助欧拉算法对获取的非线性连续状态空间表达式进行了离散化处理，得到了相应的离散迭代模型。最后，运用无迹卡尔曼滤波算法进行多次迭代辨识，获得了准确的辨识结果。该算法收敛性好，且易于与已有软件相结合，具有很好的工程应用前景。

主权项

一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法，其特征在于，包含如下步骤：(1)、获取状态分量中包含CSTR模型待辨识参数的状态空间表达式。(2)、运用欧拉算法对连续的状态空间表达式进行离散化，获得离散的状态空间表达式。(3)、初始化，包括：设定参数辨识的初值初始参数辨识误差协方差以及过程噪声和量测噪声所满足的协方差矩阵Q和R，算法迭代次数最大值L。(4)、选取k‑1时刻的sigma点，计算公式为：${\mathrm{\δ}}_{k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k-1}& {\hat{x}}_{k-1}+\mathrm{\γ}\·\sqrt{{\hat{p}}_{x,k-1}}& {\hat{x}}_{k-1}-\mathrm{\γ}\·\sqrt{{\hat{p}}_{x,k-1}}\end{array}\right]$$\mathrm{\γ}=\sqrt{n+\mathrm{\λ}},\mathrm{\λ}={\mathrm{\α}}^2\·(n+k_f)-n$式中，表示k‑1时刻的状态估计值，表示k‑1时刻的状态估计误差协方差，γ表示尺度参数，n表示的维数；常数α决定sigma点围绕均值的波动范围，通常α∈[10^‑4,1]；常数k_f是另一尺度参数，用于状态估计和参数辨识时通常取0。(5)、在上一步基础上，计算k‑1时刻的sigma点的增值点，计算公式为：${\left(f_{k-1}^{*}\right)}_i=f({\left({\mathrm{\δ}}_{k-1}\right)}_i,u_{k-1})$${\left(h_{k-1}^{*}\right)}_i=h({\left({\mathrm{\δ}}_{k-1}\right)}_i,u_{k-1})$式中，f(·)是对应具体问题系统方程的非线性函数，h(·)是对应具体问题输出方程中的非线性函数，u_k‑1是k‑1时刻输入控制矩阵，下标i表示对应于第i个sigma点的相关取值，i＝0…2n。(6)、计算k‑1时刻的状态向量均值和协方差，计算公式为：${\hat{x}}_{k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\left(f_{k-1}^{*}\right)}_i$${\hat{p}}_{x,k-1}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{w}_i^c({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-}){({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-})}^T+Q_{k-1}$式中，表示k‑1时刻的状态向量均值，表示k‑1时刻的状态协方差，Q_k‑1表示k‑1时刻系统噪声所满足的协方差矩阵，权重系数和取值的计算公式如下：$w_0^m=\frac{\mathrm{\λ}}{n+\mathrm{\λ}},w_0^c=\frac{\mathrm{\λ}}{n+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β})$$w_i^m=w_i^c=\frac{1}{2(n+\mathrm{\λ})},i=1...2n$式中β通常是包含x分布的先验知识，对高斯分布来说，其最优值一般取2。(7)、计算k‑1时刻量测向量均值和协方差，计算公式为：$y_{k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^m{\left(h_{k-1}^{*}\right)}_i$${\hat{p}}_{y,k-1}^{-}=\underset{i=0}{\overset{2n}{\mathrm{\Σ}}}{w}_i^c({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-}){({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-})}^T+R_{k-1}$式中表示k‑1时刻量测向量均值，表示k‑1量测向量的协方差，R_k‑1表示k‑1时刻的量测噪声所满足的协方差矩阵。(8)、计算交互协方差，计算公式如下：${\hat{p}}_{xy,k-1}^{-}=\underset{i=0}{\overset{2n}{\Σ}}{w}_i^c({\left(f_{k-1}^{*}\right)}_i-{\hat{x}}_{k-1}^{-}){({\left(h_{k-1}^{*}\right)}_i-y_{k-1}^{-})}^T$式中，表示k‑1时刻的交互协方差。(9)、在上一步的基础上，计算k‑1时刻的卡尔曼滤波增益，其遵循的计算公式为：$K_{k-1}={\hat{p}}_{xy,k-1}^{-}\·{\left({\hat{p}}_{y,k-1}^{-}\right)}^{-1}$式中，K_k‑1表示k‑1时刻的卡尔曼滤波增益。(10)、运用无迹卡尔曼滤波更新步，获得k时刻的状态估计值和协方差，计算公式为：${\hat{x}}_k={\hat{x}}_{k-1}^{-}+K_{k-1}\·(y_{k-1}-y_{k-1}^{-})$${\hat{p}}_{x,k}={\hat{p}}_{x,k-1}^{-}-K_{k-1}\·{\hat{p}}_{y,k-1}^{-}\·K_{k-1}^T$式中，表示k时刻的状态估计值，y_k‑1表示k‑1时刻量测输出真实值，表示k时刻的估计协方差。(11)、依据上述步骤进行多次迭代辨识，直至k≥L时，结束迭代过程，输出辨识结果。