《水塔水流量的估计》课程设计论文

数值分析课程设计论文 水塔水流量的估计水塔水流量的估计 摘摘要要: :数学建模方法是处理科学理论的一种经典方法，也是解决各类实际问题的常用 方法。本文采用曲线拟合的方法，并利用数学软件 MATLAB 对水塔流量进行计算，计 算结果与实际记录基本吻合。 关键词关键词: :建模，流量，拟合，MATLAB 1. 1.问题重述问题重述 美国某州的各用水管理机构要求各社区提供用水率（以每小时多少加仑计，英制单 位下，1 加仑=4.54596dm3，美制单位下，1 加仑=3.78533dm3）以及每天所用的总用水 量，但许多社区并没有测量流入或流出当地水塔的水量的设备，而只能以每小时测量水 塔的水位代替，其精度在 0.5%以内。更为重要的是，无论什么时候，只要水塔中的水 位下降到某一最低水位 L 时，水泵就启动向水塔重新充水直至某一最高水位 H，但也无 法得到水泵的供水量的测量数据。因此，在水泵正在工作时，不容易建立水塔中水位与 水泵工作时用水量之间的关系。水泵每天向水塔充水一次或两次，每次大约2 小时。试 估计在任何时候， 甚至包括水泵正在工作的时间内从水塔流出的流量ft， 并估计一天 的总用水量。水塔是一个垂直圆柱体，高为 40 英尺，直径为 57 英尺。 下表给出了某个小镇某一天的真实数据： 表表 1 1：某小镇某天的水塔水位（：某小镇某天的水塔水位（1m=3.2811m=3.281 英尺）英尺） 时间（秒）时间（秒）水位（英尺）水位（英尺）时间（秒）时间（秒）水位（英尺）水位（英尺）时间（秒）时间（秒）水位（英尺）水位（英尺） 0 3316 6635 10619 13937 17921 21240 25223 28543 32284 31.75 31.10 30.54 29.94 29.55 28.92 28.50 27.87 27.52 26.97 35932 39332 39435 43318 46636 49953 53936 57254 60574 64554 水泵工作 水泵工作 35.50 34.45 33.50 32.67 31.56 30.81 30.12 29.27 68535 71854 75021 79154 82649 85968 89953 93270 28.42 27.67 26.97 水泵工作 水泵工作 34.75 33.89 33.40 2. 2.问题分析问题分析 数据的单位转换： 表表 2 2 1 数值分析课程设计论文 时间（时间（h h）水位（水位（mm）时间（时间（h h）水位（水位（mm）时间（时间（h h）水位（水位（mm） 0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.6769 9.4788 9.3081 9.1253 9.0071 8.8144 8.6864 8.5030 8.3877 8.2201 9.98 10.93 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.94 水泵工作 水泵工作 10.8199 10.4998 10.2103 9.9573 9.6190 9.3904 9.1801 8.9211 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91 8.6620 8.4334 8.2201 水泵工作 水泵工作 10.5913 10.3292 10.1798 流量是单位时间流出的水的体积，可以由对应时刻的流速乘以水塔的横截面积得 到。由于水塔是正圆柱形，横截面积是常数 S，所以我们在这里研究的其实是流速的变 化。 在水泵不工作的时段，流量很容易从水位对时间的变化率 ,即流速算出，问题是如 何估计水泵供水时段的流速。水泵供水时段的流速只能靠供水时段前后的流速拟合得 到，作为用于拟合的原始数据，我们希望水泵不工作的时段流速越准确越好。这些流速 大体可由两种方法计算: 一是直接对表 2 中的水位用数值微分算出各时段的流速， 用它 们拟合其它时刻或连续时间的流速。二是先用表中数据拟合水位 -时间函数，求导数即 可得到连续时间的流速。 一般说来数值微分的精度不高，何况测量记录还是不等距的，数值微分的计算尤其 麻烦。下面我们用第二种方法处理。 有了任何时刻的流速，就不难计算一天的总用水量。其实，水泵不工作时段的用水 量可以由测量记录直接得到，如表 2 可知从 t=0 到 t=8.97(h)水位下降了 9.6769 8.22011.4568(m )，乘以水塔的横截面积 S 就是这一时段的用水量。这个数值 可以用来检验拟合的结果。 水泵第 1 次供水时段为 t=8.98(h)到 t=10.94(h),第 2 次供水时 段为 t=21(h)到 t=23 (h)。这是根据最低和最高水位分别是 8.2201m 和 10.8199m 及表 2 的水位测量记录作出的假设。其中前3 个时刻取自实测数据（精确到0.01h）,最后 1 个 时刻来自每次供水约两小时的已知条件 （从记录看， 第2次供水时段应在有记录的22.96h 之后不久结束） 。 3. 3.模型假设模型假设 2 数值分析课程设计论文 供水时段的假设供水时段的假设 水泵工作时单位时间的供水量基本为常数，这个常数大于单位时间的平均流量。流 量是单位时间流出水的体积，这里假设是水位对时间的连续函数，即h ht。为简化 处理，不影响问题的解决，假设流量与水泵是工作无关。 4. 4.流量估计流量估计 4.14.1 拟合水位拟合水位- -时间函数时间函数 从表 2 测量记录看，一天有两个供水时段（以下称第 1 供水时段和第 2 供水时段） 和 3 个水泵不工作时段(以下称第 1 用水时段 t=0 到 t=8.97，第 2 用水时段 t=10.95 到 t=20.48 和第 3 用水时段 t=23 以后)。 对第 1、2 用水时段的测量数据分别作多项式拟合， 得到水位函数h 1 h 1 t和h 2 h 2 t。为使拟合曲线比较光滑，多项式次数不要太高，一 般用 36 次。由于第3 时段只有 3 个测量记录，无法对这一时段的水位作出比较好的拟 合，可采用外推的办法解决。 4.24.2 确定流量确定流量- -时间函数时间函数 对于第 1、2 用水时段，只需将水位函数h i h i t ,i 1,2求导数即可，对于两个供水 时段的流速，则用供水时段前后（水泵不工作时段）的流速拟合得到，并且将拟合得到 的第 2 供水时段流速外推，将第3 用水时段流速包含在第 2 供水时段内，需要拟合四个 流速函数。 4.34.3 一天的总用水量一天的总用水量 总用水量等于两个水泵不工作时段和两个供水时段用水量之和， 它们都可以由流速 对时间的积分得到： V Vdt Shdt t0t0 tt 5. 5.算法设计与计算结果算法设计与计算结果 5.15.1 拟合第拟合第 1 1、2 2 时段的水位，并得出流量时段的水位，并得出流量 5.1.1 第 1 用水时段的流速 设 t、h 为已输入的时刻和水位测量记录，实现如下： 3 数值分析课程设计论文 t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83, 17.94,19.04,19.96,20.84,23.88,24.99,25.91; h= 9.6769,9.4788,9.3081,9.1253,9.0071,8.8144,8.6864,8.5030,8.3877,8.2201,10.8199,10.49 98,10.2103,9.9573,9.6190,9.3904,9.1801,8.9211,8.6620,8.4334,8.2201,10.5913,10.3292,10.1 798; f1=polyfit(t(1:10),h(1:10),5); %用 5 次多项式拟合第 1 用水时段水位h 1 h 1 t，f1 输出 5 次多项式的系数 b1=polyder(f1);% b1输出多项式（系数为f1）导数的系数，给出水位变化率 tm1=0:0.01:8.97;%将第一用水时段0 , 8.97细分 g1=-polyval(b1,tm1);%g1输出多项式b1在tm1点的函数值(取负后边为正值)， 即tm1时刻 的流速 5.1.2 第2用水时段的流速 实现如下： f2=polyfit(t(11:21),h(11:21),5);%用5次多项式拟合第2用水时段水位水位h 2 h 2 t，f2 输出5次多项式的系数 b2=polyder(f2);%b2输出多项式（系数为f2）导数的系数,给出水位变化率 tm2=10.95:0.01:20.84;%将第二用水时段10.95 , 20.84细分 g2=-polyval(b2,tm2);%g2输出多项式（b2）在tm2点的函数值(取负后边为正值)，即tm2 时刻的流速 4 数值分析课程设计论文 第1、2用水时间段的流速曲线图: 5.25.2拟合供水时段的流速拟合供水时段的流速 5.2.1第1供水时间段的流速 在第1供水时段（t = 8.9810.94）之前（即第1用水时段）和之后（第2用水时段） 各取几点，其流速已经得到，用它们拟合第1供水时段的流速。为使流速函数在供水时 段连续，只取4个点，拟合5次多项式（即曲线必过这4个点）。拟合5次多项式，实现如 下： q1=-polyval(b1,7.93,8.97);%取第1时段在t=7.93,8.97的流速 q2=-polyval(b2,10.95,12.03);%取第2时段在t=10.95,12.03的流速 dx=7.93,8.97,10.95,12.03; dy=q1,q2;%将四个点合并 d=polyfit(dx,dy,5) ;%拟合5次多项式 ex=8.97:0.01:10.95; %将第一供水时段8.97,10.95细分 ey=polyval(d,ex); %ey输出第一供水时段各时刻的流速 5 数值分析课程设计论文 5.2.2第2供水及第3用水(到t=24)时间段的流速 在第2供水时段之前取t =20,20.8两点的流速，在该时刻之后（第3用水时段）仅有3 个水位记录，我们用差分得到流速，然后用这4个数值拟合第2供水时段的流速。实现如 下： dt3=diff(t(22:24);%最后3个时刻的两两之差：t t i1 t i dh3=diff(h(22:24);%最后3个水位的两两之差：h h i1 h i dht3=-dh3/dt3;% t(22)和t(23)的流速（差商代替微商） t3=19.96,20.84,t(22),t(23);%取第2供水时段前后各两点 a=-polyval(b2,t3(1:4),dht3);%求得t3各时刻的流速 h=polyfit(t3,a,5);% 拟合5次多项式 mx=20.84:0.01:24.00;%将第2供水时段和第3用水时段细分 my=polyval(h,mx);% my输出第2供水时段（到t=24）各时刻的流速 第1供水时间段，第2供水及第3用水(到t=24)时间段的流速曲线图: 5.35.3一天的总用水量的估计一天的总用水量的估计 第1、2用水时段和第1、2供水时段流速的积分之和乘以S，就是一天总用水量。虽 6 数值分析课程设计论文 然诸时段的流量已表示为多项式函数，积分可以解析的算出，这里仍用数值积分计算。 5.3.1第1用水时间段的用水量： v 1 t1 t0 dt Sh 1 dt dt Sh 1 v 1 Sh 1 t0t0 t1t1 其中积分值 h1通过梯形公式计算： (t)dt h 1 (t)dt t h 1 h 1 t0 k1 tk1 t1 N tk (t k1)h1 (t k )h 1 2 k1 N 计算得出第1用水时间段的用水量：v 1 363.0086m3 5.3.2第2用水时间段的用水量： dt Sh 2 dt Sh 2 dt Sh 2 v 2 v 2 t2t2t2 t3t3t3 (t)dt h 2 (t)dt t h 2 h 2 t2 k1 tk1 t3 N tk (t k1)h2 (t k )h 2 2 k1 N 计算得出第2用水时间段的用水量：v 2 647.9356m3 5.3.3第1供水时间段的用水量： v 12 t2 t1 t2 t1 dt Sh 12 dt Sh 12 dt Sh 12 v 12 t1t1 N tk tk1 t2t1 (t)dt h 12 (t)dt t h 12 h 12 k1 (t k1)h12 (t k )h 12 2 k1 N 计算得出第1供水时间段的用水量：v 12 103.4852m3 5.3.4第2供水和第3用水(到t=24)时间段的用水量 dt Sh 3 dt Sh 3 dt Sh 3 v 3 v 3 t3t3t3 t4t4t4 (t)dt h 3 (t)dt t h 3 h 3 t0 k0 tk1 t N tk (t k1)h3 (t k )h 3 2 k0 N 计算得出第2供水和第3用水(到t=24)时间段的用水量v 3 146.7717m3 一天总用水量：v v 1 v 2 v 12 v 3 1261.2m3 6. 6.流速及总用水量的检验流速及总用水量的检验 7 数值分析课程设计论文 6.16.1 用水时段流速的检验方法：用水时段流速的检验方法： 计算出的各时刻的流速可用水位记录的数据来检验。用水量 V1 用第 1 用水时段水 位测量记录中下降高度H 1 9.67698.22011.4568m来计算并检验，在第一用水时段 水的实际用量为:V 1 S H 1 363.0275m3 类似地，第二用水时段用去的水的高度H 2 10.81998.2201 2.5998m 实际用水量为:V 2 S H 2 647.8756m3， 而通过算法得到的近似值：v 1 363.0068m3，v 2 647.9356m3 将v 1 、v 2 与V 1 、V 2 进行比较，两者相差无几。 6.26.2 供水时段流速的检验方法：供水时段流速的检验方法： 通过求水泵在两个供水时间段的功率p 1 与p 2 , p 1 S y S 10.59138.2201 S y 12 S 10.8199 8.2201 ，p 2 3,而两个供水时段 2.0392 水泵的功率应大致相等。 第 1、2 时段水泵的公率计算如下： p1=y12*S+(10.8199-8.2201)*S/2.039%第1供水时段水泵的功率(水量以高度计） tm3=20.84:0.01:23.88; g3=polyval(h,tm3);% g3输出第2供水时段各时刻的流速 p2=(0.01*trapz(my)+(10.5913-8.2201)*S/2%第 2 供水时段水泵的功率（水量仍以高度 计） ，p 2 368.8317第 1、2 时段水泵的公率计算结果：p 1 368.4859 通过算法得出的第 1、2 用水时段的总用水量与实际的第 1、2 用水时段的总用水量 非常接近，而两个供水时段水泵的功率相当，说明该算法切实可行！ 8 数值分析课程设计论文 7. 7.分析分析 由表 2 可以看出第 1、2 用水时间段的总用水量、水位下降的高度与实际的总用水 量、水位下降高度相差无几。所以数据拟合，数值积分取得精度总够的。由水塔一天的 流速曲线图可以看出，流速曲线与原始记录基本上吻合，零点到 10 点钟用水量最少， 10 点到下午 3 点是用水高峰期。 参考文献参考文献 1 李大潜 中国大学生数学建模竞赛M 北京 高等教育出版社 1998 2 叶其孝 大学生数学建模竞赛辅导M 长沙 湖南教育出版社 1993 3 萧树铁 数学实验M 北京 高等教育出版社 1999 9 数值分析课程设计论文 4 乐经良 数学实验M 北京 高等教育出版社 2005 5 Justin matlab 数学实验 附录附录 部分源代码：部分源代码： S=(57/3.2)/2)2*pi;S=(57/3.2)/2)2*pi; t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83, 17.94,19.04,19.96,20.84,23.88,24.99,25.91;17.94,19.04,19.96,20.84,23.88,24.99,25.91; h=h= 9.6769,9.4788,9.3081,9.1253,9.0071,8.8144,8.6864,8.5030,8.3877,8.2201,10.8199,10.499.6769,9.4788,9.3081,9.1253,9.0071,8.8144,8.6864,8.5030,8.3877,8.2201,10.8199,10.49 98,10.2103,9.9573,9.6190,9.3904,9.1801,8.9211,8.6620,8.4334,8.2201,10.5913,10.3292,10.198,10.2103,9.9573,9.6190,9.3904,9.1801,8.9211,8.6620,8.4334,8.2201,10.5913,10.3292,10.1 798;798; f1=polyfit(t(1:10),h(1:10),5);f1=polyfit(t(1:10),h(1:10),5); b1=polyder(f1);b1=polyder(f1); tm1=0:0.01:8.97;tm1=0:0.01:8.97; g1=-polyval(b1,tm1);g1=-polyval(b1,tm1);%第第 1 1 用水时段用水时段 f2=polyfit(t(11:21),h(11:21),5);f2=polyfit(t(11:21),h(11:21),5); b2=polyder(f2);b2=polyder(f2); tm2=10.95:0.01:20.84;tm2=10.95:0.01:20.84; g2=-polyval(b2,tm2);g2=-polyval(b2,tm2);%第第 2 2 用水时段用水时段 q1=-polyval(b1,7.93,8.97);q1=-polyval(b1,7.93,8.97); q2=-polyval(b2,10.95,12.03);q2=-polyval(b2,10.95,12.03); dx=7.93,8.97,10.95,12.03;dx=7.93,8.97,10.95,12.03; dy=q1,q2;dy=q1,q2; 10 数值分析课程设计论文 d=polyfit(dx,dy,5) ;d=polyfit(dx,dy,5) ; ex=8.97:0.01:10.95;ex=8.97:0.01:10.95; ey=polyval(d,ex);ey=polyval(d,ex);%第第 1 1 供水时段供水时段 dt3=diff(t(22:24);dt3=diff(t(22:24); dh3=diff(h(22:24);dh3=diff(h(22:24); dht3=-dh3/dt3;dht3=-dh3/dt3; t3=19.96,20.84,t(22),t(23);t3=19.96,20.84,t(22),t(23); a=-polyval(b2,t3(1:4),dht3);a=-polyval(b2,t3(1:4),dht3); h=polyfit(t3,a,5);h=polyfit(t3,a,5); mx=20.84:0.01:24.00;mx=20.84:0.01:24.00; my=polyval(h,mx);my=polyval(h,mx);%第第 2 2 供水及第供水及第 3 3 用水（到用水（到 t=24t=24）的时段）的时段 y1=0.01*trapz(g1);y1=0.01*trapz(g1); y2=0.01*trapz(g2);y2=0.01*trapz(g2); y12=0.01*trapz(ey);y12=0.01*trapz(ey); y3=0.01*trapz(my);y3=0.01*trapz(my); v=(y1+y2+y12+y3)*Sv=(y1+y2+y12+y3)*S%一天的总用水量一天的总用水量 v1=y1*Sv1=y1*S v2=y2*Sv2=y2*S V1=(9.6769-8.2201)*SV1=(9.6769-8.2201)*S V2=(10.8199-8.2201