居民区供水问题(共17页)

1、精选优质文档-倾情为你奉上居民区供水问题的研究与分析摘要：本文针对居民区供水问题进行研究与讨论。根据拟合原理对水塔的流量进行模拟，运用数值积分对一天的总用水量进行估计。本文建立了水塔流量模拟的数学模型，并对模型进行了分析与评价。 对于问题一，流量是单位时间流出的水的体积，由于水塔是正圆柱形，横截面积是常数，在水泵不工作的时段，流量很容易从水位对时间的变化率算出，水泵供水时段的流量只能靠供水时段前后的流量拟合得到，在本文的计算中，我们先表中所给数据拟合水位时间函数，对此函数进行求倒，得到连续时间的流量函数。对于问题二，运用问题一得出的结果，分别对第1、2时段和第1、2供水时段流量进行积分，然后求

2、和，得出的结果即是一天的总用水量。对于问题三,根据水泵的功率（单位时间泵如的水量）等于供水时段的用水量加上水位上升值除以时段长度，可以很容易的得出水泵的功率。关键词：水塔 流量 拟合 模拟 数值积分一、问题的提出    某居民区的民用自来水是由圆柱形水塔提供,水塔高12.2米,直径17.4米.水塔是由水泵根据水塔内水位高低自动加水,一般每天水泵工作两次.现在需要了解居民区用水规律与水泵的工作功率.按照设计,当水塔的水位降至最低水位,约8.2米,水泵自动启动加水;当水位升高到一个最高水位, 约10.8米,水泵停止工作.    

3、可以考虑采用用水率(单位时间的用水量)来反映用水规律,并通过间隔一段时间测量水塔里的水位来估算用水率,表1是某一天的测量记录数据,测量了28个时刻,但是由于其中有3个时刻遇到水泵正在向水塔供水,而无水位记录(表1中用/表示).    试建立合适的数学模型,推算任意时刻的用水率,一天的总用水量和水泵工作功率.    表1 原始数据(单位:时刻(小时),水塔中水位(米)时刻t00.9211.8432.9493.8714.9785.900水位9.6779.4799.3089.1258.9828.8148.686时刻t7.0067.9288.

4、9679.981110.92510.95412.032水位8.5258.3888.220/10.82010.500时刻t12.95413.87514.98215.90316.82617.93119.037水位10.2109.9369.6539.4099.1808.9218.662时刻t19.95920.83922.01522.95823.88024.98625.908水位8.4338.220/10.82010.59710.35410.180二、模型的基本假设1、 流量只取决于水位差，与水位本身无关，由于从小孔流出的流体的流速正比于水面高度的平房根，由题目给出的水塔的最低和最高水位分别为8.2m

5、和10.8m（设出口的水位为0），所以可以忽略水位对流速的影响。因为2、 水泵第一次供水时段为到，第二次供水时段为到，这是根据最低和最高水位分别是8.2m和10.8m，及表一的水位测量记录做出的假设，其中前3个时刻直接取自实测数据（精确到0.1h），最后一个时刻来自每次供水约两个小时的已知条件（从记录看，大约在20.84h之后开始供水，大约在22.9h之前结束供水）。3、 假设水泵工作时单位时间的供水量是常数，此常数大于单位时间的平均流量。4、 假设流量是关于时间的连续函数。5、 假设流量与水泵是否工作无关。6、 由于水塔截面积是常数，计算中将流量定义为单位时间流出的水的高度，即水位对时间变化

6、率的绝对值，将最后的乘以S即可。即：水位是时间的连续函数 水位对时间的变化率（流量） 任何时刻的流量： 三：符号说明：测量时刻：水塔截面积：水位测量记录：第1时段的用水量：第2时段的用水量：第1供水时段用水量：第2供水时段用水量：一天总用水量：第1时段水位下降高度：第2时段水位下降高度：第1供水时段水位下降高度：第2供水时段和第3时段水位下降高度四：问题的分析与准备1、 拟合水位时间函数从表1测量记录看，一天有两个供水时段（简称第1供水时段和第2供水时段），和三个水泵不工作时间段（以下简称第1时间段t=0到t=8.967，第2时间段t=10.954到t=20.839和第3时段t=23以后）。对

7、第1，2时段的测量数据直接分别作多项式拟合，得到水位函数。为使拟合曲线比较光滑，多项式次数不宜太高，一般在36.由于第三时段只有三个测量记录，无法对这一时段的水位做出较好的拟合。2、 确定流量时间函数对于第1、2时段只要将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第二供水时段流量外推，将第3时段流量包含在第2供水时段内，需要拟合四个流量函数。3、 一天总用水量的估计总用水量等于两个水泵不工作时段和两个供水时段用水量之和，他们都可以由流量对时间的积分得到：五：模型的建立与求解1、拟合第1、2时段的水位，并导出流量 设t，h为已输

8、入的时刻和水位测量记录（水泵启动的4个时刻不输入），为了计算方便，也为了提高精度，将水位的单位化为cm，时间精确到0.01h。t=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 0 0 1082 1050 1021 994 965 941 918 892

9、866 843 822 0 0 1059 1035 1018;第1时段各时刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）； %用3次多项式拟合第1时段水位，c1输出3次多项式的系数2）a1=polyder（c1）；% a1输出多项式（系数为c1）导数的系数3）tp1=0：0.1：9；x1=-polyval（a1，tp1）；% x1输出多项式（系数a1）在tp1点的函数值（取负后边为正值），即tp1时刻的流量。4）流量函数为：具体程序见附录shuita1.m用程序进行计算，求得在0，9内各时刻的流量值(水位变化率)如下表：ti00.921.842.953.874

10、.985.907.017.938.97h122.1119.8117.9116.1415.1214.4214.2814.6415.3816.69fig.1 第一用水时段流量曲线图第2时段各时刻的流量可如下得： 1）c2=polyfit(t(13:23),h(13:23),4) %用4次多项式拟合第2时段水位，c2输出4次多项式的系数2）a2=polyder(c2) % a2输出多项式（系数为c2）导数的系数3）tp2=10.9:0.1:21 x2=-polyval(a2,tp2) % x2输出多项式（系数为a2）在tp2点的函数值（取负后边为正值），即tp2时刻的流量4）流量函数为：具体程序见附

11、录shuita2.m用程序进行计算，求得在11,20.8内各时刻的流量值(水位变化率)如下表：ti10.9512.0312.9513.8814.9815.9016.8317.9319.0419.9620.84h131.8030.5329.2927.9926.5025.4024.5223.9123.9024.4625.56fig.2 第二用水时段流量曲线图2. 拟合供水时段的流量 在第1供水时段（t=911）之前（即第1时段）和之后（即第2时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量为使流量函数在t=9和t=11连续，我们简单地只取4个点，拟合3次多项式（即曲线必过这4个点），实

12、现如下：xx1=-polyval(a1,8 9);%取第1时段在t=8,9的流量 xx2=-polyval(a2,11 12);%取第2时段在t=11,12的流量 xx12=xx1 xx2; c12=polyfit(8 9 11 12,xx12，3);%拟合3次多项式 tp12=9：0.1：11; x12=polyval(c12,tp12);%x12输出第1供水时段各时刻的流量 拟合的流量函数为： 具体程序见附录shuita3.m用程序进行计算，求得在9,11内各时刻的流量值(水位变化率)如下表：ti911h116.73831.72 fig.3 第一供水时段流量曲线图在第2供水时段之前取t=2

13、0，20.8两点的流水量，在该时刻之后（第3时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量如下： dt3=diff（t(22：24)）； %最后3个时刻的两两之差 dh3=diff（h(22：24)）； %最后3个水位的两两之差dht3=-dh3./dt3； %t(22)和t(23)的流量t3=20 20.8 t(22) t(23)； xx3=-polyval(a2，t3(1：2)，dht3)； %取t3各时刻的流量 c3=polyfit（t3，xx3，3）； %拟合3次多项式 t3=20.8：0.1：24； x3=polyval（c3，tp3）； % x3

14、输出第2供水时段（外推至t=24）各时刻的流量拟合的流量函数为： 具体程序见shuita4.m 用程序进行计算，求得在20.8,24内各时刻的流量值(水位变化率)如下表： ti22.0122.9623.8824.9925.91h125.57423.93721.62218.478 16.08 fig.4 第2供水时段和第3用水时段流量曲线图3. 一天总用水量的估计第1、2时段和第1、2供水时段流量的积分之和，就是一天总用水量虽然诸时段的流量已表为多项式函数，积分可以解析地算出，这里仍用数值积分计算：（1）第一用水时段的用水量其中积分值通过梯形公式计算 （2）第二用水时段的用水量 （3）第一供水时

15、段的用水量 (4). 第二供水和第三用水时段的用水量 y1=0.1*trapz(x1)； %第1时段用水量（仍按高度计），0.1为积分步长y2=0.1*trapz(x2)； %第2时段用水量y12=0.1*trapz(x12)； %第1供水时段用水量y3=0.1*trapz(x3)； %第2供水时段用水量y=(y1+y2+y12+y3)*237.8*0.01; %一天总用水量计算结果：y1=146.18, y2=266.59, y12=50.25, y3=70.00，y=1267.52具体程序见shuita5.m通过程序可以计算出各时段的用水量及一天的总用水量如下：各时段用水量第一用水时段用量

16、第二用水时段用量第一供水时段用量第二供水时段用量一天总用水量ViV1V2V12V3V347.62633.95119.49166.461267.52其中四个时段都用三次多项式进行拟合，如果第一用水时段、第一、二供水时段用三次多项式，第二用水时段用四次多项式拟合，得另一计算结果：y1=146.18 y2=266.75 y12=49.00 y3= 77.34 程序见shuita6.m各时段的用水量及一天的总用水量如下:各时段用水量第一用水时段用量第二用水时段用量第一供水时段用量第二供水时段用量一天总用水量ViV1V2V12V3V347.62613.76115.33187.021263.7 fig.5

17、一天内用水流量曲线图4. 流量及总用水量的检验与功率的计算 计算出的各时刻的流量可用水位记录的数值微分来检验用水量可用第1时段水位测量记录中下降高度来计算并检验，在第一用水时段水的实际用量为 类似地，用来计算和检验，在第二用水时段实际用水量为：而计算值为：=347.62, =613.767比较接近。供水时段流量的一种检验方法如下：供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等第1、2时段水泵的功率可计算如下：p1=(y12+260)/2； %第1供水时段水泵的功率（水量仍以高度计） tp2=20.8：

18、0.1：23；xp2=polyval（c3，tp4）； % xp2输出第2供水时段各时刻的流量p2=(0.1*trapz(xp2)+260)/2.2；%第2供水时段水泵的功率（水量仍以高度计）计算结果：p1=153.7 ，p2=140.1具体程序见附录shuita7.m.六：模型的分析与改进从计算结果可以看出，第1时段用水量与水位测量记录中的下降高度146相差无几，第2时段用水量与记录中的下降高度260相差无几，所以数据拟合、数值积分的精度是足够的。对不同次数的拟合多项式，第1、2供水时段用水量相差稍大，两供水时段水泵的功率也有差别，这都说明供水时段用3次曲线通过四点的做法不够好，应该多取几点

19、做拟合，但要注意让流量曲线在不同的时段相接处保持连续。下面给出拟合第1、2用水时段水位函数，采用不同次数的多项式所得流量及总用水量的结果，用n = (n1,n2)表示这两个时段所用多项式的次数，下表是各时段的用水量、一天总用水量及两个供水时段水泵的功率。(n1,n2)V1V2V12V3Vp1p2(3,4)347.62613.76115.33187.021263.7366.81366.05(5,6)348.41612.95108.79181.761251.9363.53363.66两个供水时段水泵的功率相当，说明该算法切实可行！另外可以得知各个时段的用水量及一天的总用水量，便于掌握用水情况。 f

20、ig.6 n=(3,4)流量曲线 fig.7 n=(5,6)流量曲线由上图可以看出，流量曲线与原始记录基本上那个相吻合，0点到10点流量很低，10点到下午3点是用水高峰，全天流量平均在22（）左右，若按这个平均流量计算，一天的总用水量应为与表二的结果很接近。七：模型的推广与应用通过对水塔流量的问题进行研究，掌握其工作和控制原理并对模型进行评价，此模型不仅适用于水塔的供水水位控制系统，而且还可以应用于其他的领域，像自动蓄水装置、水坝闸门控制蓄水系统。通过数据统计，建立模型，掌握水量的需求情况，合理控制闸门，使水量保持在一定的范围内，从而达到最优的效果。八：参考文献1 赵静 但琦，数学建模与数学实

21、验，北京：高等教育出版社，2000.112 陈杰，Matlab宝典，北京：电子工业出版社，2007.13 朱仁峰，精通7，北京清华大学学研大厦：北京大学出版社，2006.5九：附录程序Shuita1.mt=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 0

22、 0 1082 1050 1021 994 965 941 918 892 866 843 822 0 0 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3)a1=polyder(c1)tp1=0:0.1:9x1=-polyval(a1,tp1)plot(tp1,x1,'r.')程序shuita2.mt=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20

23、.84 22.01 22.96 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 0 0 1082 1050 1021 994 965 941 918 892 866 843 822 0 0 1059 1035 1018;c2=polyfit(t(13:23),h(13:23),4)a2=polyder(c2)tp2=10.9:0.1:21;x2=-polyval(a2,tp2);plot(tp2,x2,'b.')hold onxlabel('hour');ylabel('cm/hou

24、r');程序Shuita3.mt=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3)a1=polyder(c1);c2

25、=polyfit(t(10.9:21),h(10.9:21),4)a2=polyder(c2);xx1=-polyval(a1,8 9);xx2=-polyval(a2,11 12);xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3)a12=polyder(c12)tp12=9:0.02:11;x12=polyval(c12,tp12)plot(tp12,x12,'r.')hold onaxis(9 11 16 32)xlabel('hour')ylabel('cm/houor')程序shuita4.mt=0 0

26、.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c2=polyfit(t(11:21),h(11:21),4);a2=polyder(c2);dt3=diff(t(22:24);dh3=d

27、iff(h(22:24);dht3=-dh3./dt3;t3=20 20.8 t(22) t(23)xx3=-polyval(a2,t3(1:2),dht3;c3=-polyfit(t3,xx3,3);a3=polyder(c3)tp3=20.8:0.05:24;x3=-polyval(c3,tp3);plot(tp3,x3,'b.')hold on axis(20.5 24.5 18 27)xlabel('hour')ylabel('cm/houor')title('流量曲线图')程序shuita5.mt=0 0.92 1.84

28、 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3);a1=polyder(c1);tp1=0:0.1:9;x1=-polyval(a1,tp1);c

29、2=polyfit(t(11:21),h(11:21),3);a2=polyder(c2);tp2=10.9:0.1:21;x2=-polyval(a2,tp2);xx1=-polyval(a1,8 9);xx2=-polyval(a2,11 12);xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3);tp12=9:0.1:11;x12=polyval(c12,tp12);dt3=diff(t(22:24);dh3=diff(h(22:24);dht3=-dh3./dt3;t3=20 20.8 t(22) t(23);xx3=-polyval(a2,t3(1:

30、2),dht3;c3=-polyfit(t3,xx3,3);tp3=20.8:0.1:24;x3=-polyval(c3,tp3);y1=0.1*trapz(x1)y2=0.1*trapz(x2)y12=0.1*trapz(x12) y3=0.1*trapz(x3)y=(y1+y2+y12+y3)*237.8*0.01程序shuita6.mt=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.9

31、1;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3);a1=polyder(c1);tp1=0:0.1:9;x1=-polyval(a1,tp1);c2=polyfit(t(11:21),h(11:21),4);a2=polyder(c2);tp2=10.9:0.1:21;x2=-polyval(a2,tp2);xx1=-polyval(a1,8 9);xx2=-polyval(a2,11 12);xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3);tp12=9:0.1:11;x12=polyval(c12,tp12);dt3=diff(t(22:24);dh3=diff(h(22:24);dht3=-dh3./dt3;t3=20 20.8 t(22) t(23);xx3=-po