一种基于无迹卡尔曼滤波算法的cstr模型参数辨识方法
【专利摘要】本发明公开了一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法。该方法首先依据CSTR连续系统模型，获得了状态分量包含待辨识参数的状态空间表达式；接着，借助欧拉算法对获取的非线性连续状态空间表达式进行了离散化处理，得到了相应的离散迭代模型。最后，运用无迹卡尔曼滤波算法进行多次迭代辨识，获得了准确的辨识结果。该算法收敛性好，且易于与已有软件相结合，具有很好的工程应用前景。
【专利说明】
一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法
技术领域
[0001] 本发明涉及一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法，属于系统建 模与参数辨识技术领域。
【背景技术】
[0002] 连续搅拌釜式反应器(CSTR)是化工生产过程中典型的、高度非线性的化学反应系 统。在化工生产的核心设备中占有相当重要的地位，在染料、医药试剂、食品及合成材料工 业中，CSTR系统得到了广泛的应用。
[0003] 由于CSTR在化工生产中的重要作用，所以有必要对该过程进行详细的研究。在对 CSTR建立模型进行分析时，模型的参数有时是未知的，所以对CSTR的模型参数辨识方法进 行研究，具有重要的工程意义。然而，已有的方法如扩展卡尔曼滤波等，这些方法辨识结果 有时会出现发散，得不到正确的结果。为了提高辨识效率和精度，研究运用无迹卡尔曼滤波 算法进行CSTR模型参数辨识，具有重要的意义。
【发明内容】
[0004] 为了有效的了解化工生产中CSTR化学反应系统，本发明提出了一种基于无迹卡尔 曼滤波算法的CSTR模型参数辨识方法，有效的实现了CSTR模型的参数辨识。
[0005] 本发明的技术解决方案是:一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方 法，其步骤如下：
[0006] (1)、获取扩维状态向量中包含CSTR模型待辨识参数的状态空间表达式；
[0007] (2)、运用欧拉算法对连续的状态空间表达式进行离散化，获得离散化的状态空间 表达式；
[0008] (3)、初始化，包括:设定参数辨识的初值V，初始参数辨识误差协方差.?。，以及过 程噪声和量测噪声所满足的协方差矩阵Q和R，算法迭代次数最大值L;
[0009] (4)、选取k-Ι时刻的sigma点，计算公式为：
[0010]
[0011]
[0012]式中，表示k-Ι时刻的状态估计值，表示k-i时刻的状态估计误差协方差， γ表示尺度参数，11表示毛^勺维数;常数α决定sigma点围绕均值的波动范围，通常ae [10-4，1 ];常数kf是另一尺度参数，用于状态估计和参数辨识时通常取0。
[0013] (5)、在上一步基础上，计算k-Ι时刻的sigma点的增值点，计算公式为：
[0014]
[0015]
[0016] 式中，f( ·)是对应具体问题系统方程的非线性函数，h( ·)是对应具体问题输出 方程中的非线性函数，Uk-^k-l时刻输入控制矩阵，下标i表示对应于第i个sigma点的相关 取值，? = 0···2η。
[0017] (6)、计算k-1时刻的状态向量均值和协方差，计算公式为：
[0018]
[0019]
[0020] 式中，ib.表示k-Ι时刻的状态向量均值，表示k-Ι时刻的状态协方差，Qk-i表示 k-Ι时刻系统噪声所满足的协方差矩阵，权重系数Mf和wf取值的计算公式如下：
[0021]
[0022]
[0023] 式中β通常是包含X分布的先验知识，对高斯分布来说，其最优值一般取2。[0024] (7)、计算k-Ι时刻量测向量均值和协方差，计算公式为：
[0025]
[0026] /-U
[0027] 式中)'a-i表示k-Ι时刻量测向量均值，爲,i-i表示k-Ι量测向量的协方差，Rk-i表示k-1 时刻的量测噪声所满足的协方差矩阵。
[0028] (8)、计算交互协方差，计算公式如下：
[0029]
[0030] 式中，^表示k-Ι时刻的交互协方差。
[0031] (9)、在上一步的基础上，计算k-Ι时刻的卡尔曼滤波增益，其遵循的计算公式为：
[0032] AV,
[0033] 式中，Kk-i表示k-Ι时刻的卡尔曼滤波增益。
[0034] (10)、运用无迹卡尔曼滤波更新步，获得k时刻的状态估计值和协方差，计算公式 为：
[0035]
[0036]
[0037] 式中，4表示k时刻的状态估计值，yk-ι表示k-Ι时刻量测输出真实值，爲』表示k时 刻的估计协方差。
[0038] (11 )、依据上述步骤进行多次迭代辨识，直至k多L时，结束迭代过程，输出辨识结 果。
[0039]本发明的有益效果是:运用无迹卡尔曼滤波对CSTR模型参数辨识，由于无迹卡尔 曼滤波采用了不同于扩展卡尔曼滤波的线性化方法，因此，克服了扩展卡尔曼滤波辨识过 程中出现发散的情况，提高了辨识效率和精度。
【附图说明】
[0040] 图1为本发明实施例的方法流程图；
[0041] 图2为实施例采用本发明所提方法的辨识结果；
[0042] 图3为实施例辨识结果的相对误差；
【具体实施方式】
[0043]下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明 而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价 形式的修改均落于本申请所附权利要求所限定的范围。
[0044] 如图1所示，CSTR模型参数辨识方法，其包含如下步骤：
[0045] (1)、获取状态分量中包含CSTR模型待辨识参数的状态空间表达式。
[0046] (2)、运用欧拉算法对连续的状态空间表达式进行离散化，获得离散化的状态空间 表达式。
[0047] (3)、初始化，包括:设定参数辨识的初值知.，初始参数辨识误差协方差/U，以及过 程噪声和量测噪声所满足的协方差矩阵Q和R，算法迭代次数最大值L。
[0048] (4)、选取k-Ι时刻的sigma点。
[0049] (5)、在上一步基础上，计算k-Ι时刻的sigma点的增值点。
[0050] (6)、计算k-Ι时刻的状态向量均值和协方差。
[0051 ] (7)、计算k-Ι时刻量测向量均值和协方差。
[0052] (8)、计算k-Ι时刻的交互协方差。
[0053] (9)、在上一步的基础上，计算k-Ι时刻的卡尔曼滤波增益。
[0054] (10)、运用无迹卡尔曼滤波更新步，获得k时刻的状态估计值和协方差。
[0055] (11 )、依据上述步骤进行多次迭代辨识，直至k多L时，结束迭代过程，输出辨识结 果。
[0056] 假设在连续搅拌釜式反应器(CSTR)内发生放热的、不可逆反应。反应物为A、生成 物为B。根据物料平衡和能量平衡关系，得到如下的反应过程模型(CSTR状态空间描述）：
[0057]
[0058]该非线性系统具有两个状态变量即：反应器中组分A的浓度CA(状态变量X1)，反应 温度T(状态变量X2);-个控制输入变量T。(输入u)。在下面的实施例仿真分析时假设P和f Μ 是待辨识的参数，它们的真值分别是1000和8750。模型中其他参数的物理意义及取值如表 1〇
[0059]表1 CSTR模型中相关参数的物理意义及取值
[0060]
[0062] 为了运用无迹卡尔曼滤波对模型中的未知参数进行辨识，首先需要获取状态分量 中包含待辨识参数的状态空间表达式，为此假设P为状态分量 X3、|为状态分量X4,则整理后 K 的扩维状态空间表达式为(带入相关的模型参数）：
[0063] 状态方程为
[0064] 输出方程为
[0065] 式中Wi(i = l，2,3,4)是系统噪声，Vi(i = l，2)是量测噪声，它们均为零均值的高斯 白噪声，且分别满足协方差矩阵Qk和Rk，即：
[0066]
[0067] 采用欧拉算法对上述状态方程进行离散化，得到了对应的离散状态方程，在此基 础上，则就可以运用无迹卡尔曼滤波进行多次迭代辨识。在辨识过程中，无迹卡尔曼滤波的 相关参数取值为：
[0068]
[0069]基于上述分析，通过多次迭代辨识，获得了准确的辨识结果。图1为实施例所用的 算法流程图，图2为实施例的参数辨识结果，图3为运用本发明所提出的方法对实施例辨识 结果的相对误差。
【主权项】
1. 一种基于无迹卡尔曼滤波算法的CSTR模型参数辨识方法，其特征在于，包含如下步 骤： (1) 、获取状态分量中包含CST財莫型待辨识参数的状态空间表达式。 (2) 、运用欧拉算法对连续的状态空间表达式进行离散化，获得离散的状态空间表达 式。 (3) 、初始化，包括:设定参数辨识的初值.?。，初始参数辨识误差协方差与.。，^及过程噪 声和量测噪声所满足的协方差矩阵Q和R，算法迭代次数最大值L。 (4) 、选取k-1时刻的Sigma点，计算公式为：式中，為_1表示k-1时刻的状态估计值，表示k-1时刻的状态估计误差协方差，丫表 示尺度参数，η表示成_1的维数;常数α决定Sigma点围绕均值為的波动范围，通常α e [ 1〇-4, 1 ];常数kf是另一尺度参数，用于状态估计和参数辨识时通常取0。 巧）、在上一步基础上，计算k-1时刻的sigma点的增值点，计算公式为：式中，f( ·)是对应具体问题系统方程的非线性函数，h( ·)是对应具体问题输出方程 中的非线性函数，uk-i是k-1时刻输入控制矩阵，下标i表示对应于第i个sigma点的相关取 值，i = 0...2n。 (6) 、计算k-1时刻的状态向量均值和协方差，计算公式为：式中，振康示k-1时刻的状态向量均值，托,W表示k-1时刻的状态协方差，Qk-康示k-1时 刻系统噪声所满足的协方差矩阵，权重系数 <和<取值的计算公式如下：式中防I常是包含X分布的先验知识，对高斯分布来说，其最优值一般取2。 (7) 、计算k-1时刻量测向量均值和协方差，计算公式为：式中.报1表示k-1时刻量测向量均值，馬.W表示k-1量测向量的协方差，Rk-I表示k-1时刻 的量测噪声所满足的协方差矩阵。 (8) 、计算交互协方差，计算公式如下：式中，托表不k-1时刻的交互协方差。 (9) 、在上一步的基础上，计算k-1时刻的卡尔曼滤波增益，其遵循的计算公式为：式中，Kk-i表示k-1时刻的卡尔曼滤波增益。 (10 )、运用无迹卡尔曼滤波更新步，获得k时刻的状态估计值和协方差，计算公式为：式中，為表示k时刻的状态估计值，yk-i表示k-1时刻量测输出真实值，把表示k时刻的估 计协方差。 (11 )、依据上述步骤进行多次迭代辨识，直至L时，结束迭代过程，输出辨识结果。
【文档编号】G06F19/00GK105975747SQ201610272832
【公开日】2016年9月28日
【申请日】2016年4月27日
【发明人】周琪, 于占东, 王焕清, 王巍
【申请人】渤海大学