水塔水流量估计模型
A题水塔水流量估计模型
摘要
文字略
一、问题的重述
二、模型的假设与符号说明
三、模型的建立与求解
实验问题
略
初始数据：
3.1问题分析
为了表示方便，我们将表1的数据全部化为国际标准单位（用excel计算），时间用小时（h）,高度用米（m）。
3.2 模型假设
1. 流量只取决于水位差，与水位本身无关。因为水塔最低和最高水位分别是8.1648m0.3024)和10.7325m(35.50)(设出口的水位为零)，而且
1.1469,约为1，所以根据物理学的Torricelli定理：从
小口流出液体的流速正比于水面高度的平方根，可忽略水位对流速的影响。
2. 将流量看做时间的连续函数，为计算简单，不防将流量定义为单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），得到结果后再乘以水塔的横截面积S即可。
3. 水塔横截面积S=(57*0.3048)^2*。
3.3 流量估计方法
首先根据表2的数据，用MATLAB作出水位-时间散点图（图1）。程序见附录1
图1
下列计算流量与时间的关系。根据数据散点图1，一种简单的处理方法是将表2的数据分为三段，然后对每一段的数据做如下处理：设某段数据为
{(),(),···，()},相距数据中点的平均流速公式
流速=（左端点的水位）/区间长度
算得即
每段数据首位点的流速采用下面的公式计算：
v()=()/(),
v()=()/().
用以上公式算得流速与时间之间的数据表（表3）如下：
由表3作出流速—时间散点图（图2）。Matlab程序见附录2
图2
下面分别用数据插值法和数据拟合两种方法来估计水塔水流量。
（1）数据插值法
由表三，对水泵不工作时段1和水泵不工作时段2采用插值方法，可以得到任意时刻的流速，从而可以得到任意时刻的流量。对于水泵不工作时段1应用前后时期的流速进行插值。由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段2合并一同进行插值处理（该段以下简称混合时段）。这样，总共需
要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理。
A．下面以第1未供水时段数据，分别用lglrcz插值，分段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。
首先编写实现lglrcz插值的函数文件，命名为lglrcz.m.程序见附录3.
用三种插值方法得出第1未供水时段流速插值曲线（见图3）。Matlab程序见附录4.
图3
B．下面以第2未供水时段数据，分别用lglrcz插值，分段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。
用三种插值方法得出第2未供水时段流速插值曲线（见图4）。Matlab程序见附录5.
图4
C．下面以第1,2未供水时段和第1供水时段数据，分别用lglrcz插值，分段线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。
用三种插值方法得出第1供水时段流速插值曲线（见图5）。Matlab程序见附录6.
图5
D．下面以第2未供水时段和混合时段数据，分别用lglrcz插值，分段
线性插值，以及三次样条插值三种方法算出流量函数和用水量（用水高度）。
用三种插值方法得出混合时段流速插值曲线（见图6）。Matlab程序见附录7.
图6
根据A，B，C，D四步的matlab的运算结果，用三种插值方法得到四段时间各自的用水量及全天的用水量。见表4
表4
（2）数据拟合法
1）拟合水位—时间函数
根据表2的数据，分别对地1,2未供水时段的测量数据直接作三次多项式拟合，得到水位与时间的关系函数。Matlab程序是对第1为供水时段的数据进行拟合,得到第1未供水时段水位与时间的三次多项式拟合曲线（图7），Matlab程序见附录8.
图7
Matlab程序对第2未供水时段的数据进行拟合,得到第2未供水时段水位与时间的三次多项式拟合曲线（图8），Matlab程序见附录9.
图8
2）确定流量—时间函数
对第1,2未供水时段的水位求导可得流量，用三次多项式拟合第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.
图11
对混合时段，在第2供水时段之前取t=19.96,20.84两点的流速，用第三未供水时段的3个记录做差分得到两个流速数据21.62,18.48，然后用4个数据做3次多项式拟合得到混合时段的流速与时间的关系曲线（图12），同数据插值法的图6相应部分较为吻合。Matlab程序见附录13.
图12
3）估计一天的用水总量
分别对供水的两个时段和不供水的两个时段积分（流速与时间的积分）并求和得到一天的用水总量约为528.82（此数据是用水总高度cm）。表5列出各段用水量，同数据插值法算得的表4数据相比，较为吻合。
4）流量与用水总量的检验
略
3.4结果简单分析
对数据插值法，由表4可以看出，使用三次样条插值法得到的第1,2未供水时段的用水高度结果145.6870和258.6547与表2中记录的下降高度146和260相差不大，说明差值结果与原始数据比较吻合。
用三次样条插值法估计出全天的用水量约为578.151
对于数据拟合法，由表5可得全天的用水总量约为
，同数据插值法得到的结果也很接近。
3.5模型总结
本实验主要进行水塔水流量的估算，第一种估算方法为数据插值方法，我们用了三种不同的插值法进行估算。第二种估算方法为数据拟合法，用多项式进行拟合，得到水塔水流量的估算。
四、模型的评价
五、模型的改进
六、参考文献
附录
下面的程序均是由MATLAB软件编写的附录1.
附录3．
附录4.
附录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,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.
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=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')
附录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.82,10.50,10.2 1,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,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 1,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.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 1,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=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); %????μú1?′1???ê±??μ?×üó???á?
x14=-polyval(a1,[7.93,8.97]); %?aoó??μ?3ìDò×?±?êy?Y
x23=-polyval(a2,tp2);
x114=-polyval(a2,[10.95:0.01:20.84]);
wgsz2=100*trapz(tp2,x114); %????μú2?′1???ê±??μ?×üó???á?
x24=-polyval(a2,[10.95,12.03]); %?aoó??μ?3ìDò×?±?êy?Y
x25=-polyval(a2,[19.96,20.84]); %?aoó??μ?3ìDò×?±?êy?Y
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,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 1,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]);
wgsz1=100*trapz(tp1,x113); %????μú1?′1???ê±??μ?×üó???á?
x14=-polyval(a1,[7.93,8.97]); %?aoó??μ?3ìDò×?±?êy?Y
x23=-polyval(a2,tp2);
x114=-polyval(a2,[10.95:0.01:20.84]);
wgsz2=100*trapz(tp2,x114); %????μú2?′1???ê±??μ?×üó???á?
x24=-polyval(a2,[10.95,12.03]); %?aoó??μ?3ìDò×?±?êy?Y
x25=-polyval(a2,[19.96,20.84]); %?aoó??μ?3ìDò×?±?êy?Y
subplot(1,2,1)
plot(tp1,x13*100)
subplot(1,2,2)
plot(tp2,x23*100)
附录12．
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 1,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=polyder(c1);
a2=polyder(c2);
x14=-polyval(a1,[7.93,8.97]);
x24=-polyval(a2,[10.95,12.03]);
dygsdsj=[7.93,8.97,10.95,12.03];
dygsdls=[x14,x24];
nhjg=polyfit(dygsdsj,dygsdls,3)
nhsj=7.93:0.1:12.03;
nhlsjg=polyval(nhjg,nhsj)
gssjl=8.97:0.01:10.95;
gsl=polyval(nhjg,[8.97:0.01:10.95])
gsysl1=100.*trapz(gssjl,gsl) %计算第1供水时段的用水总量
plot(nhsj,100*nhlsjg)
附录13．
t3=[19.96,20.84,t(22),t(23)];
ls3=[x25*100,21.62,18.48];
nhhddxsxs=polyfit(t3,ls3,3)
tp3=19.96:0.01:25.91;
xx3=polyval(nhhddxsxs,tp3)
gssj2=20.84:0.01:24;
gs2=polyval(nhhddxsxs,[20.84:0.01:24]) gsysl2=trapz(gssj2,gs2) %混合时段用水总量plot(tp3,xx3)