基于MATLAB环境数学模型参数估计课件.ppt

2023/3/8,1,第6讲：环境数学模型参数估计方法（一）,1、一元、多元线性回归参数估计2、多项式回归参数估计(自学)3、非线性回归参数估计方法 4、非线性最小二乘优化参数估计(自学),2023/3/8,2,一、参数估计重要性,环境数学模型要实现求解，确定其参数是重要的一步。模型中的参数有一些是通过试验手段获得、或者根据文献上资料获得、或者根据经验值获得。上述参数在具体的时间、空间条件下，是否真正的适合某一模型，或者说其适合的程度有多高要靠参数估计的方法评价，并且以这些从各种途径获得的模型参数作为初值，利用参数估计的方法，获得最适合的模型参数是参数估计的最主要任务。参数估计的方法多是一些数学上的优化方法，如最小二乘优化法、非线性优化规划法、梯度最优化算法、线性回归、非线性回归、经验公式等。,2023/3/8,3,二、什么是基于回归分析的参数估计？,回归分析是一种数理统计的方法，用以估计变量之间的相关关系，这种相关关系可能是线性的，也可能是非线性的。回归分析过程是根据因变量和自变量的大量观测数据，发现其大致规律，然后用一定的线性或非线性模型去拟合这些观测数据，回归分析得出因变量和自变量之间的确切数量相依关系式。反过来，如果已经知道因变量和自变量之间数学模型的基本结构，如直线型、幂指数型、多项式型，只是不知道模型中的确切参数，则根据观测值，利用回归分析的方法将模型参数确定出来的过程就是参数估计的过程。利用回归分析技术进行参数估计实际上要比完全意义上的回归分析简单一些。,2023/3/8,4,三、基于一元线性回归的参数估计,数学形式,2023/3/8,5,三、基于一元线性回归的参数估计,数学形式,2023/3/8,6,三、基于一元线性回归的参数估计,应用举例,2023/3/8,7,三、基于一元线性回归的参数估计,求解思路,2023/3/8,8,三、基于一元线性回归的参数估计,简单程序,将上述过程转化为如下的M文件。S=205080100150200250300350;mu=1.18 2.17 2.82 3.03 3.49 3.65 3.97 4.11 4.22;x=1./S;y=1./mu;X=ones(size(x,1),1),x;ab,bint,r,rint,stats=regress(y,X,0.05)mumax=1/ab(1);ks=ab(2)*mumax,运行结果,ab=0.200；12.916；bint=（略去）；r=（略去）；rint=（略去）；stats=0.9995 15412.8412 0.000000000001%可见相关系数平方接近1.0，p值接近0mumax=4.97ks=64.27,2023/3/8,9,四、基于多元线性回归的参数估计,问题描述,2023/3/8,10,四、基于多元线性回归的参数估计,应用举例,2023/3/8,11,四、基于多元线性回归的参数估计,求解思路,2023/3/8,12,四、基于多元线性回归的参数估计,编程实现,M=10000000;u=0.5;A=20;xx=500;%给出已知条件t=1803004806609001140156018002100240030003600;C=141504506246565783933022121476932;y=log(C.*sqrt(t);x1=1./t;x2=t;X=ones(size(t,1),1),x1,x2;%构造因变量自变量矩阵b012,bint,r,rint,stats=regress(y,X,0.05)%多元线性回归T=xx/u;B=b012(3)*(-1),B=(-1)*b012(2)/T2%观察两种途径求得的B是否相等？A0=exp(b012(1)-2*B*T);disp(由B算Dx，);Dx=u2/(4*B),disp(由A0算Dx，);Dx=(M/(A0*A*sqrt(4*pi)2,2023/3/8,13,四、基于多元线性回归的参数估计,运行结果,b012=1.0e+003*0.01239675394754-1.24846558474483-0.00000124753011bint=（略去）；r=（略去）；rint=（略去）stats=1.0e+007*0.00000009999996 1.09793698679399 0B=0.00124753011206B=0.00124846558474通过B计算出Dx，Dx=50.06145204457062通过A0计算出Dx，Dx=50.10180080903253,stats中的第一个元素（R2）非常接近1.0，说明多元回归效果非常好，而且stats中的第三个元素（p值）远小于0.05，印证了回归效果好的结论,2023/3/8,14,五、基于多项式回归的参数估计,问题描述,2023/3/8,15,五、基于多项式回归的参数估计,应用举例,2023/3/8,16,五、基于多项式回归的参数估计,编程求解,【求解】分别多项式回归和多元回归，看二者是否有差别t=1 2 2 3 5 6 7 8 9 10 11 12;y=92 140 213 267 128 343 208 131 188 282 203 143;plot(t,y,k-);gtext(原始曲线);hold on;%进行多项式拟合估计，通过图形观察拟合效果 p4=polyfit(t,y,4);y4=polyval(p4,t);plot(t,y4,-*);gtext(4次多项式);hold on p6=polyfit(t,y,6);y6=polyval(p6,t);plot(t,y6,-o);gtext(6次多项式);hold on p8=polyfit(t,y,8);y8=polyval(p8,t);plot(t,y8,-square);gtext(8次多项式);hold on p9=polyfit(t,y,9);y9=polyval(p9,t);plot(t,y9,-v);gtext(9次多项式);%进行多元回归分析，通过相关系数判断回归效果好坏 pr4,bint4,r4,rint4,stats4=regress(y,t.4,t.3,t.2,t,ones(size(t,1),1)%进行4次多项式回归及效果 pr6,bint6,r6,rint6,stats6=regress(y,t.6,t.5,t.4,t.3,t.2,t,ones(size(t,1),1)pr8,bint8,r8,rint8,stats8=regress(y,t.8,t.7,t.6,t.5,t.4,t.3,t.2,t,ones(size(t,1),1)pr9,bint9,r9,rint9,stats9=regress(y,t.9,t.8,t.7,t.6,t.5,t.4,t.3,t.2,t,ones(size(t,1),1)%比较多项式回归和多元分析结果是否一致，以9次多项式拟合和9次多项式回归对比ppr=p9-pr9%如果ppr=0，则说明多项式拟合和多元回归方法的结果没有差别,2023/3/8,17,五、基于多项式回归的参数估计,运行结果,2023/3/8,18,六、基于非线性回归的参数估计,问题描述,上述讨论的线性回归中的“线性”并非指y与x的关系，而是指y是系数b0、b1、b2等的线性函数,在实际科研工作中，y与参数之间的非线性关系更为常见。,2023/3/8,19,六、基于非线性回归的参数估计,应用举例-1,2023/3/8,20,六、基于非线性回归的参数估计,问题求解,【求解】本例用非线性回归的方法估计参数首先编制M函数文件，描述非线性方程。function mu=c3fun36(k,s)If nargin2;display(输入参数太多);mumax=k(1);ks=k(2);mu=mumax*s./(ks+s);调用非线性回归函数，估计参数。S=20 5080100150200250300 350;mu=1.18 2.17 2.82 3.03 3.49 3.65 3.97 4.11 4.22;k0=1 10;k,R,J=nlinfit(S,mu,c3fun36,k0);%拟合出参数kdisp(前例估算出mumax=4.97 ks=64.27);disp(本次估计出的参数:);mumax=k(1),ks=k(2)kci=nlparci(k,R,J);%获得参数k的置信区间mupred,muci=nlpredci(c3fun36,S,k,R,J);%利用新建的预测模型和原来的S估算u的预测值plot(S,mu,k-,S,mupred,b-*);legend(-观测值;*预测值);%预测值和观测值对比gtext(底物浓度S);gtext(比增长速率u),2023/3/8,21,六、基于非线性回归的参数估计,运行结果,2023/3/8,22,六、基于非线性回归的参数估计,应用举例-2,2023/3/8,23,六、基于非线性回归的参数估计,问题求解,【求解】上述解析解含有余误差函数，其手工计算一般要通过查表的方法，而MATLAB中提供了余误差函数的求解函数erfc()，可以直接实现其求解。%首先编制描述解析解模型的函数function C=c3fun39(Dx,t)c0=350;%mg/Lx=1000;%mu=0.6;%m/sC=(c0/2)*(erfc(x-u*t)./(2*sqrt(Dx*t)+exp(u*x/Dx)*erfc(x+u*t)./(2*sqrt(Dx*t);%然后调用主要函数，进行参数估算。t=60*3 914 2124 2935 3744 5056 60;C=0.000.05 6.00 80.01130.95210.31280.20313.59330.27341.11345.43349.00;Dx0=50;%给出参数初始值Dx=nlinfit(t,C,c3fun39,Dx0);disp(估计出的纵向弥散系数);Dx,2023/3/8,24,六、基于非线性回归的参数估计,求解结果,2023/3/8,25,七、基于非线性最小二乘优化的参数估计,问题描述,基于非线性最小二乘优化的算法在MATLAB中主要有两个函数lsqnonlin()和lsqcurvefit()。如x=lsqnonlin(Fun,x0)的用法就是以初始值x0为搜索起点，发现一个最好的x值，该值使函数Fun平方和最小（即达到目标minsum Fun(x).2时的x）。如果，将目标函数设置为观测值和预测值的差，将目标函数中的参数设置为自变量，同时将x0设置为参数的初始值，那么得到的结果将是使观测值和预测值之差的平方和最小情况下的参数，这个参数是最优的参数。（非原函数）k=lsqcurvefit(Fun,k0,xdata,ydata)任务是以初始值k0为搜索起点，发现一个最好的k值，该值使函数Fun的值和观测值之差的平方和最小（即达到目标minsum(Fun(k,xdata)-ydata).2时的k），其中xdata表示自变量的观测值，ydata表示因变量的观测值，k表示参数变量。(原函数）,2023/3/8,26,七、基于非线性最小二乘优化的参数估计,应用举例,2023/3/8,27,七、基于非线性最小二乘优化的参数估计,问题求解,【求解】首先编制非线性模型的M函数。function p=c3fun310(k,xdata)D=xdata;%方便理解,利用xdata变量表示自变量管径 k1=k(1);k2=k(2);k3=k(3);p=k1+k2*D.k3;调用lsqcurvefit()主函数，开始估算。xdata=0.1 0.150.30.55 0.8 11.21.6;%自变量管道直径ydata=30 5582153294452592930;%因变量价格k0=1 1 1;k=lsqcurvefit(c3fun310,k0,xdata,ydata);%拟合出参数kdisp(估计出的参数:);k1=k(1),k2=k(2),k3=k(3),k1=27.4152,k2=406.3915,k3=1.7137,2023/3/8,28,七、基于非线性最小二乘优化的参数估计,问题描述,2023/3/8,29,七、基于非线性最小二乘优化的参数估计,问题求解,【求解】首先编制非线性模型的M函数。function Zmin=c3fun311(k)D=0.1 0.15 0.30.55 0.8 11.2 1.6;%自变量管道直径观测值pdata=30 55 82153 294452 592 930;%因变量价格观测值k1=k(1);k2=k(2);k3=k(3);p=k1+k2*D.k3;%因变量价格预测值Zmin=pdata-p;%目标函数因变量观测值和预测值的差调用lsqnonlin()主函数，开始估算。k0=1 1 1;k=lsqnonlin(c3fun311,k0);%拟合出参数kdisp(估计出的参数:);k1=k(1),k2=k(2),k3=k(3),k1=27.4152,k2=406.3915,k3=1.7137,2023/3/8,30,已知河流平均流速4.2km/h，饱和溶解氧Os=10mg/L,河流起始点的BOD(L0)浓度为23mg/L,沿程几个断面的溶解氧测定数据如下，根据数据及河流溶解氧变化模式：,，利用excel的规划求解方法，估算河流耗氧速度常数Kd和复氧速度常数Ka。解：C8=C2-(C2-C7)*EXP(-1)*G3*B8/C3)+G4*C4/G5*(EXP(-1*G3*B8/C3)-EXP(-1*G4*B8/C3)1d=24*1h,八、基于EXCEL规划求解的非线性参数估计,2023/3/8,31,思考题,1、,要求用不同的估计方法,2023/3/8,32,思考题,2、,2023/3/8,33,思考题,3、,4、,用非线性回归,用非线性回归,