数学建模 | 您所在的位置:网站首页 › 泄水流量的观测方法有哪些 › 数学建模 |
上篇博客中介绍了插值,插值的一些算法以及Matalb中的一些插值工具函数,这篇博文主要介绍拟合方法以及Matlab中的拟合函数和交互式页面。 一、 拟合方法 线性最小二乘法引入 (1) 概述:与插值相似,在工程中我们常常需要根据两个变量的几组实验数据来找出这两个变量之间的函数关系的近似表达式。我们将找到的近似表达式成为经验公式。但拟合与插值不同的是得到的经验公式不要求经过每一个数据点,而是要求点之间的总偏差最小。 (2) 引例: 为了测定刀具的磨损速度,我们做这样的实验:经过一定的时间(如每隔一小时),测量一次刀具的厚度,得到一组实验数据如下:![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]()
b. 拟合效果的判断: 先在直观判断的基础上,选择几种曲线进行拟合,然后比较,看哪一条曲线的最小二乘法指标J最小 最小二乘法的Matlab实现 (1) 解方程组的方法: 在上述的记号下:![]() ![]() ![]() 输出结果: 1. f = 2. 3. 217/8 - (17*x)/56绘制出图像: 命令行窗口输入以及输出结果: 1. >> [y11,y12]=matching2 2. 3. y11 = 4. 5. 233.4286 6. 7. 8. y12 = 9. 10. 253.9286(3) 最小二乘优化 无约束规划中,若目标函数可以由若干函数的平方和构成,一般这类函数可以写成 最小二乘优化是一类比较特殊的优化问题。在Matlab优化工具箱中提供了一系列函数—— ① lsqlin函数 有约束的线性最小二乘的标准形式是: 求解 求解结果: 1. f=matching4(x0,y0) 2. 3. f = 4. 5. 217/8 - (17*x)/56与前面使用A=R\Y命令效果一致。 ② lsqcurvefit函数 该函数为非线性拟合函数。在已知拟合数据向量xdata和ydata,并且知道xdata和ydata之间的函数关系ydata=F(x,xdata),仅仅不知道函数的某些系数系数向量x,可以用lsqcurvefit函数进行对未知系数的拟合,求x使下面的式子成立: 命令行窗口输入: 1. y0=matchingData(:,2); 2. x0=matchingData(:,[3:5]); 3. parameter0=rand(2,1);%拟合参数的初始值是任意取的 4. %非线性拟合的答案是不唯一的,下面给出拟合参数的上下限 5. lb=[2,1]; 6. ub=[20;2]; 7. %上下限不唯一,可以自己设定,不设定时默认值是空矩阵 8. parameter=lsqcurvefit(@matchingf,parameter0,x0,y0,lb,ub)输出结果(不唯一) 1. Initial point is a local minimum. 2. 3. Optimization completed because the size of the gradient at the initial point 4. is less than the value of the optimality tolerance. 5. 6. 7. 8. parameter = 9. 10. 11.0000 11. 1.5000例2: 用最小二乘法拟合 解:Matlab程序如下: 1. >> load('data3.mat') 2. >> mf=@(parameter,xdata)1/sqrt(2*pi)/parameter(2)*exp(-(xdata-parameter(1)).^2/2); 3. >> parameter=lsqcurvefit(mf,rand(2,1),x0,y0)输出结果: 1. Local minimum found. 2. 3. Optimization completed because the size of the gradient is less than 4. the value of the optimality tolerance. 5. 6. 7. 8. parameter = 9. 10. 0.0000 11. 1.0000③ lsqnonlin函数 非线性最小二乘(非线性拟合)的标准形式是: 输出结果: 1. Local minimum found. 2. 3. Optimization completed because the size of the gradient is less than 4. the value of the optimality tolerance. 5. 6. 7. 8. parameter = 9. 10. 0.0000 11. 1.0000④ lsqnonneg函数 非负线性最小二乘的标准形式是 运行结果: - >> run2 - - x = - - 0 - 0.6929⑤ Matlab的曲线拟合用户图形界面解法 Matlab提供了命令cftool,该命令给出了一维数据拟合的交互式环境。具体执行步骤如下: 把数据导入到工作空间运行cftool,打开用户图形界面窗口选择适当的模型进行拟合生成一些相关的统计量进行预测具体用法参见:CSDN某位daolao的博客![]() ![]() ![]() ![]() ![]() 1.1 黄河小浪底调水调沙问题 问题的提出: 2004年6月至7月黄河进行了第三次调水调沙实验,特别是首次由小浪底、三门峡和万家寨三大水库联合调度,采用接力式防洪预泄放水,形成人造洪峰进行调沙实验获得成功,整个试验期为20多天,小浪底从6月19日开始预泄放水,直到7月13日结束并恢复正常供水。小浪底水利工程按设计拦沙量为75.5亿m3,在这之前,小浪底共积泥沙达14.15亿t.这次调水调沙试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪,在小浪底形成人造洪峰,冲刷小浪底库区沉积的泥沙,在小浪底水库开闸泄洪以后,从6月27日开始三门峡水库和万家寨水库陆续开闸放水,人造洪峰于6月29日先后到达小浪底,7月3日达到最大流量2700m3/s,使小浪底水库的排沙量也不断增加。表1是由小浪底观测站从6月29日到7月10日检测到的试验数据。 其中i表示次数,计时单位:秒/s 根据题意得出等式:某一时刻排沙量=水流量×含沙量 记第i,(i=1,2,…,24)次观测时的排沙量为y_i,水流量为v_i,含沙量为c_i 则第i次观测时的排沙量为: y_i=v_i×c_i 由微分和积分的关系可知,任意时刻的总的排沙量等于各时刻的排沙量之和,记总的排沙量为z,则有: 命令行窗口输入命令以及输出结果: 1. >>[f,total]=Mymatching(HuangHe) 2. 3. f = 4. 5. 23×1 cell 数组 6. 7. {'1.446373*x - 0.000003145275*x^2 - 2.648125e-12*x^3 + 57600.0' } 8. {'1.159796*x - 0.000003488472*x^2 - 1.151619e-12*x^3 + 114000.0' } 9. {'0.8519441*x - 0.000003637722*x^2 - 6.389389e-12*x^3 + 157500.0' } 10. {'0.5018725*x - 0.000004465786*x^2 + 8.252551e-11*x^3 + 187000.0' } 11. {'0.5780658*x + 0.000006229519*x^2 - 1.041684e-10*x^3 + 207000.0' } 12. {'0.5330865*x - 0.000007270706*x^2 + 6.622972e-11*x^3 + 235200.0' } 13. {'0.2756992*x + 0.000001312666*x^2 + 1.041959e-11*x^3 + 250000.0' } 14. {'0.4474499*x + 0.000002663045*x^2 - 4.09285e-11*x^3 + 265200.0' } 15. {'0.4483899*x - 0.000002641288*x^2 + 2.181593e-11*x^3 + 286200.0' } 16. {'0.3423239*x + 0.0000001860564*x^2 - 5.873884e-11*x^3 + 302400.0' } 17. {'0.02953683*x - 0.000007426497*x^2 + 8.910315e-11*x^3 + 312800.0' } 18. {'0.000004121271*x^2 - 0.113249*x - 4.215901e-11*x^3 + 307400.0' } 19. {'0.006792364*x - 0.000001342536*x^2 - 5.690704e-11*x^3 + 306800.0'} 20. {'7.629057e-11*x^3 - 0.000008717689*x^2 - 0.4278094*x + 300000.0' } 21. {'0.000001169569*x^2 - 0.7538882*x - 1.24219e-10*x^3 + 271400.0' } 22. {'1.87397e-10*x^3 - 0.00001492921*x^2 - 1.348305*x + 231000.0' } 23. {'0.000009357448*x^2 - 1.589005*x + 2.706168e-11*x^3 + 160000.0' } 24. {'0.00001286464*x^2 - 0.6290103*x - 2.088184e-10*x^3 + 111000.0' } 25. {'2.376449e-10*x^3 - 0.00001419822*x^2 - 0.6866208*x + 91000.0' } 26. {'0.00001660056*x^2 - 0.5828398*x - 1.773961e-10*x^3 + 54000.0' } 27. {'3.161053e-11*x^3 - 0.000006389971*x^2 - 0.1417424*x + 45500.0' } 28. {'5.715577e-11*x^3 - 0.000002293247*x^2 - 0.5168574*x + 30000.0' } 29. {'0.000005114141*x^2 - 0.3949948*x + 4.985713e-11*x^3 + 8000.0' } 30. 31. 32. total = 33. 34. 1.8440e+11
(2) 确定排沙量与水流量的关系。 由上拟合曲线图可知排沙量随着时间的变化大致分为两个阶段: 7月4日前排沙量不断增加,之后排沙量不断减少,而根据表格可以得到水流量在前一段时间随时间推进而增加,后一段时间随着时间推进而减少,所以我们大致可以将水流量与排沙量之间的关系分为两个阶段分别进行拟合。 我们对于两个阶段都分别采用一次拟合和二次拟合来进行,并且以标准差来作为最后选取的标准: 6月29日~7月4日y=250.5655v-373384 7月5日~7月10日y=0.1067v^2-180.4668v+72421.0982 实现Matlab代码如下: 1. function speedWithSand(data) 2. %准备两个阶段的数据 3. syms x 4. data1=data([1:12],:); 5. data2=data([13:23],:); 6. speed1=data1(:,1); 7. sand1=data1(:,2).*speed1; 8. speed2=data2(:,1); 9. sand2=data2(:,2).*speed2; 10. %进行第一分阶段的拟合 11. format long e 12. rmse1=ones(2,1); 13. rmse2=ones(2,1); 14. cha1=ones(2,1); 15. cha2=ones(2,1); 16. for j=1:2 17. match1{j}=polyfit(speed1,sand1,j); %拟合多项式 18. forcast1{j}=polyval(match1{j},speed1);%求出预测值 19. %以下求误差平方和以及剩余标准差 20. cha1(j)=sum((sand1-forcast1{j}).^2); 21. rmse1(j)=sqrt(cha1(j)/(10-j)); 22. end 23. celldisp(match1) 24. rmse1 25. 26. for j=1:2 27. match2{j}=polyfit(speed2,sand2,j); %拟合多项式 28. forcast2{j}=polyval(match2{j},speed2);%求出预测值 29. %以下求误差平方和以及剩余标准差 30. cha2(j)=sum((sand2-forcast2{j}).^2); 31. rmse2(j)=sqrt(cha2(j)/(10-j)); 32. end 33. celldisp(match1) 34. rmse2 35. format命令行窗口输入以及输出结果: 1. >> speedWithSand(HuangHe) 2. 3. match1{1} = 4. 5. 2.546957330612737e+02 -3.818018589089608e+05 6. 7. 8. 9. match1{2} = 10. 11. -4.973604894806963e-02 4.821949815832189e+02 -6.369512800128722e+05 12. 13. 14. 15. rmse1 = 16. 17. 1.135270235952025e+04 18. 1.111634607843465e+04 19. 20. 21. match1{1} = 22. 23. 2.546957330612737e+02 -3.818018589089608e+05 24. 25. 26. 27. match1{2} = 28. 29. -4.973604894806963e-02 4.821949815832189e+02 -6.369512800128722e+05 30. 31. 32. 33. rmse2 = 34. 35. 4.378235683047678e+04 36. 3.053573018955462e+04 |
CopyRight 2018-2019 实验室设备网 版权所有 |