一種基于無跡卡爾曼濾波算法的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日
【發明人】周琪, 于占東, 王煥清, 王巍
【申請人】渤海大學