線性擬合方法

線性擬合方法第一頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第五章實(shí)驗(yàn)數(shù)據(jù)及模型參數(shù)擬合方法第一節(jié)問題的提出第二頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第一節(jié)問題的提出在化工設(shè)計(jì)及化工模擬計(jì)算中，需要大量的物性參數(shù)及各種設(shè)備參數(shù)。這些參數(shù)有些可以通過計(jì)算得到，但大量的參數(shù)還是要通過實(shí)驗(yàn)測(cè)量得到。實(shí)驗(yàn)測(cè)量得到的常常是一組離散數(shù)據(jù)序列（xi,yi）。

如果數(shù)據(jù)序列（xi,yi）(為一般起見),i=1,2,…,m,含有不可避免的誤差（或稱“噪聲”，如圖5-1所示），或者無(wú)法同時(shí)滿足某特定的函數(shù)（如圖5-2所示）,那么，只能要求所作逼近函數(shù)ψ（x）最優(yōu)地靠近樣點(diǎn)，即向量Q=（ψ(x1),ψ(x2),…,ψ(xm)）T與Y=(y1,y2,…,ym)T的誤差或距離最小。按Q與Y之間誤差最小原則作為“最優(yōu)”標(biāo)準(zhǔn)構(gòu)造的逼近函數(shù)，稱為擬合函數(shù)。第三頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第一節(jié)問題的提出圖5-1含有噪聲的數(shù)據(jù)圖5-2無(wú)法同時(shí)滿足某特定函數(shù)的數(shù)據(jù)序列第四頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第一節(jié)問題的提出除了物性數(shù)據(jù)及設(shè)備參數(shù)需要利用數(shù)據(jù)擬合外，在化學(xué)化工中，許多模型也要利用數(shù)據(jù)擬合技術(shù)，求出最佳的模型和模型參數(shù)。如在某一反應(yīng)工程實(shí)驗(yàn)中，我們測(cè)得了如表5-1所示的實(shí)驗(yàn)數(shù)據(jù)。

現(xiàn)在要確定在其他條件不變的情況下，轉(zhuǎn)化率y和溫度T的具體關(guān)系，現(xiàn)擬用兩種模型去擬合實(shí)驗(yàn)數(shù)據(jù)，兩種模型分別是

第五頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第一節(jié)問題的提出如何求取上述模型中的參數(shù)，并判斷兩種模型的優(yōu)劣,是化學(xué)化工工作者經(jīng)常要碰到的問題，這個(gè)問題的求解將在本章下面的有關(guān)章節(jié)中進(jìn)行詳細(xì)的講解。

第六頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第二節(jié)擬合的標(biāo)準(zhǔn)第七頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第二節(jié)擬合的標(biāo)準(zhǔn)

前面已經(jīng)提到按Q與Y之間誤差最小原則作為“最優(yōu)”標(biāo)準(zhǔn)構(gòu)造的逼近函數(shù)，稱為擬合函數(shù)，而向量Q與Y之間的誤差或距離有各種不同的定義方法，一般有以下幾種。

（1）用各點(diǎn)誤差絕對(duì)值的和表示

（2）用各點(diǎn)誤差按絕對(duì)值的最大值表示

（3）用各點(diǎn)誤差的平方和表示

第八頁(yè)，共五十四頁(yè)，編輯于2023年，星期一式中R稱為均方誤差。由于計(jì)算均方誤差的最小值的原則容易實(shí)現(xiàn)而被廣泛采用。按均方誤差達(dá)到極小構(gòu)造擬合曲線的方法稱為最小二乘法。同時(shí)還有許多種其他的方法構(gòu)造擬合曲線，感興趣的讀者可參閱有關(guān)教材。本章主要講述用最小二乘法構(gòu)造擬合曲線。第二節(jié)擬合的標(biāo)準(zhǔn)

第九頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第二節(jié)擬合的標(biāo)準(zhǔn)__實(shí)例1

實(shí)驗(yàn)測(cè)得二甲醇（DME）的飽和蒸氣壓和溫度的關(guān)系，見表5-2。

序號(hào)溫度℃蒸氣壓MPa1-23.70.1012-100.174300.2544100.3595200.4956300.6627400.880表5-2DME飽和蒸氣壓和溫度的關(guān)系由表5-2的數(shù)據(jù)觀測(cè)可得，DME的飽和蒸氣壓和溫度有正相關(guān)關(guān)系，如果以函數(shù)p=a+bt來擬合,則擬合函數(shù)是一條直線。通過計(jì)算均方誤差Q(a,b)最小值而確定直線方程。（見圖5-3）圖5-3DME飽和蒸汽壓和溫度之間的線性擬合

第十頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第二節(jié)擬合的標(biāo)準(zhǔn)__實(shí)例1

擬合得到直線方程為：

相關(guān)系數(shù)R為0.97296，平均絕對(duì)偏差SD為0.0707。

第十一頁(yè)，共五十四頁(yè)，編輯于2023年，星期一擬合的標(biāo)準(zhǔn)

——

實(shí)例2如果采用二次擬合，通過計(jì)算下述均方誤差

擬合得二次方程為

相關(guān)系數(shù)R為0.99972，平均絕對(duì)偏差SD為0.00815,具體擬合曲線見圖5-4。圖5-4DME飽和蒸氣壓和溫度之間的二次擬合

第十二頁(yè)，共五十四頁(yè)，編輯于2023年，星期一擬合的標(biāo)準(zhǔn)實(shí)例2比較圖5-3和圖5-4以及各自的相關(guān)系數(shù)和平均絕對(duì)偏差可知：對(duì)于DME飽和蒸汽壓和溫度之間的關(guān)系，在實(shí)驗(yàn)溫度范圍內(nèi)用二次擬合曲線優(yōu)于線性擬合。二次擬合曲線具有局限性，由圖5-4觀察可知，當(dāng)溫度低于-30℃時(shí)，飽和壓力有升高的趨勢(shì)，但在擬合的溫度范圍內(nèi)，二次擬合的平均絕對(duì)偏差又小于一次擬合，故對(duì)物性數(shù)據(jù)進(jìn)行擬合時(shí)，不僅要看在擬合條件下的擬合效果，還必須根據(jù)物性的具體性質(zhì)，判斷在擬合條件之外的物性變化趨勢(shì)，以便使擬合公式在已做實(shí)驗(yàn)點(diǎn)數(shù)據(jù)之外應(yīng)用。第十三頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第三節(jié)單變量擬合和多變量擬合

給定一組數(shù)據(jù)（xi,yi）,i=1,2,…,m，作擬合直線p(x)=a+bx,均方誤差為

由數(shù)學(xué)知識(shí)可知，Q(a,b)的極小值需滿足：

整理得到擬合曲線滿足的方程：

5.3.1單變量擬合1.單變量擬合——

線性擬合第十四頁(yè)，共五十四頁(yè)，編輯于2023年，星期一1.單變量擬合——線性擬合該方程可用消元法或克萊姆方法解出方程（如下式所示）

第十五頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第十六頁(yè)，共五十四頁(yè)，編輯于2023年，星期一單變量擬合——

線性擬合實(shí)例例：下表為實(shí)驗(yàn)測(cè)得的某一物性和溫度之間的關(guān)系數(shù)據(jù)，表中x為溫度數(shù)據(jù)，y為物性數(shù)據(jù)。請(qǐng)用線性函數(shù)擬合溫度和物性之間的關(guān)系。

解：設(shè)擬合直線p(x)=a+bx,并計(jì)算得下表：

x7911131517192123252729y91215182124273033363942x313335373941434547

y454851545760636669

編號(hào)xyxyx212345…21Σ79111315…47567912151821…6981963108165234315…3243267334981121169225…220918389第十七頁(yè)，共五十四頁(yè)，編輯于2023年，星期一將數(shù)據(jù)代入法方程組（1-12）中，得到：解方程得：a=-1.5,b=1.5

擬合直線為：

相關(guān)系數(shù)R為1。單變量擬合——

線性擬合實(shí)例第十八頁(yè)，共五十四頁(yè)，編輯于2023年，星期一線性擬合VB清單PrivateSubCommand1_Click()Dimx(5),y(5),c,d,m,p,a,b,eerConstn=5Fori=1To5x(i)=InputBox(“x(i)=”)y(i)=InputBox(“y(i)=”)Print“x(i)=”;x(i)Print“y(i)=”;y(i)NextIc=0d=0m=0p=0

Fori=1To5c=c+x(i)d=d+x(i)^2m=m+y(i)p=p+x(i)*y(i)Nextia=(m*d-c*p)/(n*d-c^2)b=(n*p-c*m)/(n*d-c^2)‘參數(shù)計(jì)算

a=Int(a*1000+0.5)/1000b=Int(b*1000+0.5)/1000Text1.Text=Str(a)Text2.Text=Str(b)‘參數(shù)輸出Fori=1To5eer=eer+(a+b*x(i)-y(i))^2‘誤差計(jì)算

eer=Int(eer*100000+0.5)/100000Nextieer=eer/5Text3.Text=Str(eer)EndSub第十九頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第二十頁(yè)，共五十四頁(yè)，編輯于2023年，星期一有關(guān)線性擬合變型問題例如要擬合y=a+b/x2，只需在數(shù)據(jù)輸入后增加一語(yǔ)句x(i)=1/x(i)^2，而在程序后面的誤差eer的計(jì)算中則不需要修改。第二十一頁(yè)，共五十四頁(yè)，編輯于2023年，星期一2.單變量擬合二次擬合函數(shù)給定數(shù)據(jù)（xi,yi),i=1,2,…,m,用二次多項(xiàng)式函數(shù)擬合這組數(shù)據(jù)。

設(shè),作出擬合函數(shù)與數(shù)據(jù)序列的均方誤差表達(dá)式由數(shù)學(xué)知識(shí)可知，Q(a0,a1

,a2)的極小值滿足：整理左式得二次多項(xiàng)式函數(shù)擬合的滿足條件方程（5-14）：

(5-14)第二十二頁(yè)，共五十四頁(yè)，編輯于2023年，星期一解此方程得到在均方誤差最小意義下的擬合函數(shù)p(x)。式（5-14）稱為多項(xiàng)式擬合的法方程，法方程的系數(shù)矩陣是對(duì)稱的。當(dāng)擬合多項(xiàng)式n>5時(shí)，法方程的系數(shù)矩陣是病態(tài)的，在用通常的迭代方法求解線性方程時(shí)會(huì)發(fā)散，在計(jì)算中要采用一些特殊算法以保護(hù)解的準(zhǔn)確性。關(guān)于線性方程的求解方法，已在第三章中介紹。2.單變量擬合二次擬合函數(shù)第二十三頁(yè)，共五十四頁(yè)，編輯于2023年，星期一3.二次擬合函數(shù)的拓展和一次擬合一樣，二次擬合也可以有多種變型，例如套用上面的公式，可以得到關(guān)于求解此擬合函數(shù)的法方程（5-15）。值得注意的是在此法方程的構(gòu)建過程中，進(jìn)行了變量的代換。首先是擬合函數(shù)中變量的代換:。

P(x)=a0+a1x3+a2x5

（5-15）第二十四頁(yè)，共五十四頁(yè)，編輯于2023年，星期一其次是法方程的代換：將相應(yīng)擬合函數(shù)中的代換引入法方程中。同時(shí)應(yīng)注意法方程中x的4次冪是由兩個(gè)2次冪相乘得到，x的3次冪是由一個(gè)2次冪和一個(gè)1次冪相乘得到，而2次冪就是變量本身，而非兩個(gè)1次冪相乘得到。這個(gè)概念至關(guān)重要，在以后的二次擬合的各類變型中，均需利用這個(gè)概念，千萬(wàn)不要用常規(guī)的思路去進(jìn)行代入計(jì)算。

3.二次擬合函數(shù)的拓展如果我們需要求解是下面的擬合函數(shù)：

第二十五頁(yè)，共五十四頁(yè)，編輯于2023年，星期一參照上面的方法，我們很容易得到求解該擬合函數(shù)的法方程3.二次擬合函數(shù)的拓展第二十六頁(yè)，共五十四頁(yè)，編輯于2023年，星期一4.二次擬合實(shí)例請(qǐng)用二次多項(xiàng)式函數(shù)擬合下面這組數(shù)據(jù)。解：設(shè)

并計(jì)算得下表序號(hào)xyxyX2X2yX3X41-34-12936-27812-22-448-8163-13-313-114000000051-1-11-11162-2-44-881673-5-159-452781Σ01-3928-70196第二十七頁(yè)，共五十四頁(yè)，編輯于2023年，星期一4.二次擬合實(shí)例將上面數(shù)據(jù)代入式

(5-14),相應(yīng)的法方程為解方程得：

a0=0.66667,a1=-1.39286,a2=-0.13095所以：-3-2-10123-6-4-2024y=0.66667-1.39286x-0.13095x2yYX圖

5-6擬合曲線與數(shù)據(jù)序列

第二十八頁(yè)，共五十四頁(yè)，編輯于2023年，星期一二次擬合--VB程序清單PrivateSubCommand1_Click()Dimn,mAsIntegern=InputBox("n=方程次數(shù)")m=InputBox("m=實(shí)驗(yàn)次數(shù)")ReDimx(m),y(m),a0(n+1),a1(n+1),aa(m,n+1)ReDimqq(n+1,m),pp(n+1,n+1),b(n+1),g(n+1),tt(n+1,n+1)'ForI=1Tom:'讀入數(shù)據(jù)'READx(I),y(I)'NextI'Data-3,4,-2,2,-1,3,0,0,1,-1,2,-2,3,-5Fori=1Tomx(i)=InputBox("x(i)")y(i)=InputBox("y(i)")Nexti第二十九頁(yè)，共五十四頁(yè)，編輯于2023年，星期一二次擬合--VB程序清單omiga=InputBox("omiga=松弛因子")Forj=1Ton+1'為計(jì)算法方程中的系數(shù)做準(zhǔn)備

Fori=1Tomaa(i,j)=x(i)^(j-1)NextiNextjForj=1TomFori=1Ton+1qq(i,j)=aa(j,i)NextiNextj'計(jì)算法方程中的右邊項(xiàng)Fori=1Ton+1b(i)=0Forj=1Tomb(i)=b(i)+aa(j,i)*y(j)NextjNexti'開始計(jì)算法方程中的右邊項(xiàng)的系數(shù)Fori=1Ton+1Forj=1Ton+1pp(i,j)=0Fork=1Tompp(i,j)=qq(i,k)*aa(k,j)+pp(i,j)NextkNextjNexti第三十頁(yè)，共五十四頁(yè)，編輯于2023年，星期一二次擬合--VB程序清單'開始用松弛迭代法求解法方程中的變量Fori=1Ton+1a0(i)=0a1(i)=0.11NextiFori=1Ton+1g(i)=b(i)/pp(i,i)Forj=1Ton+1Ifj=iThentt(i,j)=0Elsett(i,j)=-pp(i,j)/pp(i,i)EndIfNextjNexti50eer=0Fori=1Ton+1eer=eer+Abs(a0(i)-a1(i))NextiIfeer<0.00001ThenGoTo100ElseFori=1Ton+1a0(i)=a1(i)Nexti第三十一頁(yè)，共五十四頁(yè)，編輯于2023年，星期一二次擬合--VB程序清單'關(guān)鍵的迭代計(jì)算公式Fori=1Ton+1s=g(i)Forj=1Ton+1s=s+tt(i,j)*a0(j)Nextja1(i)=omiga*s+(1-omiga)*a0(i)'Print"a1("&i&")="&a1(i)NextiGoTo50EndIf'打打印結(jié)果100Fori=1Ton+1Print"a("&i-1&")="&a1(i)NextiEndSub第三十二頁(yè)，共五十四頁(yè)，編輯于2023年，星期一3.2多變量的曲線擬合

前面介紹的曲線擬合方法只涉及單變量函數(shù)的曲線擬合，但實(shí)際在化工實(shí)驗(yàn)數(shù)據(jù)處理及模型參數(shù)擬合時(shí)，通常會(huì)碰到多變量的參數(shù)擬合問題。一個(gè)典型的例子是傳熱實(shí)驗(yàn)中努塞爾數(shù)、雷諾數(shù)及普朗特?cái)?shù)之間的擬合問題：

(1-16)

根據(jù)若干組實(shí)驗(yàn)測(cè)得的數(shù)據(jù)，如何求出式（1-16）中的參數(shù)c1、c2、c3，這是一個(gè)有2個(gè)變量的參數(shù)擬合問題，為不失一般性，我們把它表達(dá)成以下形式。給定數(shù)據(jù)序列用一次多項(xiàng)式函數(shù)擬合這組數(shù)據(jù)。

設(shè)，作出擬合函數(shù)與數(shù)據(jù)序列的均方誤差

(1-17)

(x1i,x2i,yi),i=1,2,3…,m

第三十三頁(yè)，共五十四頁(yè)，編輯于2023年，星期一多變量的曲線擬合由多元函數(shù)的極值原理，Q(a0

,a1,a2

)的極小值滿足

整理得多變量一次多項(xiàng)式函數(shù)擬合的法方程

(1-18)第三十四頁(yè)，共五十四頁(yè)，編輯于2023年，星期一通過求解方程（1-18）就可以得到多變量函數(shù)線性擬合時(shí)的參數(shù)，由于方程（1-16）不是線性方程，我們可以通過對(duì)方程（1-16）兩邊同取對(duì)數(shù)，就可以得到以下線性方程只要作如下變量代換：并將實(shí)驗(yàn)數(shù)據(jù)代入法方程（1-18）就可以求出方程（1-16）中的系數(shù)。對(duì)于變量數(shù)多于2個(gè)，并且擬合曲線模型是非線性型時(shí)，可參照本節(jié)的方法，推導(dǎo)得到法方程，通過對(duì)法方程的求解就可以求得各種擬合曲線參數(shù)。靈活運(yùn)用上面介紹的方法，可以解決大部分實(shí)驗(yàn)數(shù)據(jù)及模型參數(shù)的擬合問題。多變量的曲線擬合第三十五頁(yè)，共五十四頁(yè)，編輯于2023年，星期一多變量的曲線擬合—VB程序清單第三十六頁(yè)，共五十四頁(yè)，編輯于2023年，星期一多變量的曲線擬合實(shí)例根據(jù)某傳熱實(shí)驗(yàn)測(cè)得如下數(shù)據(jù)，請(qǐng)用方程（1-16）的形式擬合實(shí)驗(yàn)曲線。

解：利用上面的VB程序，將數(shù)據(jù)依次輸入，就可以得到方程（1-16）中的三個(gè)參數(shù)

C1=0.023C2=0.8C3=0.3

則式（1-16）就變成了常見的光滑管傳熱方程第三十七頁(yè)，共五十四頁(yè)，編輯于2023年，星期一多變量的曲線擬合實(shí)例值得注意的是程序中對(duì)c2(1)的處理，不是直接將計(jì)算結(jié)果顯示出來，而是進(jìn)行指數(shù)運(yùn)算后才顯示出來。這是由于我們?cè)谶M(jìn)行擬合計(jì)算的時(shí)候，對(duì)方程（1-16）進(jìn)行了對(duì)數(shù)運(yùn)算。如果擬合方程的形式和方程（1-16）不同，則需對(duì)上面提供的程序作適當(dāng)修改。例如以下兩個(gè)自變量的擬合函數(shù)第三十八頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第四節(jié)解矛盾方程組第三十九頁(yè)，共五十四頁(yè)，編輯于2023年，星期一第四節(jié)解矛盾方程組本節(jié)中將用最小二乘法求解線性矛盾方程的方法來構(gòu)造擬合函數(shù)，并將其推廣至任意次和任意多個(gè)變量的擬合函數(shù)，為在化學(xué)化工中實(shí)驗(yàn)數(shù)據(jù)處理及模型參數(shù)擬合提供更為一般性的方法。給定數(shù)據(jù)序列（xi,yi）,i=1,2,…,m

，做擬合直線p(x)=a0

+a1x，如果要直線p(x)過這些點(diǎn)，那么就有

p(xi)=a0+a1xi

=yi,i=1,2,…,m,即寫成矩陣形式為第四十頁(yè)，共五十四頁(yè)，編輯于2023年，星期一上述方程組中有2個(gè)未知量m個(gè)方程(m>>2)。一般地，將含有n個(gè)未知量m個(gè)方程的線性方程組其一般形式為第四節(jié)解矛盾方程組寫成矩陣形為

第四十一頁(yè)，共五十四頁(yè)，編輯于2023年，星期一一般情況下，當(dāng)方程數(shù)m多于變量數(shù)n

，且m個(gè)方程之間線性無(wú)關(guān),則方程組無(wú)解，這時(shí)方程組稱為矛盾方程組。方程組在一般意義下無(wú)解，也即無(wú)法找到n個(gè)變量同時(shí)滿足m個(gè)方程。這種情況和擬合曲線無(wú)法同時(shí)滿足所有的實(shí)驗(yàn)數(shù)據(jù)點(diǎn)相仿，故可以通過求解均方誤差極小意義下矛盾方程的解來獲取擬合曲線。由數(shù)學(xué)知識(shí)還可證明：方程組ATAX=ATY的解就是矛盾方程組AX=Y

在最小二乘法意義下的解。這樣我們只要通過求解ATAX=ATY就可以得到矛盾方程組的解，進(jìn)而得到各種擬合曲線，為擬合曲線的求解提供了另一種方法。第四節(jié)解矛盾方程組min║AX-B║22第四十二頁(yè)，共五十四頁(yè)，編輯于2023年，星期一例如，擬合直線p(x)=a0+a1x的矛盾方程組ATAX=ATY的形式如下：

化簡(jiǎn)得到與式(1-12)相同的法方程第四節(jié)解矛盾方程組第四十三頁(yè)，共五十四頁(yè)，編輯于2023年，星期一這里需要注意的是變量X和系數(shù)（a0，

a1）之間的相互轉(zhuǎn)換關(guān)系。即對(duì)于n次多項(xiàng)式曲線擬合，要計(jì)算Q(a0,a1,…,an)

的極小值問題,這與解矛盾方程組第四節(jié)解矛盾方程組與求

的極小問題是一回事?；虻谒氖捻?yè)，共五十四頁(yè)，編輯于2023年，星期一在這里故對(duì)離散數(shù)據(jù)（x