《水塔水流量的估計(jì)》課程設(shè)計(jì)論文

1、<p><b>　　水塔水流量的估計(jì)</b></p><p>　　摘 要:數(shù)學(xué)建模方法是處理科學(xué)理論的一種經(jīng)典方法，也是解決各類實(shí)際問題的常用方法。本文采用曲線擬合的方法，并利用數(shù)學(xué)軟件MATLAB對(duì)水塔流量進(jìn)行計(jì)算，計(jì)算結(jié)果與實(shí)際記錄基本吻合。</p><p>　　關(guān)鍵詞:建模，流量，擬合，MATLAB</p><p><b

2、>　　1.問題重述</b></p><p>　　美國某州的各用水管理機(jī)構(gòu)要求各社區(qū)提供用水率（以每小時(shí)多少加侖計(jì)，英制單位下，1加侖=4.54596dm3，美制單位下，1加侖=3.78533dm3）以及每天所用的總用水量，但許多社區(qū)并沒有測量流入或流出當(dāng)?shù)厮乃康脑O(shè)備，而只能以每小時(shí)測量水塔的水位代替，其精度在0.5%以內(nèi)。更為重要的是，無論什么時(shí)候，只要水塔中的水位下降到某一最低水位L時(shí)，

3、水泵就啟動(dòng)向水塔重新充水直至某一最高水位H，但也無法得到水泵的供水量的測量數(shù)據(jù)。因此，在水泵正在工作時(shí)，不容易建立水塔中水位與水泵工作時(shí)用水量之間的關(guān)系。水泵每天向水塔充水一次或兩次，每次大約2小時(shí)。試估計(jì)在任何時(shí)候，甚至包括水泵正在工作的時(shí)間內(nèi)從水塔流出的流量，并估計(jì)一天的總用水量。水塔是一個(gè)垂直圓柱體，高為40英尺，直徑為57英尺。</p><p>　　下表給出了某個(gè)小鎮(zhèn)某一天的真實(shí)數(shù)據(jù)：</p>

4、<p>　　表1：某小鎮(zhèn)某天的水塔水位（1m=3.281英尺）</p><p><b>　　2.問題分析</b></p><p><b>　　數(shù)據(jù)的單位轉(zhuǎn)換：</b></p><p><b>　　表2</b></p><p>　　流量是單位時(shí)間流出的水的體積，可以

5、由對(duì)應(yīng)時(shí)刻的流速乘以水塔的橫截面積得到。由于水塔是正圓柱形，橫截面積是常數(shù)S，所以我們在這里研究的其實(shí)是流速的變化。</p><p>　　在水泵不工作的時(shí)段，流量很容易從水位對(duì)時(shí)間的變化率,即流速算出，問題是如何估計(jì)水泵供水時(shí)段的流速。水泵供水時(shí)段的流速只能靠供水時(shí)段前后的流速擬合得到，作為用于擬合的原始數(shù)據(jù)，我們希望水泵不工作的時(shí)段流速越準(zhǔn)確越好。這些流速大體可由兩種方法計(jì)算: 一是直接對(duì)表2中的水位用數(shù)值微分

6、算出各時(shí)段的流速，用它們擬合其它時(shí)刻或連續(xù)時(shí)間的流速。二是先用表中數(shù)據(jù)擬合水位-時(shí)間函數(shù)，求導(dǎo)數(shù)即可得到連續(xù)時(shí)間的流速。</p><p>　　一般說來數(shù)值微分的精度不高，何況測量記錄還是不等距的，數(shù)值微分的計(jì)算尤其麻煩。下面我們用第二種方法處理。</p><p>　　有了任何時(shí)刻的流速，就不難計(jì)算一天的總用水量。其實(shí)，水泵不工作時(shí)段的用水量可以由測量記錄直接得到，如表2可知從t=0到t=8

7、.97(h)水位下降了，乘以水塔的橫截面積S就是這一時(shí)段的用水量。這個(gè)數(shù)值可以用來檢驗(yàn)擬合的結(jié)果。水泵第1次供水時(shí)段為t=8.98(h)到t=10.94(h),第2次供水時(shí)段為t=21(h)到t=23 (h)。這是根據(jù)最低和最高水位分別是8.2201m和10.8199m及表2的水位測量記錄作出的假設(shè)。其中前3個(gè)時(shí)刻取自實(shí)測數(shù)據(jù)（精確到0.01h）,最后1個(gè)時(shí)刻來自每次供水約兩小時(shí)的已知條件（從記錄看，第2次供水時(shí)段應(yīng)在有記錄的22.96

8、h之后不久結(jié)束）。</p><p><b>　　3.模型假設(shè)</b></p><p><b>　　供水時(shí)段的假設(shè)</b></p><p>　　水泵工作時(shí)單位時(shí)間的供水量基本為常數(shù)，這個(gè)常數(shù)大于單位時(shí)間的平均流量。流量是單位時(shí)間流出水的體積，這里假設(shè)是水位對(duì)時(shí)間的連續(xù)函數(shù)，即。為簡化處理，不影響問題的解決，假設(shè)流量與水泵是工

9、作無關(guān)。</p><p><b>　　4.流量估計(jì)</b></p><p>　　4.1擬合水位-時(shí)間函數(shù)</p><p>　　從表2 測量記錄看，一天有兩個(gè)供水時(shí)段（以下稱第1供水時(shí)段和第2供水時(shí)段）和3個(gè)水泵不工作時(shí)段(以下稱第1用水時(shí)段t=0到t=8.97，第2用水時(shí)段t=10.95到t=20.48和第3用水時(shí)段t=23以后)。對(duì)第1、2用

10、水時(shí)段的測量數(shù)據(jù)分別作多項(xiàng)式擬合，得到水位函數(shù)和。為使擬合曲線比較光滑，多項(xiàng)式次數(shù)不要太高，一般用3~6次。由于第3時(shí)段只有3個(gè)測量記錄，無法對(duì)這一時(shí)段的水位作出比較好的擬合，可采用外推的辦法解決。 </p><p>　　4.2確定流量-時(shí)間函數(shù)</p><p>　　對(duì)于第1、2用水時(shí)段，只需將水位函數(shù)求導(dǎo)數(shù)即可，對(duì)于兩個(gè)供水時(shí)段的流速，則用供水時(shí)段前后（水泵不工作時(shí)段）的流速擬合得到，并

11、且將擬合得到的第2供水時(shí)段流速外推，將第3用水時(shí)段流速包含在第2供水時(shí)段內(nèi)，需要擬合四個(gè)流速函數(shù)。</p><p>　　4.3一天的總用水量</p><p>　　總用水量等于兩個(gè)水泵不工作時(shí)段和兩個(gè)供水時(shí)段用水量之和，它們都可以由流速對(duì)時(shí)間的積分得到：</p><p>　　5.算法設(shè)計(jì)與計(jì)算結(jié)果</p><p>　　5.1擬合第1、2時(shí)段的水

12、位，并得出流量</p><p>　　5.1.1第1用水時(shí)段的流速</p><p>　　設(shè)t、h為已輸入的時(shí)刻和水位測量記錄，實(shí)現(xiàn)如下：</p><p>　　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

13、.96,20.84,23.88,24.99,25.91];</p><p>　　h=[ 9.6769,9.4788,9.3081,9.1253,9.0071,8.8144,8.6864,8.5030,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.329

14、2,10.1798];</p><p>　　f1=polyfit(t(1:10),h(1:10),5); %用5次多項(xiàng)式擬合第1用水時(shí)段水位，f1輸出5次多項(xiàng)式的系數(shù)</p><p>　　b1=polyder(f1); % b1輸出多項(xiàng)式（系數(shù)為f1）導(dǎo)數(shù)的系數(shù)，給出水位變化率</p><p>　　tm1=0:0.01:8.97; %將第一用水時(shí)段[0 , 8.

15、97]細(xì)分</p><p>　　g1=-polyval(b1,tm1); %g1輸出多項(xiàng)式b1在tm1點(diǎn)的函數(shù)值(取負(fù)后邊為正值)，即tm1時(shí)刻</p><p><b>　　的流速</b></p><p>　　5.1.2 第2用水時(shí)段的流速</p><p><b>　　實(shí)現(xiàn)如下：</b></

16、p><p>　　f2=polyfit(t(11:21),h(11:21),5); %用5次多項(xiàng)式擬合第2用水時(shí)段水位水位，f2輸出5次多項(xiàng)式的系數(shù)</p><p>　　b2=polyder(f2); %b2輸出多項(xiàng)式（系數(shù)為f2）導(dǎo)數(shù)的系數(shù),給出水位變化率</p><p>　　tm2=10.95:0.01:20.84; %將第二用水時(shí)段[10.95 , 20.84

17、]細(xì)分</p><p>　　g2=-polyval(b2,tm2); %g2輸出多項(xiàng)式（b2）在tm2點(diǎn)的函數(shù)值(取負(fù)后邊為正值)，即tm2時(shí)刻的流速</p><p>　　第1、2用水時(shí)間段的流速曲線圖:</p><p>　　5.2擬合供水時(shí)段的流速</p><p>　　5.2.1第1供水時(shí)間段的流速</p><p>

18、;　　在第1供水時(shí)段（t = 8.98~10.94）之前（即第1用水時(shí)段）和之后（第2用水時(shí)段）各取幾點(diǎn)，其流速已經(jīng)得到，用它們擬合第1供水時(shí)段的流速。為使流速函數(shù)在供水時(shí)段連續(xù)，只取4個(gè)點(diǎn)，擬合5次多項(xiàng)式（即曲線必過這4個(gè)點(diǎn)）。擬合5次多項(xiàng)式，實(shí)現(xiàn)如下：</p><p>　　q1=-polyval(b1,[7.93,8.97]); %取第1時(shí)段在t=7.93,8.97的流速</p><p&

19、gt;　　q2=-polyval(b2,[10.95,12.03]); %取第2時(shí)段在t=10.95,12.03的流速</p><p>　　dx=[7.93,8.97,10.95,12.03];</p><p>　　dy=[q1,q2]; %將四個(gè)點(diǎn)合并</p><p>　　d=polyfit(dx,dy,5) ; %擬合5次多項(xiàng)式</p>&l

20、t;p>　　ex=8.97:0.01:10.95; %將第一供水時(shí)段[8.97,10.95]細(xì)分</p><p>　　ey=polyval(d,ex); %ey輸出第一供水時(shí)段各時(shí)刻的流速</p><p>　　5.2.2第2供水及第3用水(到t=24)時(shí)間段的流速</p><p>　　在第2供水時(shí)段之前取t =20,20.8兩點(diǎn)的流速，在該時(shí)刻之后（第3用

21、水時(shí)段）僅有3個(gè)水位記錄，我們用差分得到流速，然后用這4個(gè)數(shù)值擬合第2供水時(shí)段的流速。實(shí)現(xiàn)如下： </p><p>　　dt3=diff(t(22:24)); %最后3個(gè)時(shí)刻的兩兩之差：</p><p>　　dh3=diff(h(22:24)); %最后3個(gè)水位的兩兩之差：</p><p>　　dht3=-dh3/dt3; % t(22)和t(23)的流速（差

22、商代替微商）</p><p>　　t3=[19.96,20.84,t(22),t(23)]; %取第2供水時(shí)段前后各兩點(diǎn)</p><p>　　a=[-polyval(b2,t3(1:4),dht3)]; %求得t3各時(shí)刻的流速</p><p>　　h=polyfit(t3,a,5); % 擬合5次多項(xiàng)式</p><p>　　mx=20.

23、84:0.01:24.00; %將第2供水時(shí)段和第3用水時(shí)段細(xì)分</p><p>　　my=polyval(h,mx); % my輸出第2供水時(shí)段（到t=24）各時(shí)刻的流速</p><p>　　第1供水時(shí)間段，第2供水及第3用水(到t=24)時(shí)間段的流速曲線圖:</p><p>　　5.3一天的總用水量的估計(jì)</p><p>　　第1、2

24、用水時(shí)段和第1、2供水時(shí)段流速的積分之和乘以S，就是一天總用水量。雖然諸時(shí)段的流量已表示為多項(xiàng)式函數(shù)，積分可以解析的算出，這里仍用數(shù)值積分計(jì)算。</p><p>　　5.3.1第1用水時(shí)間段的用水量：</p><p>　　其中積分值 h1通過梯形公式計(jì)算：</p><p>　　計(jì)算得出第1用水時(shí)間段的用水量：</p><p>　　5.3.2第

25、2用水時(shí)間段的用水量：</p><p>　　計(jì)算得出第2用水時(shí)間段的用水量：</p><p>　　5.3.3第1供水時(shí)間段的用水量：</p><p>　　計(jì)算得出第1供水時(shí)間段的用水量：</p><p>　　5.3.4第2供水和第3用水(到t=24)時(shí)間段的用水量</p><p>　　計(jì)算得出第2供水和第3用水(到t=

26、24)時(shí)間段的用水量</p><p><b>　　一天總用水量：</b></p><p>　　6.流速及總用水量的檢驗(yàn)</p><p>　　6.1用水時(shí)段流速的檢驗(yàn)方法：</p><p>　　計(jì)算出的各時(shí)刻的流速可用水位記錄的數(shù)據(jù)來檢驗(yàn)。用水量V1用第1用水時(shí)段水位測量記錄中下降高度來計(jì)算并檢驗(yàn)，在第一用水時(shí)段水的實(shí)際用

27、量為: </p><p>　　類似地，第二用水時(shí)段用去的水的高度</p><p><b>　　實(shí)際用水量為:，</b></p><p>　　而通過算法得到的近似值：</p><p>　　將、與、進(jìn)行比較，兩者相差無幾。</p><p>　　6.2供水時(shí)段流速的檢驗(yàn)方法：</p>&l

28、t;p>　　通過求水泵在兩個(gè)供水時(shí)間段的功率與,</p><p><b>　　,而兩個(gè)供水時(shí)段</b></p><p>　　水泵的功率應(yīng)大致相等。</p><p>　　第1、2時(shí)段水泵的公率計(jì)算如下：</p><p>　　p1=[y12*S+(10.8199-8.2201)*S]/2.039 %第1供水時(shí)段水泵的

29、功率(水量以高度計(jì)）</p><p>　　tm3=20.84:0.01:23.88; </p><p>　　g3=polyval(h,tm3); % g3輸出第2供水時(shí)段各時(shí)刻的流速</p><p>　　p2=(0.01*trapz(my)+(10.5913-8.2201))*S/2 %第2供水時(shí)段水泵的功率（水量仍以高度計(jì)）</p><p

30、>　　第1、2時(shí)段水泵的公率計(jì)算結(jié)果：</p><p>　　通過算法得出的第1、2用水時(shí)段的總用水量與實(shí)際的第1、2用水時(shí)段的總用水量非常接近，而兩個(gè)供水時(shí)段水泵的功率相當(dāng)，說明該算法切實(shí)可行！</p><p><b>　　7.分析</b></p><p>　　由表2可以看出第1、2用水時(shí)間段的總用水量、水位下降的高度與實(shí)際的總用水量、水

31、位下降高度相差無幾。所以數(shù)據(jù)擬合，數(shù)值積分取得精度總夠的。由水塔一天的流速曲線圖可以看出，流速曲線與原始記錄基本上吻合，零點(diǎn)到10點(diǎn)鐘用水量最少，10點(diǎn)到下午3點(diǎn)是用水高峰期。</p><p><b>　　參考文獻(xiàn)</b></p><p>　　[1] 李大潛 中國大學(xué)生數(shù)學(xué)建模競賽[M] 北京 高等教育出版社 1998</p><p>　　[2

32、] 葉其孝 大學(xué)生數(shù)學(xué)建模競賽輔導(dǎo)[M] 長沙 湖南教育出版社 1993</p><p>　　[3] 蕭樹鐵 數(shù)學(xué)實(shí)驗(yàn)[M] 北京 高等教育出版社 1999</p><p>　　[4] 樂經(jīng)良 數(shù)學(xué)實(shí)驗(yàn)[M] 北京 高等教育出版社 2005</p><p>　　[5] Justin matlab 數(shù)學(xué)實(shí)驗(yàn)</p><p><b>　

33、　附錄</b></p><p><b>　　部分源代碼：</b></p><p>　　S=[((57/3.2)/2)^2]*pi;</p><p>　　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

34、,17.94,19.04,19.96,20.84,23.88,24.99,25.91];</p><p>　　h=[ 9.6769,9.4788,9.3081,9.1253,9.0071,8.8144,8.6864,8.5030,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

35、,10.5913,10.3292,10.1798];</p><p>　　f1=polyfit(t(1:10),h(1:10),5); </p><p>　　b1=polyder(f1);</p><p>　　tm1=0:0.01:8.97;</p><p>　　g1=-polyval(b1,tm1); %第1用水時(shí)段</p>

36、<p>　　f2=polyfit(t(11:21),h(11:21),5);</p><p>　　b2=polyder(f2);</p><p>　　tm2=10.95:0.01:20.84;</p><p>　　g2=-polyval(b2,tm2); %第2用水時(shí)段</p><p>　　q1=-polyval(b1,[7.9

37、3,8.97]);</p><p>　　q2=-polyval(b2,[10.95,12.03]);</p><p>　　dx=[7.93,8.97,10.95,12.03];</p><p>　　dy=[q1,q2];</p><p>　　d=polyfit(dx,dy,5) ;</p><p>　　ex=8.97:

38、0.01:10.95;</p><p>　　ey=polyval(d,ex); %第1供水時(shí)段</p><p>　　dt3=diff(t(22:24));</p><p>　　dh3=diff(h(22:24));</p><p>　　dht3=-dh3/dt3;</p><p>　　t3=[19.96,20.84,

39、t(22),t(23)];</p><p>　　a=[-polyval(b2,t3(1:4),dht3)];</p><p>　　h=polyfit(t3,a,5);</p><p>　　mx=20.84:0.01:24.00;</p><p>　　my=polyval(h,mx); %第2供水及第3用水（到t=24）的時(shí)段</p&g

40、t;<p>　　y1=0.01*trapz(g1); </p><p>　　y2=0.01*trapz(g2);</p><p>　　y12=0.01*trapz(ey);</p><p>　　y3=0.01*trapz(my);</p><p>　　v=(y1+y2+y12+y3)*S %一天的總用水量</p>

41、<p><b>　　v1=y1*S</b></p><p><b>　　v2=y2*S</b></p><p>　　V1=(9.6769-8.2201)*S</p><p>　　V2=(10.8199-8.2201)*S</p><p>　　p1=[y12*S+(10.8199-8.2

42、201)*S]/2.039 </p><p>　　tm3=20.84:0.01:23.88;</p><p>　　g3=polyval(h,tm3);</p><p>　　p2=(0.01*trapz(my)+(10.5913-8.2201))*S/2 </p><p>　　A=-polyval(b1,t(1:10));</p>

43、<p>　　B=-polyval(b2,t(11 :21));</p><p>　　plot(t(1:10),A,’o’,tm1,g1’b’,dx,dy,’o’,ex,ey,’r’,t(11 :21),B,’o’,tm2,g2,’k’,t3,a,’o’,mx,my,’m’);</p><p>　　title(’水塔一天的流速曲線圖’);</p>

44、<p>　　legend(’第1用水時(shí)段該時(shí)刻對(duì)應(yīng)的流速’,’第1用水時(shí)段’，’第1供水時(shí)段該時(shí)刻對(duì)應(yīng)的流速’,’第1供水時(shí)段’,’第2用水時(shí)段該時(shí)刻對(duì)應(yīng)的流速’,’第2用水時(shí)段’,’第2供水及第3用水時(shí)段(到t=24)該時(shí)刻對(duì)應(yīng)的流速’,’第2供水及第3用水時(shí)段(到t=24)’);</p><p>　　xlabel(‘時(shí)間(h)’);</p><p>　　ylabel(‘流速