本發(fā)明涉及一種多速率采樣系統(tǒng)的滾動時域估計方法，尤其涉及一種多速率采樣連續(xù)攪拌釜式反應(yīng)器的滾動時域估計方法。

背景技術(shù)：
連續(xù)攪拌釜式反應(yīng)器（ContinuousStirredTankReactor，簡稱為CSTR）是聚合化學(xué)反應(yīng)中廣泛使用的一種反應(yīng)器。CSTR在化工生產(chǎn)的核心設(shè)備中占有相當重要的地位，在染料、醫(yī)藥試劑、食品及合成材料工業(yè)中，CSTR得到了廣泛的應(yīng)用。在早期反應(yīng)釜的自動控制中，一般采用單元組合儀表組成的位置式控制裝置，但是化學(xué)反應(yīng)過程一般都有很強的非線性和時滯性，采用這種簡單的控制方式很難達到理想的控制精度。隨著計算機技術(shù)的發(fā)展，越來越多的化學(xué)反應(yīng)采用計算機控制系統(tǒng)。計算機控制系統(tǒng)需要使用傳感器采集離散信息，設(shè)計離散控制器。在聚合化學(xué)反應(yīng)中，反應(yīng)釜內(nèi)的溫度和壓力能夠通過傳感器測量快速得到。反應(yīng)物的分子量及濃度則采用凝膠滲透色譜法得到。相比于溫度和壓力的測量周期，分子量和濃度的測量周期將會更長。所以，CSTR中的對溫度和濃度的測量采樣頻率是不相同的，即多速率采樣。對于多速率采樣，在不同的采樣時刻得到的測量輸出數(shù)據(jù)量是不同的，并會導(dǎo)致測量信息不完整。因此需要利用有限的測量信息來估計系統(tǒng)狀態(tài)，反應(yīng)系統(tǒng)的實時性。

技術(shù)實現(xiàn)要素：
本發(fā)明的目的在于提供一種適用于具有多速率采樣連續(xù)攪拌釜式反應(yīng)器的滾動時域估計方法。本發(fā)明解決其技術(shù)問題所采用的技術(shù)方案是：一種具有多速率采樣連續(xù)攪拌釜式反應(yīng)器的滾動時域估計方法，具體步驟如下：（1）、對CSTR裝置進行狀態(tài)空間建模，并確定其正常工作的工作范圍和工作范圍內(nèi)的穩(wěn)態(tài)工作點。（2）、在工作范圍內(nèi)的穩(wěn)態(tài)工作點附近將CSTR的狀態(tài)空間模型線性化。（3）、設(shè)定系統(tǒng)的標準采樣周期和各個傳感器的采樣周期，將線性化的CSTR狀態(tài)空間模型離散化，對于未采樣的測量輸出用預(yù)測值補償，得到多速率的CSTR線性離散化狀態(tài)空間模型。（4）、設(shè)定滾動時域窗口長度N和權(quán)重矩陣，將多速率CSTR的滾動時域估計問題轉(zhuǎn)化為等價的最小化問題。（5）、通過一階最優(yōu)性原理求解步驟（4）設(shè)定的最小化問題，具體步驟如下：S1-1：初始化，設(shè)定測試時間長度K，在可行域的區(qū)間范圍內(nèi)，任意初始化k時刻的先驗估計值，k-N時刻到k時刻的測量輸出序列；S1-2：根據(jù)一階最優(yōu)性原理，以k時刻的先驗估計值為初始迭代點，對步驟（4）中的最小化問題求一階偏導(dǎo)數(shù)，得到k-N時刻最優(yōu)估計值；S1-3：根據(jù)滾動優(yōu)化原理，計算當前k時刻的最優(yōu)估計值；S1-4：根據(jù)k時刻的最優(yōu)估計值更新k+1時刻的先驗估計值；S1-5：判斷終止條件：如果k=K，結(jié)束，得到信噪比估計最優(yōu)值；否則，k=k+1，轉(zhuǎn)到S1-2。本發(fā)明的技術(shù)構(gòu)思為：本發(fā)明考慮了CSTR內(nèi)聚合反應(yīng)過程中對反應(yīng)釜內(nèi)的溫度和反應(yīng)物濃度測量頻率的差異，給出了一種具有多速率采樣的CSTR線性離散化狀態(tài)空間模型，設(shè)計了基于滾動時域估計方法的CSTR狀態(tài)估計器，給出了反應(yīng)釜內(nèi)溫度和反應(yīng)物濃度的最優(yōu)估計值。從上述技術(shù)方案可以看出，本發(fā)明的有益效果主要表現(xiàn)在：多速率CSTR滾動時域估計方法，與已有的估計方法相比，滾動時域估計方法能夠滾動優(yōu)化和在線計算，并對傳感器沒有采樣的時刻用預(yù)測值替代起到很好的補償作用，從而能夠更加準確地給出反應(yīng)釜內(nèi)各參數(shù)的值。附圖說明圖1是本發(fā)明實施例中CSTR示意圖。圖2是本發(fā)明實施例中求解最小化問題的流程圖。圖3是本發(fā)明實施例中，采用本發(fā)明方法的效果圖。具體實施方式為使本發(fā)明的目的、技術(shù)方案和優(yōu)點更加清晰，下面結(jié)合附圖和實施例對本發(fā)明的技術(shù)方案作進一步描述。參照圖1～圖3，一種具有多速率采樣CSTR的滾動時域估計方法，將本發(fā)明提出的滾動時域估計方法用于一階并行反應(yīng)CSTR的狀態(tài)估計，其目的是估計出反應(yīng)釜內(nèi)的反應(yīng)物濃度和溫度。一階并行反應(yīng)，即往CSTR中加入某物質(zhì)A，進入反應(yīng)釜后，A發(fā)生化學(xué)反應(yīng)生成物B，同時由于反應(yīng)釜內(nèi)的強反應(yīng)，部分物質(zhì)B生成物質(zhì)C，并達到化學(xué)平衡的過程。接下來介紹具體實施步驟：（1）、對附圖1所示一階并行反應(yīng)CSTR裝置進行狀態(tài)空間建模，并確定其正常工作的工作范圍和工作范圍內(nèi)的穩(wěn)態(tài)工作點。如附圖1所示，建立一階并行反應(yīng)CSTR的狀態(tài)空間模型如下所示：dCAdt=FV(CA0-CA)-k1CAdCBdt=-FVCB+k1CA-k2CBdθdt=FV(θ0-θ)+kwARρCPV(θk-θ)-k1CAΔHRAB+k2CBΔHRBCρCP---(1)]]>式中，CA為物質(zhì)A的濃度，CA0為物質(zhì)A初始濃度，CAs為物質(zhì)A在穩(wěn)態(tài)時濃度，CB為物質(zhì)B的濃度，CBs為物質(zhì)B在穩(wěn)態(tài)時濃度，θ為反應(yīng)釜內(nèi)溫度，θ0為反應(yīng)釜內(nèi)初始溫度，θs為穩(wěn)態(tài)時反應(yīng)釜內(nèi)溫度，θk為冷卻劑溫度，F(xiàn)/V為稀釋率，V為體積流量，AR為反應(yīng)堆表面面積，CP為熱容量，kw為傳熱系數(shù)，ρ為密度，為物質(zhì)A到物質(zhì)B反應(yīng)焓，為物質(zhì)B到物質(zhì)C反應(yīng)焓，反應(yīng)速率系數(shù)k1和k2由反應(yīng)釜內(nèi)溫度決定i=1,2，k0為頻率因子，EA1和EA2為活化能，R為理想氣體常數(shù)。CSTR的工作范圍和穩(wěn)態(tài)工作點取值如表1所示：表1.CSTR模型參數(shù)和穩(wěn)態(tài)工作點（2）、根據(jù)一階并行反應(yīng)CSTR的狀態(tài)空間模型，在工作范圍內(nèi)的穩(wěn)態(tài)工作點附近將CSTR的狀態(tài)空間模型線性化，得到線性化狀態(tài)空間模型如下所示：x·(t)=Ax(t)+w(t)---(2)]]>式中，x(t)=[x1(t)x2(t)x3(t)]T，x1(t)為物質(zhì)A在t時刻的濃度，x2(t)為物質(zhì)B在t時刻的濃度，x3(t)為在t時刻反應(yīng)釜內(nèi)的溫度，w(t)為在t時刻反應(yīng)釜內(nèi)的有界擾動即||w(t)||≤0.3，A為線性化后得到的系統(tǒng)參數(shù)A=-FV-k10EA1Rθs2k1CAsk1-FV-k2-EA1Rθs2k1CAs+EA2Rθs2k2CBs-k1ΔHRABρCP-k2ΔHRBCρCP-FV-kwARρCPV+EA1k1CAsΔHRAB+EA2k2CBsΔHRBCRθs2ρCP.]]>=-0.938800.04590.625-0.9388-0.0125-0.93352.4449-0.8894]]>（3）、設(shè)定系統(tǒng)的標準采樣周期T0=1min，CSTR測量輸出量為物質(zhì)B的濃度和反應(yīng)釜內(nèi)溫度，對物質(zhì)A濃度的測量周期為T1=2min，對反應(yīng)釜內(nèi)溫度的測量周期為T2=1min。以標準采樣周期T0將線性化的CSTR狀態(tài)空間模型離散化，得到離散狀態(tài)空間模型為x(k+1)=A‾x(k)+w(k)---(3)]]>式中，A‾=eAT0=0.38720.02220.01820.24440.38970.0007-0.06850.97110.4008.]]>由于對物質(zhì)A濃度的測量為2倍的標準采樣周期，從而會導(dǎo)致未采樣時刻得不到測量值，這種情況下可將未采樣時刻的數(shù)據(jù)視為測量信息丟失，并采用預(yù)測值替代進行輸出補償。由此可得CSTR的測量輸出方程如下所示：y(k)=θ(k)[Cx(k)+v(k)]+[I-θ(k)]y‾(k)---(4)]]>式中，y(k)=[y1(k)y2(k)]T，y1(k)為物質(zhì)A在k時刻濃度的測量值，y2(k)為在k時刻反應(yīng)釜內(nèi)溫度的測量值，v(k)為在k時刻反應(yīng)釜內(nèi)的有界擾動即||v(k)||≤0.2，為k時刻系統(tǒng)狀態(tài)x(k)的預(yù)測值，C=100001]]>為測量輸出權(quán)重矩陣，I=1001]]>為單位矩陣，θ(k)=θ1(k)00θ2(k),]]>θ1(k)=1表示物質(zhì)A在k時刻有測量值，θ1(k)=0表示物質(zhì)A在k時刻沒有測量值，θ2(k)=1。（4）、設(shè)定滾動時域窗口長度N=5和權(quán)重矩陣μ=0.1，多速率CSTR的滾動時域估計問題如下所示：minx^(k-N)J(k)---(5)]]>約束條件：J(k)=||x^(k-N)-x‾(k-N)||μ2+Σi=k-Nk||y(i)-θ(i)[Cx^(i)+v(i)]+[I-θ(i)]y‾(i)||2]]>x^(i+1)=Ax^(i),i=k-N,···,k-1]]>x‾(k-N)=Ax^(k-N-1),k=N+1,N+2,···]]>式中，J(k)為性能指標，為k時刻x(k)的最優(yōu)估計值。（5）、通過一階最優(yōu)性原理求解步驟（4）設(shè)定的最小化問題（5），具體步驟如下：S1-1：初始化，設(shè)定測試時間長度K=80，在可行域的區(qū)間范圍內(nèi)，任意初始化k時刻的先驗估計值x‾(0)=110T]]>和系統(tǒng)狀態(tài)初始值x(0)=[0.50.50.5]T，k-5時刻到k時刻的測量輸出序列yk-5k=y(k-5)...y(k-1)y(k);]]>S1-2：根據(jù)一階最優(yōu)性原理，以k時刻的先驗估計值為初始迭代點，對步驟（4）中的最小化問題求一階偏導(dǎo)數(shù)，得到k-5時刻的最優(yōu)估計值如下所示x^(k-5)=(μI3+G1)-1(μx‾(k-5)+G2yk-5k)]]>式中，G1=F5TF5-F15T(I5-Θ(k))F15,]]>G2=F5T-F15T(I5-Θ(k)),]]>I3為3維的單位矩陣，I5為12維的單位矩陣，Θ(k)=diag{θ(k-5),…,θ(k-1),θ(k)}，F(xiàn)5=CCA...CA5=1000010.38720.02220.0182-0.06850.97110.40080.15410.0350.01440.18340.76610.16010.06720.0310.00860.24730.45810.06810.0330.02190.00470.20310.25010.03210.01780.01380.00250.13760.13310.0168,]]>F15=C0...0=100001000000000000000000000000000000;]]>S1-3：根據(jù)滾動優(yōu)化原理，計算當前k時刻的最優(yōu)估計值S1-4：根據(jù)k時刻的最優(yōu)估計值更新k+1時刻的先驗估計值x‾(k-4)=Ax^(k-5);]]>S1-5：判斷終止條件：如果k=K，結(jié)束，得到信噪比估計最優(yōu)值；否則，k=k+1，轉(zhuǎn)到S1-2。采用所述步驟，計算80個采樣時刻得到的結(jié)果如圖3所示，其中，圖3中的(a)為CSTR物質(zhì)A的濃度和其估計值，橫坐標為采樣次數(shù)，縱坐標為物質(zhì)A的濃度和其估計值；圖3中的(b)為CSTR物質(zhì)B的濃度和其估計值，橫坐標為采樣次數(shù)，縱坐標為物質(zhì)B的濃度和其估計值；圖3中的(c)為CSTR內(nèi)溫度和其估計值，橫坐標為采樣次數(shù)，縱坐標為溫度和其估計值。從圖3可以看出，針對具有多速率采樣的CSTR，滾動時域估計方法對未采樣部分信息缺失也能夠很好的估計出系統(tǒng)的狀態(tài)。并且滾動時域估計方法具有滾動優(yōu)化和在線計算的優(yōu)點，從而能夠更加準確地給出CSTR反應(yīng)物濃度和反應(yīng)器溫度。以上闡述的是本發(fā)明給出的實例表現(xiàn)出的優(yōu)良估計效果。需要指出的是，本發(fā)明不只限于上述實施例，對于其他化工聚合反應(yīng)多速率采樣估計問題，采用本發(fā)明給出的方法設(shè)計滾動時域估計器，均能給出系統(tǒng)的狀態(tài)估計值。