MATLAB水塔流量的估计

水塔水流量的估计
摘要:数学建模方法是处理科学理论的一种经典方法，也是解决各类实际问题的常用方法。

本文采用曲线拟合的方法，并利用数学软件MATLAB对水塔流量进行计算，计算结果与实际记录基本吻合。

关键词:建模，流量，拟合，MATLAB
1.问题重述
美国某州的各用水管理机构要求各社区提供用水率（以每小时多少加仑计，英制单位下，1加仑=，美制单位下，1加仑=）以及每天所用的总用水量，但许多社区并没有测量流入或流出当地水塔的水量的设备，而只能以每小时测量水塔的水位代替，其精度在%以内。

更为重要的是，无论什么时候，只要水塔中的水位下降到某一最低水位L时，水泵就启动向水塔重新充水直至某一最高水位H，但也无法得到水泵的供水量的测量数据。

因此，在水泵正在工作时，不容易建立水塔中水位与水泵工作时用水量之间的关系。

水泵每天向水塔充水一次或两次，每次大约2小时。

试估计在任何时候，甚至包括水泵正在工作的时间内从水塔流出的流量()
f t，并估计一天的总用水量。

水塔是一个垂直圆柱体，高为40英尺，直径为57英尺。

下表给出了某个小镇某一天的真实数据：
表1
2.问题分析
数据的单位转换：
流量是单位时间流出的水的体积，由于水塔是正圆柱形，横截面积是常数，在水泵不工作的时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。

水泵供水时段的流量只能靠供水时段前后的流量拟合得到，作为用于拟合的原始数据，我们希望水泵不工作的时段流量越准确越好。

这些流量大体可由两种方法计算: 一是直接对表2中的水位用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量。

二是先用表中数据拟合水位-时间函数，求导数即可得到连续时间的流量。

一般说来数值微分的精度不高，何况测量记录还是不等距的，数值微分的计算尤其麻烦。

下面我们用第二种方法处理。

有了任何时刻的流量，就不难计算一天的总用水量。

其实，水泵不工作时段的用水量可以由测量记录直接得到，如表2可知从t=0到t=(h)水位下降了(m)，乘以水塔的截面积就是这一时段的用水量。

这个数值可以用来检查拟合的结果。

3.模型假设
供水时段的假设
水泵第1次供水时段为t=9到t=11(h),第2次供水时段为t=到t=23(h)。

这是根据最低和最高水位分别是和及表2的水位测量记录作出的假设。

其中前3个时刻取自实测数据（精确到）,最后1个时刻来自每次供水约两小时的已知条件（从记录看，第2次供水时段应在有记录的之后不久结束）。

水泵工作时单位时间的供水量基本为常数，这个常数大于单位时间的平均流量。

流量是单位时间流出水的体积，这里假设
流量是对时间的连续函数，即()t h h =。

为简化处理，不影响问题的解决，假设流量与水泵是工作无关。

由于水塔截面积是常数S ,为简单起见，计算中将流量定义为单位时间流出的水的高度，即水位对时间变化率的绝对值（水位是下降的），最后给出结果时再乘以S 即可。

即：
水位对时间的变化率(流量): 任何时刻的流量: ()()S t h t f ⋅'=
4.流量估计
拟合水位-时间函数
从表2 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段）和3个水泵不工作时段(以下称第1用水时段t=0到t=，第2用水时段t=到t=和第3用水时段t=23以后)。

对第1、2用水时段的测量数据分别作多项式拟合，得到水位函数()t h h 11=和()t h h 22=。

为使拟合曲线比较光滑，多项式次数不要太高，一般用3~6次。

由于第3时段只有3个测量记录，无法对这一时段的水位作出比较好的拟合，可采用外推的办法解决。

确定流量-时间函数
对于第1、2用水时段，只需将水位函数()2,1,==i t h h i i 求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3用水时段流量包含在第2供水时段内，需要拟合四个流量函数。

一天的总用水量
总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。

5.算法设计与计算结果
dt
t dh h )
(=
'⎰⎰'-='=t
t t t dt
h S dt V V 0
拟合第1、2时段的水位，并得出流量 1第1时段的流速
设t 、h 为已输入的时刻和水位测量记录，实现如下： t=[0,,,,,,,,,,,,,,,,,,,,,,,]; h=[ ,,,,,,,,,,,,,,,,,,,,,,,];
f1=polyfit(t(1:10),h(1:10),5); %用5次多项式拟合第1用水时段水位()t h h 11=，f1输出5次多项式的系数
b1=polyder(f1); % b1输出多项式（系数为f1）导数的系数，给出水位变化率 tm1=0::; %将第一用水时段[0 , ]细分
g1=-polyval(b1,tm1); %g1输出多项式b1在tm1点的函数值(取负后边为正值)，即tm1时刻的流量(水位下降的速率)。

2第2时段的流速 实现如下：
f2=polyfit(t(11:21),h(11:21),5); %用5次多项式拟合第2用水时段水位水位()t h h 22=，f2输出5次多项式的系数
b2=polyder(f2); %b2输出多项式（系数为f2）导数的系数,给出水位变化率 tm2=::; %将第二用水时段[ , ]细分
g2=-polyval(b2,tm2); %g2输出多项式（b2）在tm2点的函数值(取负后边为正值)，即tm2时刻的流量(水位下降的速率)
第1、2用水时段(水位变化率)曲线图:
拟合供水时段的流量
1
在第1供水时段（t = 9~11）之前（即第1用水时段）和之后（第2用水时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量。

为使流量函数在t =9 和t =11连续，只取4个点，拟合5次多项式（即曲线必过这4个点）。

拟合5次多项式，实现如下：
q1=-polyval(b1,[,]); %取第1时段在t=,的流量
q2=-polyval(b2,[,]); %取第2时段在t=,的流量
dx=[,,,];
dy=[q1,q2]; %将四个点合并
d=polyfit(dx,dy,5) ; %拟合5次多项式
ex=::; %将第一供水时段[,]细分
ey=polyval(d,ex); %ey输出第一供水时段各时刻的流量
2
在第2供水时段之前取t =20,两点的水流量，在该时刻之后（第3用水时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量。

实现如下：
dt3=diff(t(22:24)); %最后3个时刻的两两之差
dh3=diff(h(22:24)); %最后3个水位的两两之差
dht3=-dh3/dt3; % t(22)和t(23)的流量（差商代替微商）
t3=[,,t(22),t(23)]; %取第2供水时段前后各两点
a=[-polyval(b2,t3(1:2),dht3)]; %求得t3各时刻的流量
h=polyfit(t3,a,5); % 拟合5次多项式
mx=::; %将第2供水时段和第3用水时段细分
my=polyval(h,mx); % my输出第2供水时段（外推至t=24）各时刻的流量
第1供水时段，第2供水时段和第3用水时段(水位变化率)曲线图:
一天的总用水量的估计
第1、2用水时段和第1、2供水时段流量的积分之和，就是一天总用水量。

虽然诸时段的流量已表示为多项式函数，积分可以解析的算出，这里仍用数值积分计算。

（1）第1用水时间段的用水量：
其中积分值 h1 通过梯形公式计算：
计算得出0272.3631=v
（2）第2用水时间段的用水量：
计算得出9356.6472=v （3）第1供水时间段的用水量：
计算得出3987.10612=v
（4）第2供水和第3用水时段的用水量
计算得出7717.1463=v
（5）一天总用水量
111111
1010Sh dt h S dt h S dt V v t t t t t t ='='='=⎰⎰⎰∑∑⎰⎰=-='+'∆='='=-N
k k k N
k t t
t t t h t h t dt t h dt t h h k
k 1
1111
1112)
()()()(1
1
222
223
2
3232Sh dt h S dt h S dt V v t
t t t t t ='='='=⎰⎰⎰∑∑⎰⎰=-='+'∆='='=-N
k k k N
k t t
t t t h t h t dt t h dt t h h k
k 1
212
1
2222)()()()(1
3
2
12121212121
1
2121Sh dt h S dt h S dt V v t t t t t t ='='='=⎰⎰⎰∑∑⎰⎰=-='+'∆='='=-N
k k k N
k t t
t t t h t h t dt t h dt t h h k
k 1
12112
1
1212122)()()()(1
2
1
333
334
3
43
4
3Sh dt h S dt h S dt V v t t t t t t ='='='=⎰⎰⎰∑∑⎰⎰=-='+'∆='='=-N
k k k N
k t t t
t t h t h t dt t h dt t h h k
k 0
313
3332)()()()(1
=v
+
v
v
v
v
1.
+
3
1264
+
12
2
1=
一天总用水量的曲线图如下：
6流量及总用水量的检验
用水时段流量的检验方法：
计算出的各时刻的流量可用水位记录的数据来检验。

用水量V1用第1用水时段水位测量记录中下降高度H1=来计算并检验，在第一用水时段水的实际用量为: V1=S* H1 类似地，第二用水时段用去的水的高度H2=实际用水量为:V2=S* H2
让v1、v2与V1、V2比较。

供水时段流量的检验方法：
供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等。

第1、2时段水泵的公率可以计算如下：
p1=[y12*S+ %第1供水时段水泵的功率(水量以高度计）
tm3=::;
g3=polyval(h,tm3); % g3输出第2供水时段各时刻的流量
p2=*trapz(my)+ %第2供水时段水泵的功率（水量仍以高度计）
计算结果如下表：
表4
通过算法得出的第1、2用水时段的总用水量与实际的第1、2用水时段的总用水量非常接近，而两个供水时段水泵的功率相当，说明该算法切实可行！
7分析
由表2可以看出第1、2用水时间段的总用水量、水位下降的高度与实际的总用水量、水位下降高度相差无几。

所以数据拟合，数值积分取得精度总够的。

由一天总用水量的曲线图可以看出，流量曲线与原始记录基本上吻合，零点到10点钟流量最少，10点到下午3点是用水高峰期。

附录
部分源代码：
S=[((5710
plot(tm1,g1,'r',ex,ey,'k',tm2,g2,'r',mx,my,'k')。