過程控制技術(shù)-第5章課件2

1、5.3 測試法建模5.3.1 測試法建模的方法 機理法建模能夠解決實際生產(chǎn)中一部分過程的建模問題。但是，還有很多生產(chǎn)過程由于工藝的復(fù)雜性、產(chǎn)品本身在加工中的變化性（如物理或化學(xué)變化），使得用機理法建模在技術(shù)上遇到了極大的困難，人們不得不考慮用其它的方法建模。 從學(xué)科角度看，建立數(shù)學(xué)模型應(yīng)該屬于系統(tǒng)辨識（System Identification）與參數(shù)估計（Parametric Estimation）的范疇。簡單地說，系統(tǒng)辨識主要是對被研究對象的結(jié)構(gòu)進(jìn)行判斷，解決“是什么”的問題，例如一階慣性環(huán)節(jié)，二階系統(tǒng)；而參數(shù)估計則對支撐結(jié)構(gòu)的參數(shù)進(jìn)行估計，解決“是多少”的問題。 事實上，有很多比較復(fù)雜的

2、過程，我們對其工作機理并不清楚，更難以用數(shù)學(xué)和物理的方法加以具體描述，此時用測試法建模是一個不得已而為之的方法。與前述的機理法建模相比，測試法建模不需要深入了解過程的工作機制，通常的做法是將其看作一個“黑箱”，通過從外部施加適當(dāng)?shù)妮斎胄盘?，測得過程的輸出信號，通過對這些輸入和輸出信號的處理和研究，獲得其動態(tài)特性和數(shù)學(xué)模型。因此，問題主要歸納為： 1）施加何種輸入信號才能最大限度地激勵被測過程，使得動態(tài)特性得以充分表現(xiàn)，并通過輸出信號顯露出來？ 2）對獲得的數(shù)據(jù)或波形，通過什么方法和技術(shù)才能估算出適用于控制用的動態(tài)模型？ 一般說來，模型有非參數(shù)模型（Nonparametric Model）和參數(shù)

3、模型（Parametric Model）之分。建立非參數(shù)模型的方法通常有：時域法（Time-domain Method）、頻域法（Frequency- 5.3 測試法建模5.3.1 測試法建模的方法 -domain Method）和統(tǒng)計相關(guān)法（Statistical Correlation Method）等，這類建模不需要事先確定模型的結(jié)構(gòu)，可用于廣泛的被控過程；獲得參數(shù)模型的方法主要有：最小二乘法（Least Square Method）、極大似然法（Maximum Likelihood Method）和梯度校正法（Gradient Correction Method）等，這類建模需要假設(shè)一

4、定的模型結(jié)構(gòu)，通過極小化模型與過程之間的誤差準(zhǔn)則，來確定相應(yīng)的模型參數(shù)。 用時域法測定被控過程的數(shù)學(xué)模型：對過程施加階躍信號，或者方波信號，測取響應(yīng)曲線，并由此確定過程的傳遞函數(shù)。該方法具有測試簡單、需用設(shè)備少的優(yōu)點，但測試精度不高，其獲得的模型可用于一般工業(yè)過程控制； 用頻域法測定被控過程的數(shù)學(xué)模型：對過程施加不同頻率的正弦波輸入信號，獲得相應(yīng)的輸出幅值與相位，由此可得到該過程的頻率特性，由頻率特性獲得傳遞函數(shù)。該方法需用專門的頻率發(fā)生和測試設(shè)備，模型精度比用時域法高； 用統(tǒng)計相關(guān)法測定被控過程的數(shù)學(xué)模型：對被控過程施加偽隨機信號，采用統(tǒng)計相關(guān)法獲得過程的動態(tài)特性。它的特點是，可在生產(chǎn)狀態(tài)下

5、施加隨機信號，并測取相關(guān)數(shù)據(jù)，精度較高，但需獲得較多數(shù)據(jù)，并借助計算機協(xié)助處理。 最小二乘法又稱最小平方法，是估計離散時間數(shù)學(xué)模型參數(shù)的一常用種方法。隨著計算機技術(shù)在控制中的應(yīng)用，最小二乘法在過程辨識的實踐中被越來越廣泛地采用。 本章主要討論用時域響應(yīng)法和最小二乘法獲取數(shù)學(xué)模型。5.3.2 時域響應(yīng)曲線法 1響應(yīng)曲線的測取 -domain Method）和統(tǒng)計相關(guān)法（Statisti時域響應(yīng)曲線法是對被控過程施加階躍信號，如果被控過程不允許長期施加階躍信號，則改用矩形脈沖信號，然后測取響應(yīng)曲線，并由此求取輸入和輸出之間的傳遞函數(shù)。 （1）階躍響應(yīng)曲線的測取 當(dāng)被控過程穩(wěn)定之后，對調(diào)節(jié)閥施加一個

6、幅值合適的階躍信號，用記錄儀或數(shù)據(jù)采集系統(tǒng)記錄被控量的變化曲線，例如被控量為溫度時，就記錄溫度響應(yīng)曲線，直到變化曲線進(jìn)入穩(wěn)定狀態(tài)為止。下面幾點是在這一過程中值得注意的： a) 在施加階躍信號之前，被控過程應(yīng)處在較為穩(wěn)定的工作狀態(tài)，在下一輪施加輸入信號時，應(yīng)等前一過程結(jié)束、并恢復(fù)穩(wěn)態(tài)一段時間之后進(jìn)行。 b) 施加階躍信號的幅度：通常為正常輸入信號的 515，以不影響正在進(jìn)行的生產(chǎn)為好。同時，幅度也不能太小，因為過小的輸入容易被其它信號淹沒，在響應(yīng)曲線上難以表現(xiàn)出來。 c) 多次、全面測試，消除偶然性，獲得真實結(jié)果：試驗不僅應(yīng)在相同條件下，重復(fù)幾次，以獲得兩次及其以上的較為接近的響應(yīng)曲線，而且也應(yīng)

7、選取不同負(fù)荷、不同輸入值，測得相應(yīng)的響應(yīng)曲線，以獲得全面的動態(tài)特性。 （2）矩形脈沖相應(yīng)曲線的測取 在有些情況下，用階躍信號輸入時，可能會危及安全生產(chǎn)，或者影響產(chǎn)品的質(zhì)量和數(shù)量。此時，輸入可考慮用矩形脈沖信號代替階躍信號，測得過程的矩形脈沖響應(yīng)曲線。由于通過階躍響應(yīng)曲線求傳遞函數(shù)被人們所熟悉，所以，往往將矩形脈沖響應(yīng)曲線再轉(zhuǎn)化為階躍響應(yīng)曲線，進(jìn)而按階躍響應(yīng)曲線法確定傳遞函數(shù)。 圖 5-12a) 為矩形脈沖輸入信號，它可以分解為圖 b) 所示的兩個階躍信號的疊加，即時域響應(yīng)曲線法是對被控過程施加階躍信號，如果被控過程不允許長其中 a 為脈沖寬度， 。如果被控過程是線性的，則其矩形脈沖響應(yīng)曲線 可

8、分解為階躍響應(yīng)曲線 和 ，即圖 5-12 矩形脈沖及其響應(yīng)分解圖 見圖5-12 c)所示。這里， 和 分別為 和 的響應(yīng)。于是，起源于零點的階躍響應(yīng)為： （5-22） 其中， 為矩形脈沖響應(yīng)， 為起源于 點的階躍響應(yīng)（注意此時的符號為正）。式（5-22）為由矩形脈沖響應(yīng)求階躍響應(yīng)的公式，可用作圖法逐步求出： t 在 0a 時， ，當(dāng) 時， ，此時的 前面已經(jīng)有， 為已知，所以可求得 ，如此類推，一直進(jìn)行到進(jìn)入穩(wěn)態(tài)。 其中 a 為脈沖寬度， 由無自平衡能力過程的矩形脈沖響應(yīng)曲線轉(zhuǎn)化為階躍響應(yīng)曲線，也可通過作圖法實現(xiàn)，這可作為練習(xí)，留給讀者來完成（見思考題與習(xí)題 5-5）。 2由過程階躍響應(yīng)曲線確

9、定傳遞函數(shù) 由階躍響應(yīng)曲線確定傳遞函數(shù)，通常需要確定傳遞函數(shù)的結(jié)構(gòu)及其參數(shù)兩部分。結(jié)構(gòu)形式是指被控過程的傳遞函數(shù)形式，生產(chǎn)過程主要是：一階慣性環(huán)節(jié)、二階慣性環(huán)節(jié)、或等 n 階慣性環(huán)節(jié)，并且這些環(huán)節(jié)時常含有純滯后，其表達(dá)形式為 ， ， ， ， 對于無自平衡能力的過程，也有類似的形式 ， ， ， ， 傳遞函數(shù)的參數(shù)是伴隨結(jié)構(gòu)形式出現(xiàn)的待定常數(shù)，如一階慣性、具有純時延、有自平衡能力的傳遞函數(shù)含有：K、T 和 三個需要確定的參數(shù)。 由無自平衡能力過程的矩形脈沖響應(yīng)曲線轉(zhuǎn)化為階 關(guān)于傳遞函數(shù)結(jié)構(gòu)形式的確定，主要有兩方面的考慮：一是根據(jù)對被控過程的經(jīng)驗和知識（即通常所說的先驗知識）來確定；二是根據(jù)控制的要

10、求，盡量將一個原本較復(fù)雜的過程用低階的傳遞函數(shù)來近似描述，因此產(chǎn)生的誤差只要處在可接受的范圍即可。下面的討論，集中在參數(shù)的確定上。 （1）由階躍響應(yīng)曲線求一階慣性加純時延環(huán)節(jié)的參數(shù) 這里傳遞函數(shù)形式為（5-23） 它有三個參數(shù)需要確定，即放大系數(shù) K、時間常數(shù) T 和純時延 。原本就是一階慣性加純時延過程的階躍響應(yīng)曲線見圖 5-5，其 T 和 從圖中很容易確定，放大系數(shù)為 （其中 為階躍輸入信號的幅值）。下面討論原本是二階及其以上過程，且響應(yīng)曲線呈 “s” 形，如何用式（5-23）來近似描述。 現(xiàn)有階躍響應(yīng)曲線如圖 5-13 所示，試圖用式（5-23）來近似描述，需確定 K、T 和 三個參數(shù)。

11、 顯然，放大系數(shù)可用下式求得（5-24） 關(guān)于傳遞函數(shù)結(jié)構(gòu)形式的確定，主要有兩方面的 圖 5-13 由階躍響應(yīng)曲線確定 T 和 圖5-14 縱坐標(biāo)的標(biāo)么化其中 為輸入幅值，為已知量。 關(guān)于 T 和 的確定，有兩種方法：一是作圖法，二是計算法，下面分別介紹。 用作圖法求 T 和 ：首先找到響應(yīng)曲線上凹和下凹的交接點 - 拐點 D，過 D 點作曲線的切線，切線與時間軸 t 相交于 A 點，與 相較于 C 點，該點在時間軸上的投影為 B，則OA為 ，AB 為 T。 作圖法的問題是：曲線的拐點有時不容易找到，并且作切線時，有一定的隨意性。所以，用作圖法求 T 和 ，可能會因人而異，有一定的誤差。 計算

12、法求 T 和 ：將階躍響應(yīng)曲線的縱坐標(biāo) 標(biāo)么化，即：實際值基準(zhǔn)值，這里取基準(zhǔn)值為 ，于是 的標(biāo)么值（Unit Value）為 圖 5-13 由階躍響應(yīng)曲線確定 T 和 前面的圖 5-13 則變?yōu)閳D 5-14。該標(biāo)么化處理，并不改變響應(yīng)曲線的橫坐標(biāo)和形狀，僅方便求得參數(shù)。我們的目的是要通過圖5-14所示的階躍響應(yīng)曲線，容易求出式（5-23）中的 T 和 。 式（5-23）的階躍響應(yīng)標(biāo)么化后，輸出為 （5-25） 為求 T 和 ，在圖 5-14 曲線上取兩點：E 和 F ，且 ，則有 由此解出前面的圖 5-13 則變?yōu)閳D 5-14。該標(biāo)么化處理，并不改當(dāng)然，為了計算上的方便，也可取 ， ，代入上兩

13、式，有 ， 算出 T 和 后，可檢驗一下用式（5-25）與實測曲線的誤差大小，如果誤差可接受，則所求的傳遞函數(shù)式（5-23）可用。否則，應(yīng)考慮用其它型傳遞函數(shù)（例如高階傳遞函數(shù)）來描述。 具體方法為，另取三點： 、 和 ，具體為， ， 由式（5-25）算得， ， 并分別與圖5-13中 、 、 對應(yīng)的縱坐標(biāo)比較即可。 （2）由階躍響應(yīng)曲線求二階慣性及其以上環(huán)節(jié)的參數(shù) 當(dāng)你用一階慣性環(huán)節(jié)近似被控過程傳遞函數(shù)，檢驗發(fā)現(xiàn)誤差不能滿足原定的精度時，可考慮二階及其以上的慣性環(huán)節(jié)傳遞函數(shù)。 當(dāng)然，為了計算上的方便，也可取 設(shè)有階躍響應(yīng)曲線如圖 5-15，現(xiàn)在，欲用二階慣性環(huán)節(jié)的傳遞函數(shù) （5-26） 來近似

14、描述它，其中， k 、 和 為待定的參數(shù)。 圖 5-15 階躍響應(yīng)曲線圖 5-16 具有純時延的階躍響應(yīng)曲線 當(dāng)輸入為 時，該傳遞函數(shù)的響應(yīng)為 （5-27） 設(shè)有階躍響應(yīng)曲線如圖 5-15，現(xiàn)在，欲 在圖 5-15 所示的階躍響應(yīng)曲線上，找出兩點：A( ， ) 和 B( ， ),并將這兩點分別代入式（5-27），有 其近似解為 研究表明，由式（5-27）表示的階躍響應(yīng)，應(yīng)有 ，并且當(dāng) 時，被控過程應(yīng)為一階慣性環(huán)節(jié) ，且時間常數(shù)為 當(dāng) 時，被控過程可為二階等容慣性環(huán)節(jié) ，且 在圖 5-15 所示的階躍響應(yīng)曲線上，找出兩當(dāng) 時，被控過程應(yīng)為二階以上慣性環(huán)節(jié)，可用等 n 容慣性環(huán)節(jié)來描述 其中 n

15、和 T 分別按下式計算 如果算得 n 不為整數(shù)，應(yīng)取最接近的整數(shù)。n 與 的關(guān)系也可用表 5-1 表示。 表5-1 多容過程的 n 與 之間的關(guān)系 n1234567891012140.3170.4600.5340.5840.6180.6400.6660.6840.6990.7120.7340.751當(dāng) 時，被控過程應(yīng)為二 如果階躍響應(yīng)曲線有明顯的純時延，如圖 5-16 所示，則應(yīng)在式（5-26）右邊乘上一個純時延環(huán)節(jié)： ，變?yōu)?其中 見圖5-16。具體用上面的公式求 、 時，應(yīng)在 和 中減去 時間段。 對于以上傳遞函數(shù)中的放大系數(shù) K，仍可用式（5-24）求取。 （3）由無自平衡過程的階躍響應(yīng)

16、曲線求過程參數(shù) 當(dāng)階躍響應(yīng)的曲線如圖 5-17 所示時，該過程的傳遞函數(shù)則具有無自平衡特性，其特點是，隨著 ，響應(yīng)曲線的變化速率逐漸趨于某一常數(shù)。 圖 5-17 無自平衡過程階躍響應(yīng)曲線 如果階躍響應(yīng)曲線有明顯的純時延，如圖 5- 根據(jù)該響應(yīng)曲線，該過程可用（5-28） 來近似。 當(dāng)階躍信號 作用于輸入端時，其輸出為 下面將討論 T 和 的確定。作階躍響應(yīng)直線部分的延長線（見圖中虛線段），與 t 軸相交于點 ，該線與時間軸 t 夾角為 ，于是 由于 ，所以 其中， 和 見圖 5-17。 根據(jù)該響應(yīng)曲線，該過程可用（5-28） 來近似。 由圖5-17可知，用式（5-28）近似原過程的最大誤差發(fā)生

17、在 這一段曲線上，即響應(yīng)的起始段。為此可再加一個慣性環(huán)節(jié)來減小誤差，即采用下列傳遞函數(shù)來描述過程 （5-29） 其中 T 的確定如上所示，即不變，而 和 的確定如下： 顯然，用式（5-29）描述響應(yīng)曲線表示的過程，比式（5-28）要精確些。 5.4 基于最小二乘法的過程辨識 應(yīng)該說，用最小二乘法建模仍然是一種測試法建模，但含有較多的處理技巧與方法，這里將其單列為一節(jié)，主要是考慮內(nèi)容稍多、篇幅較大。最小二乘法是系統(tǒng)辨識中的一種常用參數(shù)估計方法，它具有原理明了、算法簡捷、收斂較快、相對容易理解的特點，因而被廣泛用于參數(shù)估計之中。最小二乘法包括：最小二乘的批處理法、遞推法、漸消記憶法和增廣法等。 由

18、圖5-17可知，用式（5-28）近似原過程5.4.1 離散時間系統(tǒng)模型 數(shù)學(xué)模型分為連續(xù)時間系統(tǒng)模型和離散時間系統(tǒng)模型。前面討論的是連續(xù)時間系統(tǒng)模型，隨著計算機的普及與應(yīng)用，離散時間系統(tǒng)模型被越來越重視，最小二乘法采用的是離散時間系統(tǒng)模型。在這兩類模型中按是否有隨機擾動，每類又可分為確定性和隨機性兩種形式。 確定性離散系統(tǒng)（Deterministic Discrete Systems）的單輸入/輸出方程形式： （5-30） 式中 d 為純時延，且 ， 和 分別為k時刻的輸出和輸入。 隨機離散系統(tǒng)（Stochastic Discrete Systems）輸入-輸出差分模型一般有下面幾種： 自回歸

19、滑動平均（Auto-Regressive Moving Average，ARMA）模型 （5-31）式中 、 與式（5-30）中的相同， 5.4.1 離散時間系統(tǒng)模型 數(shù)學(xué)模型分為連 為白噪聲（White Noise）序列，且 ， 式（5-31）右邊第 1 項稱為滑動平均項，左邊項稱為自回歸項。式（5-31）也可寫為 這里， 分別被稱為過程模型和噪聲模型，后者也被稱為成形濾波器。而可以看作是白噪聲經(jīng)線性環(huán)節(jié)的輸出，它一般為有色噪聲（Coloured Noise）。 為白噪聲（White Noise）序列，且 自回歸積分滑動平均（Auto-Regressive Integrated Moving

20、 Average，ARIMA）模型或者這里， 、 、 和 與式（5-31）中相同。與式（5-31）相比，這里假定 。當(dāng) 時， 多項式中前 項的系數(shù)為零。 最小二乘模型（Least Square Model） 或者式中， 、 、d 和 有與前面相同的含義。 滑動平均（Moving Average，MA）模型 自回歸積分滑動平均（Auto-Reg與式（5-31）相比，這里有： 。 5.4.2 批處理最小二乘法 考慮最小二乘模型 式中， ， ， 為白噪聲。 設(shè)已知 、 ，現(xiàn)在的任務(wù)是根據(jù)可量測的輸入和輸出，確定參數(shù)： ， 由輸入-輸出模型有 （5-32） 令： 為觀測向量； 為待估參數(shù)向量，則式（5

21、-32）可寫為另一形式（5-33） 與式（5-31）相比，這里有： 現(xiàn)有 N 次觀測數(shù)據(jù)組并且 當(dāng) 時，由式（5-32）有： 引入下列符號， 現(xiàn)有 N 次觀測數(shù)據(jù)組并且 當(dāng) 由式（5-33），有下列矩陣形式： 由于真實的參數(shù)向量 并不知道，不妨用 來表示它的估計值，于是，基于 的輸出估計為 式中， 現(xiàn)在定義殘差 （也是一隨機變量）為實際輸出與估計輸出之差： （5-34） 對于 ，則有其中 。 現(xiàn)在的任務(wù)是：求使目標(biāo)函數(shù)為最小的 （記為 ）。 由式（5-33），有下列矩陣形式： 由于真實的參數(shù)向量 展開上式，有由求極值的方法，對 求一階導(dǎo)，并令其為零 從而有（5-35） 由于二階導(dǎo)所以由式（5-

22、35）求得 的為極小值。式（5-35）即為批處理法的最小二乘估計。 最小二乘估計的統(tǒng)計特性討論： 由于 為白噪聲序列，故有 ， 展開上式，有由求極值的方法，對 求 （1）無偏性：使 J 為最小的參數(shù)估計向量 的數(shù)學(xué)期望為參數(shù)真值向量，即 這是因為它揭示最小二乘參數(shù)估計是圍繞參數(shù)真值波動的統(tǒng)計性質(zhì)。 （2）估計誤差（偏差）協(xié)方差 主對角線上各元表現(xiàn)參數(shù)估計的散度，非對角線上各元反映參數(shù)估計 分量相互影響程度或相關(guān)性大小。 上式的成立，是因為 （1）無偏性：使 J 為最小的參數(shù)估計向量 （3）最小方差估計：設(shè) 為 的任一其他線性無偏估計，則 即最小二乘估計是最小方差估計，也就是參數(shù)估計 離參數(shù)真值

23、 最近。 由于 為 的任一線性無偏差估計，所以 可表示為 式中， 。由于， ，即 從而因為 （3）最小方差估計：設(shè) 為 所以 （4）一致收斂性：若 存在且正定，則 是一致收斂的，即 由于所以所以 （4）一致收斂性：若 存在且正定，則 由 的無偏性知 所以 例 5-1 現(xiàn)有最小二乘模型 其中 、 、 ， 為零均值白噪聲，實驗獲得輸入輸出數(shù)據(jù)如表 5-2，試用批處理最小二乘法確定多項式 和 的參數(shù)。 解：首先組成矩陣 ，然后按式（5-31）求出 。下面用Matlab 語言來做。 u = -1; -1; -1; -1; -1; 1; 1; -1; -1; 1; -1; 1; 1; -1; 1; 1;

24、 1; 1; -1; 1; -1; 1; -1; -1; -1; 1; -1; -1; 1; 1; 1 ; % 將u 輸入到工作空間 y = 0.0582; -0.6395; -1.8510; -2.0242; -1.4363; -0.9188; 0.3053; 2.2938; 1.0769; -2.2943; -1.9656; 0.4587; 1.3710; 1.8783; 0.2454; -1.1223; 0.7848; 2.4983; 2.2147; -0.2424; -1.5523; -0.5707; 0.5078; 0.7394; -1.4378; -2.6328; -0.5359;

25、 1.4520; -0.4325; -1.2545; 1.1510 ; % 將y 輸入到工作空間 由 的無偏性知 所以 表5-2 實驗獲得的輸入輸出數(shù)據(jù) tuytuytuy0.0005.5-1-1.965611.01-0.57070.5-10.05826.010.458711.5-10.50781.0-1-0.63956.511.371012.0-10.73941.5-1-1.85107.0-11.878312.5-1-1.43782.0-1-2.02427.510.245413.01-2.63282.5-1-1.43638.01-1.122313.5-1-0.53593.01-0.91888

26、.510.784814.0-11.45203.510.30539.012.498314.51-0.43254.0-12.29389.5-12.214715.01-1.25454.5-11.0769101-0.242415.511.15105.01-2.294310.5-1-1.5523-表5-2 實驗獲得的輸入輸出數(shù)據(jù) tuytuytuy0.00 phi = 0; -1*y(1:end-1), 0; 0; -1*y(1:end-2), 0; u(1:end-1), 0; 0; u(1:end-2); % 組建 Theta = inv (phi*phi) * phi*y；% 計算 運行結(jié)果為：

27、= -0.5076，0.6075，0.6854，0.7947 這與真值：-0.5, 0.6, 0.7, 0.8 相差不大。 另外，也可以用Matlab的辨識工具箱獲得結(jié)果。 T = iddata (y, u, 0.5); % 處理數(shù)據(jù) G = arx (T, 2,2,1); % 調(diào)用辨識函數(shù) 運行結(jié)果為： Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t)A(q) = 1 - 0.5077 q-1 + 0.6075 q-2 B(q) = 0.6863 q-1 + 0.7947 q-2Estimated using ARX from da

28、ta set T Loss function 0.00146694 and FPE 0.00190159 phi = 0;Sampling interval: 0.5其中 即為 ，符號 q 就是本書中的符號 z。由此可見，兩結(jié)果幾乎相同。 5.4.3 遞推最小二乘法 前面介紹的批處理法需要大量的計算機內(nèi)存（隨著N增大），而且當(dāng) 變化時， 不能自動跟蹤其變化，實時性不好。而遞推最小二乘法（Recursive Least Square Method）則能解決這一問題。 引理：設(shè)A、C 和 均為非奇異方陣，則有 證明：用 左乘上式右端，若結(jié)果為單位矩陣，則說明其等式關(guān)系成立。 當(dāng) 時，有 Sampl

29、ing interval: 0.5其中 （5-36） 遞推最小二乘算法 最小二乘模型描述的系統(tǒng) 基于 及以前的輸出 y 、 及以前的輸入 u ，未知參數(shù)向量 的最小二乘估計 的遞推公式為 （5-37） （5-38） （5-39） （5-36） 遞推最小二乘算法 式中， 為觀測向量， 為未知參數(shù)向量 的估計值， 為增益向量。 證明： 設(shè) 是基于 和 的估計，根據(jù)批處理最小二乘法有 對于 時的最小二乘估計為 式中 ， 于是式中， 為觀測向量， 為未知參數(shù)向量 的估計值， （5-40） 令并視 、 和 分別為式（5-36）中的 A、B、D，再考慮 在此為標(biāo)量，由引理有 再令 ，則有 又令（5-41）

30、 （5-40） 令并視 、 即為式（5-39），且式（5-41）可寫為它就是式（5-39）。 將式（5-39）代入式（5-40），并注意交換標(biāo)量和矩陣的位置，有即為式（5-39），且式（5-41）可寫為它就是式（5-39即為式（5-37）。 遞推最小二乘算法的說明： （1）從式（5-37）看，新的參數(shù)向量估計值 為先前的參數(shù)向量估計值 加修正項 在不斷更新過程中， 、 和 的行列數(shù)不變，但它們的舊數(shù)據(jù)不斷被新數(shù)據(jù)替換； （2） 為增益向量， 為誤差的協(xié)方差陣，一般 與 成正比，協(xié)方差越大，說明估計值與真值相差越大，增益向量也會越大，所產(chǎn)生的校正作用也越大； （3）初值 和 的確定。 方法1：若

31、已有 組數(shù)據(jù)，則可批處理它們，并將結(jié)果作為初值，即 ， 方法2： ， ，其中 。 5.4.4 具有遺忘因子的遞推最小二乘法 遞推最小二乘法有一個缺點：常常出現(xiàn)“數(shù)據(jù)飽和”。隨著k的增加， 和 變得越來越小，式（5-37）中的修正項對 的修正能力變得越來越弱，即新近加入的輸入/輸出數(shù)據(jù)對參數(shù)向量估計值的更新作用不大。這樣導(dǎo)致的結(jié)果是：參數(shù)估計值難以接近真值；當(dāng)參數(shù)真值時變時，該算法無法跟蹤這種變化，從而使實時參數(shù)辨識失敗。 即為式（5-37）。 遞推最小二乘算法的說明： 解決該問題的方法之一是用具有遺忘因子（Forgetting Factor）的遞推最小二乘法。 取性能指標(biāo)函數(shù)： 式中， 為加權(quán)

32、對角陣： N為觀測數(shù)據(jù)組數(shù)， 為遺忘因子： 。 按與前面相同的思路，可推出具有遺忘因子的遞推最小二乘估計公式： 解決該問題的方法之一是用具有遺忘因子（Fo式中 。 幾點說明： （1）當(dāng) 時，該估計公式組即為遞推最小二乘算法公式組； （2） 的選取范圍一般在： ，參數(shù)變化快時， 取小點；變化慢時，取大些； （3）初值的選取與前面的遞推最小二乘法相同。 5.4.5 遞推增廣最小二乘法 以上用的是最小二乘模型，在很多情況下， ，所以需考慮更一般的情況，即ARMA模型： （5-42） 其中， , , 式（5-42）可變?yōu)?式中 將其表示為 式中在 中，由于 不可測，所以只能用其估計值 來替換，設(shè) 式中

33、 用與前面類似的方法，并考慮用 代替 中的 ，則有下列遞推公式 將其表示為 式中在 中，由于 不可測 的估計可用下列方法之一來進(jìn)行： （1）預(yù)測形式： （2）濾波形式：并且當(dāng) 時，有 。 前面介紹了最小二乘參數(shù)估計算法和特性，但是估計參數(shù)是否收斂到真實值尚未提及，由于這方面的內(nèi)容涉及較多的定理及證明，這里只能給出結(jié)論，證明可參考有關(guān)資料。 被估參數(shù)收斂定理：若 是 N 階持續(xù)激勵（Persistently Exciting of Order N），即充分豐富（Sufficient Rich），則式（5-33）式（5-35）所示的算法能保證指數(shù)收斂。 的估計可用下列方法之一來進(jìn)行： 例 5.2

34、現(xiàn)有自回歸滑動平均模型 已知 ， 為白噪聲：零均值，方差為 0.1，取采樣周期 ，實驗獲得輸入 / 輸出數(shù)據(jù)如表 5-3，試用遞推增廣最小二乘法確定對象參數(shù)： 、 、 、 、 和 ，并將參數(shù)估計變化過程用圖形表示出來。 表5-3 輸入/輸出數(shù)據(jù)u1.00001.00000.00000.0000-1.0000-1.00001.00001.0000y0.11650.08021.02512.34462.13011.5708-0.1742-1.6821u0.00000.00001.00001.00000.00000.0000-1.0000-1.0000y-0.90941.11281.32800.928

35、11.69902.56812.41131.4639u1.00001.00000.00000.00001.00001.00000.00000.0000y-0.1749-1.9893-1.10591.08121.44970.93041.50192.6898u-1.0000-1.00001.00001.00000.00000.00001.00001.0000y2.35731.5460-0.1111-1.7658-0.90110.89141.47701.1226u0.00000.0000-1.0000-1.00001.00001.00000.00000.0000y1.68002.40182.30391

36、.4445-0.1399-1.9036-0.96571.1687 例 5.2 現(xiàn)有自回歸滑動平均模型 已 解：由遞推增廣最小二乘法公式，有參數(shù)估計算法： (1) 初始化 、 等，并給 u 和 y 賦值； (2) ； (3) 計算 ； (4) 計算 ； (5) 保存參數(shù) 、 、 、 、 和 ； (6) 計算 ； (7) 求濾波形式的 ； (8) 移位處理 u、y、 和 ； (9) k 步終了嗎? 否， ，去第 2 步；是，去第 10 步； (10) 打印結(jié)果。 按算法編寫程序，經(jīng)調(diào)試，得出 、 、 、 、 和 隨拍數(shù)變化的規(guī)律，具體參見圖5-18。具體程序為（用 Matlab 編寫）： clea

37、rN=40; 解：由遞推增廣最小二乘法公式，有參數(shù)估計算 u=1,1,0,0,-1,-1,1,1,0,0,1,1,0,0,-1,-1,1,1,0,0,1,1,0,0,-1,-1,1,1,0,0,1,1,. 0,0,-1,-1,1,1,0,0; %錄入輸入數(shù)據(jù) y=0.1165,0.0802,1.0251,2.3446,2.1301,1.5708,-0.1742,-1.6821,-0.9094, 1.1128,1.3280,0.9281,1.6990,2.5681,2.4113,1.4639,-0.1749,-1.9893, -1.1059,1.0812,1.4497,0.9304,1.5019

38、,2.6898,2.3573,1.5460,-0.1111, -1.7658,-0.9011, 0.8914,1.4770,1.1226,1.6800,2.4018,2.3039,1.4445, -0.1399,-1.9036,-0.9657,1.1687; % 錄入輸出數(shù)據(jù) I=eye(6); P=1e9*I; Q0=zeros(6,1); h0=1;t=0; %初始化 y1=0; y2=0; u1=0; u2=0; u3=0; x1=0; x2=0; for j=1:N f=-y1;-y2;u2;u3;x1;x2; %構(gòu)成 K=P*f/(1+f*P*f); Q=Q0+K*(y(j)-f*Q0); % a1(j)=Q(1);a2(j)=Q(2);b0(j)=Q(3);b1(j)=Q(4);c1(j)=Q(5);c2(j)=Q(6); P=(eye(6)-K*f)*P; x(j)=y(j)-f*Q; %求 u=1,1,0,0,-1,-1,1, y2=y1;y1=y(j);u3=u2;u2=u1;u1=u(j);x2=x1