水塔水流量估计模型

1、A题水塔水流量估计模型摘要文字略插值法算得各时段及全天的用水总量（以高度计）时段第1未供水 时段第2未供水 时段第1供水时 段混合时段全天Iglrcz 插值法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水泵开动859693475394353550899533397433183445

3、9327033403.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.

4、8512485.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.1648m0.3024）和10.7325m（35.50其匚訓期）（设出口的水位为零），而且、八烬飞沁厂1.1469,约为1,所以根据物理学的Torricelli定理：从小口流出液体的流速正比于水面高度的平方根，可忽略水位对流速的影响。2. 将流量看做时间的连续函数，为计算简单，不防将流量定义为单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），得到结果后再乘以水塔的横截面积S即可。3. 水塔横截面积 S=（57*0.3048）A2* 一一.TF +_。3.3流量估计方法首先根据表2的数据，用MATLA作出水位-时间散点图（图1）。程序见附录1水位

6、-时间散点.削下列计算流量与时间的关系。根据数据散点图1,一种简单的处理方法是将表2的数据分为三段，然后对每一段的数据做如下处理：设某段数据为（九旳,（乱1）,，（备Vn），相距数据中点的平均流速公式流速=（左端点的水位算得即-出每段数据首位点的流速采用下面的公式计算:v（ _）=（豔:“魁1土：細/（练鈿）， v（ .）=（也：妁点慎皿）/（音T*）.用以上公式算得流速与时间之间的数据表（表3）如下:表3流速与时间关系数据表时间（h）流速（cm/h）029.8913.4229.030.4621.7414.4326.361.3818.4815.4426.092.39516.2216.3724.

7、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.6325.9113.312.4931.52由表3作出流速一时间散点图（图2）。Matlab程序见附录2流速一时间共系散点團图2下面分别用数据插值法和数据拟合两种方法来估计水塔水流量。（1）数据插值法由表三，对水泵不

8、工作时段1和水泵不工作时段2采用插值方法，可以得到 任意时刻的流速，从而可以得到任意时刻的流量。对于水泵不工作时段1应用前 后时期的流速进行插值。由于最后一段水泵不工作时段数据太少，我们将它与水 泵工作时段2合并一同进行插值处理（该段以下简称混合时段）。这样，总共需 要对四段数据（第1, 2未供水时段，第1供水时段，混合时段）进行插值处理。A. 下面以第1未供水时段数据，分别用Iglrcz插值，分段线性插值，以及 三次样条插值三种方法算出流量函数和用水量（用水高度）。首先编写实现Iglrcz插值的函数文件，命名为Iglrcz.m.程序见附录3.用三种插值方法得出第1未供水时段流速插值曲线（见图

9、 3）。Matlab程序 见附录4.B. 下面以第2未供水时段数据，分别用lglrcz插值，分段线性插值，以及 三次样条插值三种方法算出流量函数和用水量（用水高度）。用三种插值方法得出第2未供水时段流速插值曲线（见图 4）。Matlab程序 见附录5.图4C. 下面以第1,2未供水时段和第1供水时段数据，分别用Iglrcz插值，分 段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。用三种插值方法得出第1供水时段流速插值曲线（见图 5）。Matlab程序见 附录6.图5D. 下面以第2未供水时段和混合时段数据，分别用Iglrcz插值，分段线性插值，以及三次样条插值三种方法算出

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

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

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

13、段之前取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相差不 大，说明差值结果与原始数据比较吻合。用三次样条插值法估计出全天的用水量约为578151X 如8X10 二 1374343皿（卅对于数据拟合法，由表5可得全天的用水总量约为528.82 X 2375X10二1257533,%侨，同数据插值法得到的结果也很接近。3.5模型总结本实验主要进行水塔水流量的估算，第一种估算方法为数据插值方法

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

16、read('建模.xlsx'I5流速与时间关系,'A2:A32');b=xlsread('建模.xlsx'I5流速与时间关系,'B2:B32');plot(a,b,'b+').xlabel('x轴');ylabel('y轴');J附录3.y=lglrcz(xO,yO,x)n=len gth(xO);m=le ngth(x);for i=1:mz=x(i);s=0.0;for k=1: np=1.0;for j=1: nif j=k p=p*(z-xO(j)/(xO(k)-xO(j);

17、endends=p*yO(k)+s;endy(i)=s;end附录4.ft=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.29t0=0:0.1:8.97; lgl=lglrcz(t,v,tO); lgljf=0.1*trapz(lgl)xx=i nterp1(t,v,tO);xxjf=0.1*trapz(xx)sy=in terp1(t,v,tO, syjf=0.1*trapz(sy) plot(t,v

18、,'igi'xx'spl in e');'r' ,tO,xx, 'g' ,tO,sy,'b')'sy'gtext(gtext(附录5.0.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,4;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;t0=10.95:0.1:20.84;lgl=lglrcz(t,v,t0);lgljf

19、=0.1*trapz(lgl)xx=i nterp1(t,v,t0);xxjf=0.1*trapz(xx)sy=in terp1(t,v,t0,'spli ne');syjf=0.1*trapz(sy)plot(t,v,'*' ,t0,lgl,'r',t0,xx,'g' ,t0,sy,'b')gtext('lgl')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

20、,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,23.86,22.17;t0=8.97:0.1:10.95;lgl=lglrcz(t,v,t0);lgljf=0.1*trapz(lgl)xx=i nterp1(t,v,t0);xxjf

21、=0.1*trapz(xx)sy=in terp1(t,v,t0,'spli ne' );syjf=0.1*trapz(sy)plot(t,v,'*' ,t0,lgl,'r' ,t0,xx,'g' ,t0,sy,'b' )gtext('lgl' )gtext('xx' )gtext('sy' )附录 7.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.

22、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=interp1(t,v,t0, 'spline' );syjf=0.1*trapz(sy)plot(t,v,'*' ,t0,lgl,'r&

23、#39; ,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,1 3.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.8

24、2,10.50,10.2; 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,1 3.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.2;c2=

25、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.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,1 3.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.2; c1=polyfit(

26、t(1:10),h(1:10),3);c2=polyfit(t(11:21),h(11:21),3);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);wgsz1=100*trapz(tp1,x113);%?u 1?' 1? e± ?卩?xuo ? d?x14=-polyval(a1,7.93,8.97);%?aoo ? ?3 i Dox ?± ?e y?Yx23=-polyv

27、al(a2,tp2);x114=-polyval(a2,10.95:0.01:20.84);wgsz2=100*trapz(tp2,x114);%?u 2?' 1? e± ?卩?xuo ? d?x24=-polyval(a2,10.95,12.03);%?aoo ?卩?3 i DoX ?± ?e y?Yx25=-polyval(a2,19.96,20.84);%?aoo ?卩?3 i DoX ?± ?e y?Ysubplot(1,2,1)plot(tp1,x13*100) subplot(1,2,2) plot(tp2,x23*100)附录11.t=0,0

28、.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,1 3.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.2; 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

29、.01:8.97; tp2=10.95:0.01:20.84;x13=-polyval(a1,tp1);x113=-polyval(a1,0:0.01:8.97);%?“ 1?' 1? e± ?卩?Xuo ? d?%?aoo ?卩?3 i DoX ?± ?e y?Y%?u 2?' 1? e± ?卩?Xuo ? d?%?aoo ?卩?3 i DoX ?± ?e y?Y%?aoo ?卩?3 i DoX ?± ?e y?Ywgsz1=100*trapz(tp1,x113);x14=-polyval(a1,7.93,8.97);x23=-polyval(a2,tp2);x114=-polyval(a2,10.95:0.01:20.84);wgsz2