用灰色理论分析的水箱水流量问题

水箱水流量问题
一、问题重述
许多供水单位由于没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位，估计在任意时刻t流出水箱的流量f(t)。

给出原始数据表以及单位换算之后（长度m，时间h）的表。

其中泵水用-1
二、模型假设
（1）影响水箱流量的唯一因素是该区公众对水的普通需要；
（2）水泵的灌水速度为常数；
（3）从水箱中流出水的最大流速小于水泵的灌水速度；
（4）每天的用水量分布都是相似的；
（5）水箱的流速可用光滑曲线来近似；
（6）当水箱的水容量达到3
⨯停止泵水。

51410g
677.610g
⨯开始泵水；达到3
三、符号说明
t：时间
V：水箱中水的体积
f（t）任意时刻t水箱的流量
四、模型建立与求解
（一）插值拟合解法：
1.首先将表格中的数据转化为标准单位制，其中时间用小时，体积用立方米。

1E=0.3024m
1L=7.481G
2.（1）水塔中水的体积计算：
24
V D h π
=
式中：D 为水塔的直径；h 为水塔中的水位高度。

（2）水塔中水流速度的估计。

水流速度应该是水塔中水的体积对时间的导数，但由于没有每一时刻水体积的具体数学表达式，这里我们用插商近似导数。

由于在两个时间段，水泵向水塔供水，无法确定水位的高度，因此计算水流速度要分3段进行计算。

第一段从0省道32284s ，第二段从39435s 到75021s ，第三段从85968s 到93270s
3.为了得到流速的连续曲线可以使用三次样条插值处理，然后做出时间-流速的散点图。

图1：时间流速散点图
4.曲线拟合
用MATLAB 进行三次样条插值拟合的曲线如下
图2：时间-流速曲线
5.对水流速度进行数值积分可以得到24小时的用水量为1358.43m
（二）GM （2,1）模型解法：
1.
-1
2.设时间序列为n t ，体积序列为n V 将11i i
i i
V V t t ++--组成新的一个数组n Q 并存储在文件名为data2.txt 的纯文本文件
中。

然后对序列n Q 建立GM （2,1）模型。

计算出
可以看出用GM（2,1）对本问题进行估计误差较大，所以不建议使用。

GM(2,1)模型流速图
其中蓝线是GM（2,1）模型给出的流速图。

五、附录
计算的MATLAB程序如下：
首先把原始数据粘贴到纯文本文件data.txt中，并把其中的泵水用-1标记。

clc,clear
a=load('data.txt');
t0=a(:,[1,3]);to=t0(:);%Ìá³öʱ¼äÊý¾Ý£¬²¢Õ¹¿ª³ÉÁÐÏòÁ¿
h0=a(:,[2,4]);h0=h0(:);%Ìá³ö¸ß¶ÈÊý¾Ý£¬²¢Õ¹¿ª³ÉÁÐÏòÁ¿
hs=0.3024;%µ¥Î»»»ËãÊý¾Ý
D=57*hs;%Ë®ËþÖ±¾¶£¬µ¥Î»m
h=h0/100*hs;%¸ß¶ÈÊý¾Ý£¬µ¥Î»m
t=t0/3600;%ʱ¼äת»¯³ÉСʱ
V=pi/4*D^2*h;%¼ÆËã¸÷¸öʱ¿ÌµÄÌå»ý
dv=gradient(V,t);%¼ÆËã¸÷ʱ¿ÌµÄÊýÖµµ¼Êý
no1=find(h0==-1);%ÕÒ³öÔ-ʼÎÞЧÊý¾ÝµÄµØÖ·
no2=[no1(1)-1:no1(2)+1,no1(3)-1:no1(4)+1];%ÕÒ³öµ¼ÊýÊý¾ÝÎÞЧµÄµØÖ·tt=t;tt(no2)=[];%ɾ³ýµ¼ÊýÊý¾ÝÎÞЧµÄµØÖ·µÄʱ¼ä
dv2=-dv;dv2(no2)=[];%¸ø³ö¸÷ʱ¿ÌµÄÁ÷ËÙ
plot(tt,dv2,'*')%»-³öÁ÷ËÙµÄÉ¢µãͼ
pp=csape(tt,dv2);%¶ÔÁ÷ËÙ½øÐвåÖµ
tt0=0:0.1:tt(end);%¸ø³ö²åÖµµã
fdv=ppval(pp,tt0);%¼ÆËã¸÷²åÖµµãµÄÁ÷ËÙÖµ
hold on,plot(tt0,fdv)%»-³ö²åÖµÇúÏß
I=trapz(tt0(1:241),fdv(1:241))%¼ÆËã24hÄÚ×ÜÁ÷Á¿µÄÊýÖµ»ý·Ö
clc,clear
a=load('data2.txt');%ÌáÈ¡Êý¾Ý
tt=a(:,1);
x0=a(:,2);
x0=x0/7.481
x0=x0'
n=length(x0);
x1=cumsum(x0)
a_x0=diff(x0)'
z=0.5*(x1(2:end)+x1(1:end-1))';
B=[-x0(2:end)',-z,ones(n-1,1)];
u=B\a_x0
x=dsolve('D2x+a1*Dx+a2*x=b','x(0)=c1,x(26)=c2');
x=subs(x,{'a1','a2','b','c1','c2'},{u(1),u(2),u(3),x1(1),x1(27)}); yuce=subs(x,'t',0:n-1);
x=vpa(x,6)
ezplot(x,[1,26])
x0_hat=[yuce(1),diff(yuce)]
epsilon=x0-x0_hat
plot(tt,x0_hat,tt,x0)
delta=abs(epsilon./x0)。