线性拟合方法

《线性拟合方法》由会员分享，可在线阅读，更多相关《线性拟合方法（54页珍藏版）》请在装配图网上搜索。

1、第五章第五章 实验数据及模型参数拟合方法实验数据及模型参数拟合方法第五章第五章 实验数据及模型参数拟合方法实验数据及模型参数拟合方法第一节第一节 问题的提出问题的提出第一节第一节 问题的提出问题的提出 在化工设计及化工模拟计算中，需要大量的物性参数及各种设备参数。这些参数有些可以通过计算得到，但大量的参数还是要通过实验测量得到。实验测量得到的常常是一组离散数据序列（xi ,yi）。 如果数据序列（xi ,yi）(为一般起见), i=1,2, ,m ,含有不可避免的误差（或称“噪声” ，如图5-1所示），或者无法同时满足某特定的函数（如图5-2所示）,那么，只能要求所作逼近函数（x）最优地靠近样

2、点，即向量Q=（(x1), (x2), , (xm)）T与Y=(y1,y2, ,ym)T的误差或距离最小。按Q与Y之间误差最小原则作为“最优”标准构造的逼近函数，称为拟合函数。第一节第一节 问题的提出问题的提出-202468101214161820050100150200Y X 024681005101520Y X 图5-1 含有噪声的数据图5-2 无法同时满足某特定函数的数据序列第一节第一节 问题的提出问题的提出 除了物性数据及设备参数需要利用数据拟合外，在化学化工中，许多模型也要利用数据拟合技术，求出最佳的模型和模型参数。如在某一反应工程实验中，我们测得了如表5-1所示的实验数据。 序号

3、1 2 3 4 5 6 7 8 温度 T 10 20 30 40 50 60 70 80 转化率 y 0.1 0.3 0.7 0.94 0.95 0.68 0.34 0.13 现在要确定在其他条件不变的情况下，转化率y和温度T的具体关系，现拟用两种模型去拟合实验数据，两种模型分别是 2111TcTbay2222)45(Tbacy第一节第一节 问题的提出问题的提出 如何求取上述模型中的参数，并判断两种模型的优劣,是化学化工工作者经常要碰到的问题，这个问题的求解将在本章下面的有关章节中进行详细的讲解。 2111TcTbay2222)45(Tbacy第二节第二节 拟合的标准拟合的标准第二节第二节 拟

4、合的标准拟合的标准 前面已经提到按Q与Y之间误差最小原则作为“最优”标准构造的逼近函数，称为拟合函数，而向量Q与Y之间的误差或距离有各种不同的定义方法，一般有以下几种。 （1）用各点误差绝对值的和表示 （2）用各点误差按绝对值的最大值表示 （3）用各点误差的平方和表示 miiiyxR11)(iimiyxR)(max122212Y-Q(x)R )(或imiiyxRR 式中R称为均方误差。由于计算均方误差的最小值的原则容易实现而被广泛采用。按均方误差达到极小构造拟合曲按均方误差达到极小构造拟合曲线的方法称为最小二乘法线的方法称为最小二乘法。同时还有许多种其他的方法构造拟合曲线，感兴趣的读者可参阅有

5、关教材。本章主要讲述用最小二乘法构造拟合曲线。第二节第二节 拟合的标准拟合的标准 第二节第二节 拟合的标准拟合的标准_实例实例1 1 实验测得二甲醇（DME）的饱和蒸气压和温度的关系，见表5-2。 序号温度 蒸气压 MPa1-23.70.1012-100.174300.2544100.3595200.4956300.6627400.880表5-2 DME饱和蒸气压和温度的关系-30-20-10010203040500.00.20.40.60.81.0pt 由表5-2的数据观测可得，DME的饱和蒸气压和温度有正相关关系，如果以函数p=a+bt来拟合,则拟合函数是一条直线。通过计算均方误差Q (

6、a , b )最小值而确定直线方程。 （见图5-3 ）图5-3 DME饱和蒸汽压和温度之间的线性拟合 第二节第二节 拟合的标准拟合的标准_实例实例1 1 2121)()(),(imiiimiipbtaptpbaQ 拟合得到直线方程为： 相关系数R为0.97296，平均绝对偏差SD为0.0707。 tp0.01210.30324拟合的标准拟合的标准 实例实例2 2 如果采用二次拟合，通过计算下述均方误差 拟合得二次方程为 相关系数R为0.99972，平均绝对偏差SD为0.00815,具体拟合曲线见图5-4。21221012210)()(),(imiiimiiiptataaptpaaaQ20001

7、50009570248450t.t.p-30-20-10010203040500.00.20.40.60.81.0y=0.24845+0.00957 x+0.00015 x2压力, P(MPa)温 度, t()图5-4 DME饱和蒸气压和温度之间的二次拟合 拟合的标准拟合的标准 实例实例2 2 比较图5-3和图5-4以及各自的相关系数和平均绝对偏差可知：n对于DME饱和蒸汽压和温度之间的关系，在实验温度范围内用二次拟合曲线优于线性拟合。n二次拟合曲线具有局限性，由图5-4观察可知，当温度低于-30时，饱和压力有升高的趋势，但在拟合的温度范围内，二次拟合的平均绝对偏差又小于一次拟合，故对物性数据

8、进行拟合时，不仅要看在拟合条件下的拟合效果，还必须根据物性的具体性质，判断在拟合条件之外的物性变化趋势，以便使拟合公式在已做实验点数据之外应用。 第三节第三节 单变量拟合和多变量拟合单变量拟合和多变量拟合 给定一组数据（xi,yi）,i=1, 2 , , m ， 作拟合直线p (x)=a + bx , 均方误差为 2121)()(),(imiiimiiybxayxpbaQ由数学知识可知，Q (a , b)的极小值需满足： 0)(2),(1imiiybxaabaQ0)(2),(1iimiixybxabbaQ整理得到拟合曲线满足的方程： mimimiiiiimimiiiyxbxaxybxma111

9、211)()()(5.3.1 单变量拟合单变量拟合 1. 单变量拟合单变量拟合 线性拟合线性拟合1.单变量拟合单变量拟合线性拟合线性拟合 该方程可用消元法或克莱姆方法解出方程（如下式所示） )()()(/()( 2112111211211121121112111mimiiimimiimiiiimiimiimiiimiimiimiimiimiimiimiimiiimiimiixxmyxyxmbxxmyxxxyxxxmxyxxya单变量拟合单变量拟合 线性拟合实例线性拟合实例例：例： 下表为实验测得的某一物性和温度之间的关系数据，表中x为温度数据，y为物性数据。请用线性函数拟合温度和物性之间的关系

10、。 解：设拟合直线p(x)=a+bx ,并计算得下表： x7911131517192123252729y91215182124273033363942x313335373941434547 y454851545760636669 编号xyxyx21234521791113154756791215182169819631081652343153243267334981121169225220918389 将数据代入法方程组（1-12）中，得到： 解方程得：a = -1.5 , b = 1.5 拟合直线为： 267338191838956756721bax .p( x ) 5 .151相关系数R为

11、1。单变量拟合单变量拟合 线性拟合实例线性拟合实例线性拟合线性拟合VBVB清单清单Private Sub Command1_Click()Dim x(5),y(5),c,d,m,p,a,b,eerConst n=5For i =1 To 5 x(i)=InputBox(“x(i)=”) y(i)=InputBox(“y(i)=”)Print “x(i)=”;x(i)Print “y(i)=”;y(i)Next I c=0 d=0 m=0 p=0 For i=1 To 5 c=c+x(i) d=d+x(i)2 m=m+y(i) p=p+x(i)*y(i)Next i a=(m*d-c*p)/(n

12、*d-c2) b=(n*p-c*m)/(n*d-c2) 参数计算 a=Int(a*1000+0.5)/1000 b=Int(b*1000+0.5)/1000Text1.Text=Str(a)Text2.Text=Str(b) 参数输出For i=1 To 5 eer=eer+(a+b*x(i)-y(i)2 误差计算 eer=Int(eer*100000+0.5)/100000Next i eer=eer/5Text3.Text=Str(eer)End Sub有关有关线性线性拟合变型问题拟合变型问题 例如要拟合y=a+b/x2，只需在数据输入后增加一语句x(i)=1/x(i)2，而在程序后面的误

13、差eer 的计算中则不需要修改。2.单变量拟合单变量拟合 二次拟合函数二次拟合函数 给定数据（xi ,yi), i=1, 2 , , m ,用二次多项式函数拟合这组数据。 设 ,作出拟合函数与数据序列的均方误差表达式2210 x ax a ap ( x ) 21221012210)()(),(imiiimiiiyxaxaayxpaaaQ由数学知识可知，Q( a0 ,a1 ,a2 )的极小值满足 ：miiiiimiiiiimiiiixyxaxaaaQxyxaxaaaQyxaxaaaQ12221021221011221000)(20)(20)(2整理左式得二次多项式函数拟合的满足条件方程（5-14

14、）： miiimiiimiimiimiimiimiimiimiimiimiiyxyxyaaaxxxxxxxxm121121014131213121121(5-14 ) 解此方程得到在均方误差最小意义下的拟合函数p ( x )。式（5-14）称为多项式拟合的法方程，法方程的系数矩阵是对称的。当拟合多项式n 5时，法方程的系数矩阵是病态的，在用通常的迭代方法求解线性方程时会发散，在计算中要采用一些特殊算法以保护解的准确性。关于线性方程的求解方法，已在第三章中介绍。2.单变量拟合单变量拟合 二次拟合函数二次拟合函数3.二次拟合函数的拓展二次拟合函数的拓展 和一次拟合一样，二次拟合也可以有多种变型，例

15、如 套用上面的公式，可以得到关于求解此拟合函数的法方程（5-15）。值得注意的是在此法方程的构建过程中，进行了变量的代换。首先是拟合函数中变量的代换: 。 P(x)=a0+a1x3+a2x5 253,xxxxmiiimiiimiimiimiimiimiimiimiimiimiiyxyxyaaaxxxxxxxxm1513121011018151815131513（5-15）其次是法方程的代换：将相应拟合函数中的代换引入法方程中。同时应注意法方程中x的4次幂是由两个2次幂相乘得到，x的3次幂是由一个2次幂和一个1次幂相乘得到，而2次幂就是变量本身，而非两个1次幂相乘得到。这个概念至关重要，在以后的

16、二次拟合的各类变型中，均需利用这个概念，千万不要用常规的思路去进行代入计算。 3.二次拟合函数的拓展二次拟合函数的拓展如果我们需要求解是下面的拟合函数： 5 . 1110)273(273lnxbxaaymiiimiiimiimiimiimiimiimiimiimiimiiyxxyyaaaxxxxxxxxm15 . 1112101315 . 015 . 115 . 015 . 1115 . 11ln)273(273lnln)273()273()273()273()273(2731)273(2731参照上面的方法，我们很容易得到求解该拟合函数的法方程3.二次拟合函数的拓展二次拟合函数的拓展4.二次

17、拟合实例二次拟合实例请用二次多项式函数拟合下面这组数据。解解：设 并计算得下表序 号 1 2 3 4 5 6 7 x -3 -2 -1 0 1 2 3 y 4 2 3 0 -1 -2 -5 2210 x ax a ap ( x ) 序号xyxyX2X2yX3X41-34-12936-27812-22-448-8163-13-313-114000000051-1-11-11162-2-44-881673-5-159-45278101-3928-701964.4.二次拟合实例二次拟合实例将上面数据代入式 (5-14) ,相应的法方程为7196028a390280a12807a210210210aa

18、aaaa解方程得： a0 =0.66667 , a1 = -1.39286 , a2 = -0.13095所以：2130950392861666670 x.x- . - .p (x ) -3-2-10123-6-4-2024y=0.66667-1.39286 x-0.13095 x2yY X 图图 5-6 拟合曲线与数据序列拟合曲线与数据序列 二次拟合二次拟合-VB-VB程序清单程序清单Private Sub Command1_Click()Dim n, m As Integern = InputBox(n=方程次数)m = InputBox(m=实验次数)ReDim x(m), y(m),

19、a0(n + 1), a1(n + 1), aa(m, n + 1)ReDim qq(n + 1, m), pp(n + 1, n + 1), b(n + 1), g(n + 1), tt(n + 1, n + 1)For I = 1 To m: 读入数据READ x(I), y(I)Next IData -3, 4, -2, 2, -1, 3, 0, 0, 1, -1, 2, -2, 3, -5For i = 1 To mx(i) = InputBox(x(i)y(i) = InputBox(y(i)Next i二次拟合二次拟合-VB-VB程序清单程序清单omiga = InputBox(o

20、miga=松弛因子)For j = 1 To n + 1为计算法方程中的系数做准备 For i = 1 To m aa(i, j) = x(i) (j - 1) Next iNext jFor j = 1 To m For i = 1 To n + 1 qq(i, j) = aa(j, i) Next iNext j计算法方程中的右边项For i = 1 To n + 1b(i) = 0For j = 1 To m b(i) = b(i) + aa(j, i) * y(j)Next jNext i开始计算法方程中的右边项的系数For i = 1 To n + 1For j = 1 To n +

21、 1 pp(i, j) = 0 For k = 1 To m pp(i, j) = qq(i, k) * aa(k, j) + pp(i, j) Next kNext jNext i二次拟合二次拟合-VB-VB程序清单程序清单开始用松弛迭代法求解法方程中的变量For i = 1 To n + 1 a0(i) = 0 a1(i) = 0.11Next iFor i = 1 To n + 1 g(i) = b(i) / pp(i, i) For j = 1 To n + 1If j = i Then tt(i, j) = 0Else tt(i, j) = -pp(i, j) / pp(i, i)E

22、nd IfNext jNext i50 eer = 0 For i = 1 To n + 1eer = eer + Abs(a0(i) - a1(i) Next iIf eer 2 )。一般地，将含有n个未知量m个方程的线性方程组其一般形式为第四节解矛盾方程组第四节解矛盾方程组mnmnmmnnnnyxaxaxayxaxaxayxaxaxa22112222212111212111 mnmnmmnnyyyxxxaaaaaaaaa2121212222111211 写成矩阵形为 一般情况下，当方程数m多于变量数n ，且m个方程之间线性无关, 则方程组无解，这时方程组称为矛盾方程组。方程组在一般意义下无

23、解，也即无法找到n个变量同时满足m个方程。这种情况和拟合曲线无法同时满足所有的实验数据点相仿，故可以通过求解均方误差 极小意义下矛盾方程的解来获取拟合曲线。 由数学知识还可证明：方程组ATAX = ATY的解就是矛盾方程组AX = Y 在最小二乘法意义下的解。这样我们只要通过求解ATAX = ATY就可以得到矛盾方程组的解，进而得到各种拟合曲线，为拟合曲线的求解提供了另一种方法。第四节解矛盾方程组第四节解矛盾方程组minAX-B22例如，拟合直线p (x ) = a0 +a1x的矛盾方程组 ATAX = ATY的形式如下： 化简得到与式(1-12)相同的法方程第四节解矛盾方程组第四节解矛盾方程

24、组mmmmyyyxxxaaxxxxxx2121102121 111111 111miiimiimiimiimiiyxyaaxxxm11101211 这里需要注意的是变量X和系数（a0 ， a1）之间的相互转换关系。即 对于n次多项式曲线拟合，要计算Q ( a0 ,a1 , , an ) 的极小值问题,这与解矛盾方程组第四节解矛盾方程组第四节解矛盾方程组mxxxAaaX111 ,2110miinini - yx+ a+ x + a a1210) (mnmnmnnnnyxaxaayxaxaayxaxaa102221011110 mnyyyaaaA2110 与求 的极小问题是一回事。minaaa12

25、ini110y-xx或在这里 nmmnnxxxxxxA111 2211故对离散数据（xi,yi）,i=1, 2 , , m ；所作的n次拟合曲线y= xa + +xa +anini10，可通过解下列方程组求得： mTnTyyy A aaaAA2110 （1-21） 第四节解矛盾方程组第四节解矛盾方程组n如果拟合函数有n个自变量并进行一次拟合，则其拟合函数为： （1-22） miinimiinimiiimiiimiinnminiminiminiminiminiminiminiminiminimiimiimiiminimiimiiyxyxyxyxyaaaaaxxxxxxxxxxxxxxxm1111

26、211121012121111121111111131211121xnnnkkxaxaxaxaxaay1122110通过m（mn）次实验，测量得到了m组 ),(, 1,2, 1ininikiiixxxxxy的实数据，则可得到上面n个自变量拟合函数的法方程第四节解矛盾方程组第四节解矛盾方程组n只要对法方程（1-22）稍加修改，就可以得到有n个自变量的任意次方的拟合函数的法方程，通过法方程的求，就可以得到拟合函数中的各项系数。 miiinmiiinmiiimiiimiinnmininmiiinmiiinmiinmiininmiinmiiinmiinmiinimiiimiimiimiinmiimii

27、yxyxyxyxyaaaaxxxxxxxxxxxxxxxxxxxxxm1,1, 112,1, 1111012,1, 2,1, 1,1,1, 112, 11, 1, 11, 11, 11, 2,.112, 11, 11,1, 21, 1（1-23） 第四节解矛盾方程组第四节解矛盾方程组利用解矛盾方程的方法，用二次多项式函数拟合下面数据。解：记二次拟合曲线为 形成法方程第四节解矛盾方程组第四节解矛盾方程组实例实例1 1x -3 -2 -1 0 1 2 3 y 4 2 3 0 -1 -2 -5 2210)(x ax a a xf 721210 yyy A aaaAATT71i471i371i271i

28、371i271i71i271i2772222112722217217 111 111iiiiiiiiTxxxxxxxxxxxxxxxxxxxxAA19602802802807931421111001111421931941014932101231111111 第四节解矛盾方程组第四节解矛盾方程组实例实例1 1而第四节解矛盾方程组第四节解矛盾方程组实例实例1 1739132103249410149321012311111117127171iiiiiiiiTyxyxyYA739119602802802807210aaa得到：解方程得到：a0 = 0.66667 , a1 = -1.39286 ,

29、a2 = -0.13095 f(x)=0.66667-1.39286x-0.13095x2 n例 1.5：给出一组数据，见下表。用解矛盾方程的思路将下面数据拟合成 的经验公式。3 bx a f ( x ) x-3 -2 -1 2 4y14.3 8.3 4.7 8.3 22.7第四节解矛盾方程组第四节解矛盾方程组实例实例2n解：列出法方程： YAbaAATT4954363655111111111151651351335343332313534333231iiiiiiTxxxxxxxxxxxxxAA10623 .587 .223 . 87 . 43 . 83 .14648182711111YAT而：第四节解矛盾方程组第四节解矛盾方程组实例实例2n故法方程为：n解方程得：a = 10.675 , b = 0.137 n拟合曲线为：10623 .58495436365ba3137. 0678.10)(xxf第四节解矛盾方程组第四节解矛盾方程组实例实例2