MATLAB--水塔流量的估计.doc

MATLAB--水塔流量的估计

《数值分析》课程设计论文
水塔水流量的估计

摘 要:数学建模方法是处理科学理论的一种经典方法，也是解决各类实际问题的常用方法。本文采用曲线拟合的方法，并利用数学软件MATLAB对水塔流量进行计算，计算结果与实际记录基本吻合。

关键词:建模，流量，拟合，MATLAB

1.问题重述
美国某州的各用水管理机构要求各社区提供用水率（以每小时多少加仑计，英制单位下，1加仑=4.54596dm3，美制单位下，1加仑=3.78533dm3）以及每天所用的总用水量，但许多社区并没有测量流入或流出当地水塔的水量的设备，而只能以每小时测量水塔的水位代替，其精度在0.5%以内。更为重要的是，无论什么时候，只要水塔中的水位下降到某一最低水位L时，水泵就启动向水塔重新充水直至某一最高水位H，但也无法得到水泵的供水量的测量数据。因此，在水泵正在工作时，不容易建立水塔中水位与水泵工作时用水量之间的关系。水泵每天向水塔充水一次或两次，每次大约2小时。试估计在任何时候，甚至包括水泵正在工作的时间内从水塔流出的流量，并估计一天的总用水量。水塔是一个垂直圆柱体，高为40英尺，直径为57英尺。
下表给出了某个小镇某一天的真实数据：
表1
某小镇某天的水塔水位（1m=3.281英尺）
时间（秒）
水位（英尺）
时间（秒）
水位（英尺）
时间（秒）
水位（英尺）
0
31.75
35932
水泵工作
68535
28.42
3316
31.10
39332
水泵工作
71854
27.67
6635
30.54
39435
35.50
75021
26.97
10619
29.94
43318
34.45
79154
水泵工作
13937
29.55
46636
33.50
82649
水泵工作
17921
28.92
49953
32.67
85968
34.75
21240
28.50
53936
31.56
89953
33.89
25223
27.87
57254
30.81
93270
33.40
28543
27.52
60574
30.12

32284
26.97
64554
29.27

2.问题分析

数据的单位转换：
表2
时间（h）
水位（m）
时间（h）
水位（m）
时间（h）
水位（m）
0
9.6769
9.98
水泵工作
19.04
8.6620
0.92
9.4788
10.93
水泵工作
19.96
8.4334
1.84
9.3081
10.95
10.8199
20.84
8.2201
2.95
9.1253
12.03
10.4998
22.01
水泵工作
3.87
8.9864
12.95
10.2103
22.96
水泵工作
4.98
8.8144
13.88
9.9573
23.88
10.5913
5.90
8.6864
14.98
9.6190
24.99
10.3292
7.01
8.5200
15.90
9.3904
25.91
10.1798
7.93
8.3877
16.83
9.1801

8.97
8.2201
17.94
8.9211

流量是单位时间流出的水的体积，由于水塔是正圆柱形，横截面积是常数，在水泵不工作的时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量拟合得到，作为用于拟合的原始数据，我们希望水泵不工作的时段流量越准确越好。这些流量大体可由两种方法计算: 一是直接对表2中的水位用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量。二是先用表中数据拟合水位-时间函数，求导数即可得到连续时间的流量。
一般说来数值微分的精度不高，何况测量记录还是不等距的，数值微分的计算尤其麻烦。下面我们用第二种方法处理。
有了任何时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，如表2可知从t=0到t=8.97(h)水位下降了9.6769-8.2201=1.4568 (m)，乘以水塔的截面积就是这一时段的用水量。这个数值可以用来检查拟合的结果。
3.模型假设
供水时段的假设
水泵第1次供水时段为t=9到t=11(h),第2次供水时段为t=20.8到t=23 (h)。这是根据最低和最高水位分别是8.2201m和10.8199m及表2的水位测量记录作出的假设。其中前3个时刻取自实测数据（精确到0.01h）,最后1个时刻来自每次供水约两小时的已知条件（从记录看，第2次供水时段应在有记录的22.96h之后不久结束）。水泵工作时单位时间的供水量基本为常数，这个常数大于单位时间的平均流量。流量是单位时间流出水的体积，这里假设流量是对时间的连续函数，即。为简化处理，不影响问题的解决，假设流量与水泵是工作无关。
由于水塔截面积是常数S ,为简单起见，计算中将流量定义为单位时间流出的水的高度，即水位对时间变化率的绝对值（水位是下降的），最后给出结果时再乘以S即可。即：
水位对时间的变化率(流量):
任何时刻的流量:
4.流量估计
4.1拟合水位-时间函数
从表2 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段）和3个水泵不工作时段(以下称第1用水时段t=0到t=8.97，第2用水时段t=10.95到t=20.48和第3用水时段t=23以后)。对第1、2用水时段的测量数据分别作多项式拟合，得到水位函数和。为使拟合曲线比较光滑，多项式次数不要太高，一般用3~6次。由于第3时段只有3个测量记录，无法对这一时段的水位作出比较好的拟合，可采用外推的办法解决。
4.2确定流量-时间函数
对于第1、2用水时段，只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3用水时段流量包含在第2供水时段内，需要拟合四个流量函数。
4.3一天的总用水量
总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。

5.算法设计与计算结果
5.1拟合第1、2时段的水位，并得出流量
1第1时段的流速
设t、h为已输入的时刻和水位测量记录，实现如下：
t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.94,19.04,19.96,20.84,23.88,24.99,25.91];
h=[ 9.6769,9.4788,9.3081,9.1253,8.9864,8.8144,8.6864,8.5200,8.3877,8.2201,10.8199,10.4998,10.2103,9.9573,9.6190,9.3904,9.1801,8.9211,8.6620,8.4334,8.2201,10.5913,10.3292,10.1798];
f1=polyfit(t(1:10),h(1:10),5); %用5次多项式拟合第1用水时段水位，f1输出5次多项式的系数
b1=polyder(f1); % b1输出多项式（系数为f1）导数的系数，给出水位变化率
tm1=0:0.01:8.97; %将第一用水时段[0 , 8.97]细分
g1=-polyval(b1,tm1); %g1输出多项式b1在tm1点的函数值(取负后边为正值)，即tm1时刻的流量(水位下降的速率)。
2第2时段的流速
实现如下：
f2=polyfit(t(11:21),h(11:21),5); %用5次多项式拟合第2用水时段水位水位，f2输出5次多项式的系数
b2=polyder(f2); %b2输出多项式（系数为f2）导数的系数,给出水位变化率
tm2=10.95:0.01:20.84; %将第二用水时段[10.95 , 20.84]细分
g2=-polyval(b2,tm2); %g2输出多项式（b2）在tm2点的函数值(取负后边为正值)，即tm2时刻的流量(水位下降的速率)
第1、2用水时段(水位变化率)曲线图:
5.2拟合供水时段的流量
1
在第1供水时段（t = 9~11）之前（即第1用水时段）和之后（第2用水时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量。为使流量函数在t =9 和t =11连续，只取4个点，拟合5次多项式（即曲线必过这4个点）。拟合5次多项式，实现如下：
q1=-polyval(b1,[7.93,8.97]); %取第1时段在t=7.93,8.97的流量
q2=-polyval(b2,[10.95,12.03]); %取第2时段在t=10.95,12.03的流量
dx=[7.93,8.97,10.95,12.03];
dy=[q1,q2]; %将四个点合并
d=polyfit(dx,dy,5) ; %拟合5次多项式
ex=8.97:0.01:10.95; %将第一供水时段[8.97,10.95]细分
ey=polyval(d,ex); %ey输出第一供水时段各时刻的流量
2
在第2供水时段之前取t =20,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=[19.96,20.84,t(22),t(23)]; %取第2供水时段前后各两点
a=[-polyval(b2,t3(1:2),dht3)]; %求得t3各时刻的流量
h=polyfit(t3,a,5); % 拟合5次多项式
mx=20.84:0.01:24.00; %将第2供水时段和第3用水时段细分
my=polyval(h,mx); % my输出第2供水时段（外推至t=24）各时刻的流量
第1供水时段，第2供水时段和第3用水时段(水位变化率)曲线图:

5.3一天的总用水量的估计
第1、2用水时段和第1、2供水时段流量的积分之和，就是一天总用水量。虽然诸时段的流量已表示为多项式函数，积分可以解析的算出，这里仍用数值积分计算。
（1）第1用水时间段的用水量：

其中积分值 h1 通过梯形公式计算：

计算得出
（2）第2用水时间段的用水量：

计算得出
（3）第1供水时间段的用水量：

计算得出
（4）第2供水和第3用水时段的用水量

计算得出

（5）一天总用水量

一天总用水量的曲线图如下：

6流量及总用水量的检验
6.1用水时段流量的检验方法：
计算出的各时刻的流量可用水位记录的数据来检验。用水量V1用第1用水时段水位测量记录中下降高度H1=9.6769-8.2201=1.4568m来计算并检验，在第一用水时段水的实际用量为: V1=S* H1
类似地，第二用水时段用去的水的高度H2=10.8199-8.2201=2.5998m,
实际用水量为:V2=S* H2
让v1、v2与V1、V2比较。
6.2供水时段流量的检验方法：
供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等。第1、2时段水泵的公率可以计算如下：
p1=[y12*S+(10.8199-8.2201)*S]/2.001 %第1供水时段水泵的功率(水量以高度计）
tm3=20.84:0.01:23.88;
g3=polyval(h,tm3); % g3输出第2供水时段各时刻的流量
p2=(0.01*trapz(my)+(10.8199-8.2201))*S/2.110 %第2供水时段水泵的功率（水量仍以高度计）
计算结果如下表：
表4
第1、2用水时段的（近似）总用水量
第1、2用水时段的（实际）总用水量
第1供水时段的水泵功率
第2供水时段的水泵功率
v1=363.0272
V1=363.0275
p1=376.9397
p2=376.6016
v2=647.9356
V2=647.8756
通过算法得出的第1、2用水时段的总用水量与实际的第1、2用水时段的总用水量非常接近，而两个供水时段水泵的功率相当，说明该算法切实可行！

7分析
由表2可以看出第1、2用水时间段的总用水量、水位下降的高度与实际的总用水量、水位下降高度相差无几。所以数据拟合，数值积分取得精度总够的。由一天总用水量的曲线图可以看出，流量曲线与原始记录基本上吻合，零点到10点钟流量最少，10点到下午3点是用水高峰期。

附录
部分源代码：
S=[((57/3.2)/2)^2]*pi;
t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.94,19.04,19.96,20.84,23.88,24.99,25.91];
h=[ 9.6769,9.4788,9.3081,9.1253,8.9864,8.8144,8.6864,8.5200,8.3877,8.2201,10.8199,10.4998,10.2103,9.9573,9.6190,9.3904,9.1801,8.9211,8.6620,8.4334,8.2201,10.5913,10.3292,10.1798];
f1=polyfit(t(1:10),h(1:10),5);
b1=polyder(f1);
tm1=0:0.01:8.97;
g1=-polyval(b1,tm1);
f2=polyfit(t(11:21),h(11:21),5);
b2=polyder(f2);
tm2=10.95:0.01:20.84;
g2=-polyval(b2,tm2);
q1=-polyval(b1,[7.93,8.97]);
q2=-polyval(b2,[10.95,12.03]);
dx=[7.93,8.97,10.95,12.03];
dy=[q1,q2];
d=polyfit(dx,dy,5) ;
ex=8.97:0.01:10.95;
ey=polyval(d,ex);
dt3=diff(t(22:24));
dh3=diff(h(22:24));
dht3=-dh3/dt3;
t3=[19.96,20.84,t(22),t(23)];
a=[-polyval(b2,t3(1:4),dht3)];
h=polyfit(t3,a,5);
mx=20.84:0.01:24.00;
my=polyval(h,mx);
y1=0.01*trapz(g1); %第1用水时段用水量(按高度计), 0.01为积分步长
y2=0.01*trapz(g2);
y12=0.01*trapz(ey);
y3=0.01*trapz(my);
v=(y1+y2+y12+y3)*S
v1=y1*S
v2=y2*S
V1=(9.6769-8.2201)*S
V2=(10.8199-8.2201)*S
p1=[y12*S+(10.8199-8.2201)*S]/2.001
tm3=20.84:0.01:23.88;
g3=polyval(h,tm3);
p2=(0.01*trapz(my)+(10.8199-8.2201))*S/2 .110
plot(tm1,g1,'r',ex,ey,'k',tm2,g2,'r',mx,my,'k')
13