水塔水流量估计模型

1、A题 水塔水流量估计模型摘要文字略插值法算得各时段及全天的用水总量（以高度计）时段第1未供水时段第2未供水时段第1供水时段混合时段全天lglrcz插 值法145.6231258.866457.0597127.8244589.3736分段线性插值法147.143258.969749.6051115.0578570.7756三次样条插值法145.687258.654753.7915120.0178578.151数据拟合法算得各时段及全天的用水总量（以高度计）时段第1未供水时段第2未供水时段第1供水时段混合时段全天用水高度145.67260.6649.8272.67528.82一、问题的重述二、模型

2、的假设与符号说明三、模型的建立与求解实验问题 略初始数据：表1 某小镇某一天水塔水位统计表时间/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵开动35932水泵开动82649水泵开动39332水泵开动85969347539435355089953339743318344593270

3、33403.1问题分析 为了表示方便，我们将表1的数据全部化为国际标准单位（用excel计算），时间用小时（h）,高度用米（m）。表2 一天内水塔水位记录时间/(h)水位/(m)09.601212.9544444410.13040.9211111119.4046413.875833339.858241.8430555569.23529614.982222229.5770082.9497222229.05385615.903888899.3350883.8713888898.91172816.826111119.1082884.9780555568.74540817.931666678.8512

4、485.98.618419.03758.5942087.0063888898.4520819.959444448.3674087.9286111118.32204820.839166678.1557288.9677777788.15572822.015水泵开动9.981111111水泵开动22.95805556水泵开动10.92555556水泵开动23.8802777810.508410.9541666710.735224.9869444410.27252812.0327777810.4176825.9083333310.100163.2 模型假设1. 流量只取决于水位差，与水位本身无关。因为

5、水塔最低和最高水位分别是8.1648m(27×0.3024)和10.7325m(35.50×0.3024)(设出口的水位为零)，而且(10.7352/8.1648)1.1469,约为1，所以根据物理学的Torricelli定理：从小口流出液体的流速正比于水面高度的平方根，可忽略水位对流速的影响。2. 将流量看做时间的连续函数，为计算简单，不防将流量定义为单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），得到结果后再乘以水塔的横截面积S即可。3. 水塔横截面积S=(57*0.3048)2*4=237.8（平方米）。3.3 流量估计方法 首先根据表2的数据，用M

6、ATLAB作出水位-时间散点图（图1）。程序见附录1 图1下列计算流量与时间的关系。根据数据散点图1，一种简单的处理方法是将表2的数据分为三段，然后对每一段的数据做如下处理：设某段数据为(x0, y0),(x1,y1),···，(xn,yn),相距数据中点的平均流速公式 流速=（左端点的水位-右端点的水位）/区间长度算得即 vxi+1+xi2=yi-yi+1xi+1-xi每段数据首位点的流速采用下面的公式计算：v(x0)=(3y0-4y1+y2)/(x2-x0),v(xn)=(-3yn+4yn+1-yn-2)/(xn-xn-2).用以上公式算得流速与时间之间的数据

7、表（表3）如下： 表3 流速与时间关系数据表时间（h）流速（cm/h）029.8913.4229.030.4621.7414.4326.361.3818.4815.4426.092.39516.2216.3724.733.4116.317.3823.644.42515.3218.4923.425.4413.0419.5256.4515.4520.423.867.46513.9820.8422.178.4516.3522.02水泵开动8.9719.2922.96水泵开动9.98水泵开动23.8827.0910.93水泵开动24.4321.6210.9533.525.4518.4811.4929.

8、6325.9113.312.4931.52由表3作出流速时间散点图（图2）。Matlab程序见附录2图2下面分别用数据插值法和数据拟合两种方法来估计水塔水流量。（1）数据插值法由表三，对水泵不工作时段1和水泵不工作时段2采用插值方法，可以得到任意时刻的流速，从而可以得到任意时刻的流量。对于水泵不工作时段1应用前后时期的流速进行插值。由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段2合并一同进行插值处理（该段以下简称混合时段）。这样，总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理。A下面以第1未供水时段数据，分别用lglrcz插值，分段线性插值，以及三次样

9、条插值三种方法算出流量函数和用水量（用水高度）。首先编写实现lglrcz插值的函数文件，命名为lglrcz.m.程序见附录3.用三种插值方法得出第1未供水时段流速插值曲线（见图3）。Matlab程序见附录4. 图3B下面以第2未供水时段数据，分别用lglrcz插值，分段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。用三种插值方法得出第2未供水时段流速插值曲线（见图4）。Matlab程序见附录5. 图4C下面以第1,2未供水时段和第1供水时段数据，分别用lglrcz插值，分段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。用三种插值方法得出第1供水时段流

10、速插值曲线（见图5）。Matlab程序见附录6. 图5 D下面以第2未供水时段和混合时段数据，分别用lglrcz插值，分段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。用三种插值方法得出混合时段流速插值曲线（见图6）。Matlab程序见附录7.图6根据A，B，C，D四步的matlab的运算结果，用三种插值方法得到四段时间各自的用水量及全天的用水量。见表4表4插值法算得各时段及全天的用水总量（以高度计）时段第1未供水第2未供水第1供水段混合时段全天lglrcz插值法145.6231258.866457.0597127.8244589.3736分段线性插值法147.14325

11、8.969749.6051115.0578570.7756三次样条插值法145.687258.654753.7915120.0178578.151（2）数据拟合法 1）拟合水位时间函数 根据表2的数据，分别对地1,2未供水时段的测量数据直接作三次多项式拟合，得到水位与时间的关系函数。Matlab程序是对第1为供水时段的数据进行拟合,得到第1未供水时段水位与时间的三次多项式拟合曲线（图7），Matlab程序见附录8.图7Matlab程序对第2未供水时段的数据进行拟合,得到第2未供水时段水位与时间的三次多项式拟合曲线（图8），Matlab程序见附录9.图82） 确定流量时间函数 对第1,2未供水时

12、段的水位求导可得流量，用三次多项式拟合第1,2未供水时段的流速与时间关系曲线。 得第1,2未供水时段的流速与时间关系曲线（三次多项式拟合），见图9，matlab程序见附录10.图9上述作的拟合曲线与插值所得的曲线相比发现，用三次多项式拟合的效果不是很好，若改用5次多项式拟合则得效果很好的第1,2未供水时段的流速与时间关系曲线（图10），matlab程序见附录11.图10 第1供水时段的流速则用前后时刻的流速拟合得到。为使流速函数在t=9和t=11处连续，则只取4个点，用3次多项式拟合得第1供水时段流速与时间的关系曲线图（图11），同数据插值法的图5相应部分较为吻合。Matlab程序见附录12.

13、图11对混合时段，在第2供水时段之前取t=19.96,20.84两点的流速，用第三未供水时段的3个记录做差分得到两个流速数据21.62,18.48，然后用4个数据做3次多项式拟合得到混合时段的流速与时间的关系曲线（图12），同数据插值法的图6相应部分较为吻合。Matlab程序见附录13.图123）估计一天的用水总量 分别对供水的两个时段和不供水的两个时段积分（流速与时间的积分）并求和得到一天的用水总量约为528.82（此数据是用水总高度cm）。表5列出各段用水量，同数据插值法算得的表4数据相比，较为吻合。表5 数据拟合法算得各个时段及全天用水量（高度计）时段第1未供水时段第2未供水时段第1供水

14、时段混合时段全天用水高度145.67260.6649.8272.67528.824）流量与用水总量的检验 略3.4结果简单分析对数据插值法，由表4可以看出，使用三次样条插值法得到的第1,2未供水时段的用水高度结果145.6870和258.6547与表2中记录的下降高度146和260相差不大，说明差值结果与原始数据比较吻合。用三次样条插值法估计出全天的用水量约为578.151×237.8×10=1374843.08（升）对于数据拟合法，由表5可得全天的用水总量约为528.82×237.8×10=1257533.96（升），同数据插值法得到的结果也很接近。3

15、.5模型总结本实验主要进行水塔水流量的估算，第一种估算方法为数据插值方法，我们用了三种不同的插值法进行估算。第二种估算方法为数据拟合法，用多项式进行拟合，得到水塔水流量的估算。四、模型的评价五、模型的改进六、参考文献附录下面的程序均是由MATLAB软件编写的附录1. a=xlsread(建模.xlsx','单位换算后统计表','A2:A29');b=xlsread('建模.xlsx','单位换算后统计表','B2:B29');plot(a,b,'r*')xlabel('x轴'

16、);ylabel('y轴');附录2. a=xlsread('建模.xlsx','流速与时间关系 ','A2:A32');b=xlsread('建模.xlsx','流速与时间关系','B2:B32');plot(a,b,'b+')xlabel('x轴');ylabel('y轴');附录3function y=lglrcz(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0;

17、for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end附录4.t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97;v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29;t0=0:0.1:8.97;lgl=lglrcz(t,v,t0);lgljf=0.1*trapz(lgl)xx=interp1(t,v,t0);xxj

18、f=0.1*trapz(xx)sy=interp1(t,v,t0,'spline');syjf=0.1*trapz(sy)plot(t,v,'*',t0,lgl,'r',t0,xx,'g',t0,sy,'b')gtext('lgl')gtext('xx')gtext('sy')附录5. t=10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84;v=33.50,29.63,31.5

19、2,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17;t0=10.95:0.1:20.84;lgl=lglrcz(t,v,t0);lgljf=0.1*trapz(lgl)xx=interp1(t,v,t0);xxjf=0.1*trapz(xx)sy=interp1(t,v,t0,'spline');syjf=0.1*trapz(sy)plot(t,v,'*',t0,lgl,'r',t0,xx,'g',t0,sy,'b')gtext('lgl'

20、)gtext('xx')gtext('sy')附录6. t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84;v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00

21、,23.86,22.17;t0=8.97:0.1:10.95;lgl=lglrcz(t,v,t0);lgljf=0.1*trapz(lgl)xx=interp1(t,v,t0);xxjf=0.1*trapz(xx)sy=interp1(t,v,t0,'spline');syjf=0.1*trapz(sy)plot(t,v,'*',t0,lgl,'r',t0,xx,'g',t0,sy,'b')gtext('lgl')gtext('xx')gtext('sy')附录7.

22、t=10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84,23.88,24.43,25.45,25.91;v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17,27.09,21.62,18.48,13.30;t0=20.84:0.1:25.91;lgl=lglrcz(t,v,t0);lgljf=0.1*trapz(lgl)xx=interp1(t,v,t0);xxjf=0.1*trapz(xx)sy=interp

23、1(t,v,t0,'spline');syjf=0.1*trapz(sy)plot(t,v,'*',t0,lgl,'r',t0,xx,'g',t0,sy,'b')gtext('lgl')gtext('xx')gtext('sy')附录8. t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,2

24、3.88,24.99,25.66;h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),3);tp1=0:0.1:8.9;x1=polyval(c1,tp1);plot(tp1,x1)附录9. t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14

25、.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66;h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c2=polyfit(t(11:21),h(11:21),3);tp2=10.9:0.1:20.9;x2=polyval(c2,tp2);plot(tp2,x2)附录10. t=0,0.92,1.84,2.95,3.87,4.9

26、8,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66;h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),3);c2=polyfit(t(11:21),h(11:21),3);a1=p

27、olyder(c1);a2=polyder(c2);tp1=0:0.01:8.97;tp2=10.95:0.01:20.84;x13=-polyval(a1,tp1);x113=-polyval(a1,0:0.01:8.97);wgsz1=100*trapz(tp1,x113); %¼ÆËãµÚ1δ¹©Ë®Ê±¶ÎµÄ×ÜÓÃË®Á¿x

28、14=-polyval(a1,7.93,8.97); %ΪºóÃæµÄ³ÌÐò×¼±¸Êý¾Ýx23=-polyval(a2,tp2);x114=-polyval(a2,10.95:0.01:20.84);wgsz2=100*trapz(tp2,x114); %¼ÆËãµÚ2δ¹©Ë&

29、#174;ʱ¶ÎµÄ×ÜÓÃË®Á¿x24=-polyval(a2,10.95,12.03); %ΪºóÃæµÄ³ÌÐò×¼±¸Êý¾Ýx25=-polyval(a2,19.96,20.84); %ΪºóÃ

30、;æµÄ³ÌÐò×¼±¸Êý¾Ýsubplot(1,2,1)plot(tp1,x13*100)subplot(1,2,2)plot(tp2,x23*100)附录11. t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66;h=9.68,9

31、.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),5);c2=polyfit(t(11:21),h(11:21),5);a1=polyder(c1);a2=polyder(c2);tp1=0:0.01:8.97;tp2=10.95:0.01:20.84;x13=-polyval(a1,tp1);x113=-polyval(a1,0:0.01:8.97);

32、wgsz1=100*trapz(tp1,x113); %¼ÆËãµÚ1δ¹©Ë®Ê±¶ÎµÄ×ÜÓÃË®Á¿x14=-polyval(a1,7.93,8.97); %ΪºóÃæµÄ³ÌÐò×¼&#

33、177;¸Êý¾Ýx23=-polyval(a2,tp2);x114=-polyval(a2,10.95:0.01:20.84);wgsz2=100*trapz(tp2,x114); %¼ÆËãµÚ2δ¹©Ë®Ê±¶ÎµÄ×ÜÓÃË®Á¿x24=-polyval(a2,10.95,12.03); %ΪºóÃæµÄ³ÌÐò×¼±¸Êý