水塔流量估計

水塔流量估計摘要1關(guān)鍵字2一.問題重述2二.問題分析2三.模型假設(shè)3四．問題求解3〔1〕數(shù)據(jù)差值法6〔2〕數(shù)值擬合法11五.結(jié)果簡單分析18六．模型檢驗18七．模型改良18八.實驗小結(jié)18九．參考文獻18附錄18摘要針對水塔流量估計問題，由于水塔流量與水位變化密切相關(guān)，所以通過水位變化來確定其流量；用流速和時間的關(guān)系來確定水位的變化，即水位差；由物理學(xué)中Torricelli定律確定流速，繼而解決具體問題。我們對某一天水塔水位統(tǒng)計數(shù)據(jù)進行分析，將其分為水泵供水時段和水泵不工作時段，然后通過插值與擬合模型分別進行計算。計算結(jié)果如下：圖1.1插值法算得各時段句全天的用水量〔用高度計〕第一未供水時段第二未供水時段第一供水時段混合時段全天拉格朗日插值法145.6231258.866454.0689127.8244586.3828分段線性插值法147.1430258.969749.6051115.0578570.7756三次樣條插值法145.6870258.654753.3334120.0178577.6929圖1.2數(shù)據(jù)擬合法算得各個時段及全天的用水量〔用高度計〕時段第一未供水時段第二未供水時段第一供水時段混合時段全天用水高度146.5150257.760546.131776.3076526.7148關(guān)鍵字：流量插值與擬合torricelli定律一.問題重述:某一地區(qū)的用水管理機構(gòu)要求各社區(qū)提供以每小時多少加侖計的用水量以及每天所用水的總量。許多社區(qū)沒有測量流入或流出水塔的水量裝置，他們只能通過測量水塔每小時的水位高度來代替，其誤差不超過0.05.但是，當(dāng)水塔中的水位下降到最低水位L時水泵就自動啟動向水塔輸水直到最高水位H，在此期間無法測量水泵的供水量。這樣，當(dāng)水泵正在輸水時就不容易建立水塔中水位和用水量之間的關(guān)系，水泵每天輸水一次或兩次，每次約兩小時。當(dāng)水位降至約27.00ft〔英尺〕時，水泵開始啟動向水塔供水，當(dāng)水位升到35.50ft時，水泵自動停止工作。某一小鎮(zhèn)的水塔高度為40ft，直徑為57ft的正圓柱，下表記錄的是某一天的水塔水位的真實數(shù)據(jù)〔注：1ft=0.3024m〕。表1.1某一天水塔水位統(tǒng)計表時間/s水位/0.01ft時間/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵開動35932水泵開動82649水泵開動39332水泵開動859683475394353550899533397433183445932703340試估計任何時刻〔包括水泵正在輸水時間〕從水塔流出的水流量和一天的用水總量。二.問題分析：流量是單位時間內(nèi)流出水的體積，由于水塔是正圓柱形，橫截面積時常數(shù)，所以在水泵不工作時段，流量很容易根據(jù)水位相對時間的變化率算出。問題的難點在于如何估計水泵供水時段的流量。水泵供水時段的流量只能依據(jù)供水時段前后的流量經(jīng)差值或擬合得到。作為用于差值或擬合的原始數(shù)據(jù)，我們希望水泵不工作時段的流量越準(zhǔn)確越好。事實上，水泵不工作時段的用水量可以由測量記錄直接得到嗎，由表中記錄的下降水位乘以水塔的橫截面積就該時段的用水量。這個數(shù)值可以用來檢驗數(shù)據(jù)差值或擬合的結(jié)果。為了表示方便，我們將表中的數(shù)據(jù)全部化為國際標(biāo)準(zhǔn)單位〔如下表〕，時間用小時〔h〕，高度用米〔m〕。表2.1一天內(nèi)水塔水位記錄時間〔h〕水位〔m〕時間〔h〕水位〔m〕09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.848.228.978.2222.02水泵開動9.98水泵開動22.96水泵開動10.93水泵開動23.8810.5910.9510.8224.9910.3512.0310.5025.9110.18三.模型假設(shè)〔1〕假設(shè)流量只取決于水位差，與水位本身無關(guān)。由水塔的最低和最高水位分別是8.1648m〔27x0.3024〕和10.7352m〔35.50x0.3024〕，設(shè)出口的水位為零，從小孔流出的流速正比于水面高度的平方根，可忽略水位對流速的影響?！?〕假設(shè)流量是時間的連續(xù)函數(shù)。為計算簡單，不妨將流量定義成單位時間流出水的高度，即水位對時間變化率的絕對值〔水位時下降的〕，得到結(jié)果后在乘以水塔橫截面積即可。〔3〕水塔橫截面積為S=〔57*0.3048〕^2*π/4=237.8〔平方米〕。四．問題求解首先根據(jù)表2.1的數(shù)據(jù)，用MATLAB軟件作出水位—時間散點圖4.1。程序見附錄1圖4.1水位—時間散點圖下面計算流量與時間的關(guān)系。根據(jù)數(shù)據(jù)散點圖1.0,一種簡單的處理方法是先將表2.0中的數(shù)據(jù)分為三段，然后對每一段數(shù)據(jù)做如下處理：設(shè)某段數(shù)句為{〔x0，y0〕，〔x1,y1〕,…,〔xn，yn〕}，相鄰數(shù)據(jù)中點的平均流速采用公式：流速=〔左端點的水位-右端點的水位〕/區(qū)間長度，即：vEQ〔〕=(yi-y(i+1))/(x〔i+1〕-xi)每段數(shù)據(jù)首尾點的柳樹采用下面的公式計算：v(x0)=(3y0-4y1+y2)/(x2-x0),v(xn)=(-3yn+4y(n-1)-y(n-2)).用以上公式算得流速與時間之間的數(shù)據(jù)表4.0如下：表4.1流速與時間關(guān)系數(shù)據(jù)表時間〔h〕流速〔cm/h〕時間〔h〕流速〔cm/h〕029.8913.4229.030.4621.7414.4326.361.3818.4815.4426.092.39516.2216.3724.733.4116.3017.3823.644.42515.3218.4923.425.4413.0419.5025.006.4515.4520.4023.867.46513.9820.8422.178.4516.3522.02水泵開動8.9719.2922.96水泵開動9.98水泵開動23.8827.0910.93水泵開動24.4321.6210.9533.5025.4518.4811.4929.6325.9113.3012.4931.52由表4.1作出流速—時間散點圖4.2，程序見附錄2：圖4.2流速—時間關(guān)系散點圖下面用數(shù)據(jù)差值和數(shù)據(jù)擬合兩種方法來估計水塔水流量：〔1〕數(shù)據(jù)差值法由表4.1，對水泵不工作時段1和水泵不工作時段2采用差值方法。可以得到任意時刻的流速，從而可以知道任意時刻的流量。對于水泵工作時段數(shù)據(jù)太少，我們將它與水泵工作時段2合并一同進行差值。由于最后一段水泵不工作時段數(shù)據(jù)太少，我們將它與水泵工作時段2合并一同進行差值處理〔該時段一下簡稱混合時段〕。這樣，總共需要對四段數(shù)據(jù)〔第1，2未供水時段，第1供水時段，混合時段〕進行差值處理。A.第1未供水時段數(shù)據(jù)差值求解，分別用拉格朗日插值法，分段線性差值及三次樣條差值三種方法算出流量函數(shù)和用水量〔用水高度〕。圖4.3中曲線lglr，fdxx，scyt分別表示用拉格朗日差值法，分段線性差值法及是三次樣條差值法對第1未供水時段數(shù)據(jù)差值得到的流速曲線，程序見附錄3圖4.3第1未供水時段流速曲線B.第2未供水時段數(shù)據(jù)差值求解，分別用拉格朗日插值法，分段線性差值及三次樣條差值三種方法算出流量函數(shù)和用水量〔用水高度〕。圖4.4中曲線lglr，fdxx，scyt分別表示用拉格朗日差值法，分段線性差值法及是那次樣條差值法對第2未供水時段數(shù)據(jù)差值得到的流速曲線，程序代碼見附錄4：圖4.4第2未供水時段流速曲線C.第1供水時段數(shù)據(jù)差值求解，分別用拉格朗日插值法，分段線性差值及三次樣條差值三種方法算出流量函數(shù)和用水量〔用水高度〕。圖4.5中曲線lglr，fdxx，scyt分別表示用拉格朗日差值法，分段線性差值法及是那次樣條差值法對第2未供水時段數(shù)據(jù)差值得到的流速曲線，程序代碼見附錄5：圖4.5第2未供水時段流速曲線D.混合時段數(shù)據(jù)差值求解，分別用拉格朗日插值法，分段線性差值及三次樣條差值三種方法算出流量函數(shù)和用水量〔用水高度〕。圖4.2中曲線lglr，fdxx，scyt分別表示用拉格朗日差值法，分段線性差值法及是那次樣條差值法對第2未供水時段數(shù)據(jù)差值得到的流速曲線:程序執(zhí)行后顯示結(jié)果，程序代碼見附錄6：圖4.6混合時段流速曲線綜合得到：表4.2差值法算得個時段及全天的用水總量〔以高度計〕時段及方法第1未供水時段第2未供水時段第1供水時段混合時段全天拉格朗日插值法145.6231258.866454.0689127.8244586.3828分段線性差值法147.1430258.969749.6051115.0578570.7756三次樣條差值法145.6870258.654753.3334120.0178577.6929圖4.6是用分段線性及三次樣條差值方法算得全天水塔水流速與時間的關(guān)系曲線，圖中曲線scyt表示三次樣條差值曲線，曲線fdxx表示分段線性差值的曲線，程序代碼見附錄7：圖4.7全天水塔流速與時間的關(guān)系曲線(2)數(shù)值擬合法1〕擬合水位—時間函數(shù)根據(jù)表4.1中的數(shù)據(jù)，分別對第1、2未供水時段的測量數(shù)據(jù)直接做三次多項式擬合，得到水位與實踐的關(guān)系函數(shù)。程序中，變量x1存放了以0.1為時間步長第1未供水時段各個時刻的水位高度。圖4.7繪出的是第1未供水時段水位與時間的三次多項式擬合曲線，程序代碼見附錄8：圖4.8第1未供水時段水位與時間的三次多項式擬合曲線下列圖4.9，Matlab程序?qū)Φ?未供水時段的數(shù)據(jù)進行擬合,得到的第2未供水時段水位與時間三次多項式擬合曲線，Matlab程序見附錄9.圖4.9第2未供水時段水位與時間的三次多項式擬合曲線2〕確定流量—時間函數(shù)對第1,2未供水時段的水位求導(dǎo)可得流量，用三次多項式擬合第1,2未供水時段的流速與時間關(guān)系曲線，程序代碼見附錄10：圖4.10第1,2未供水時段的流速與時間關(guān)系曲線〔三次多項式擬合〕與圖4.6中相應(yīng)局部段比照發(fā)現(xiàn)，用三次多項式擬合效果不是很好，假設(shè)改用5次多項式擬合那么得效果較好的第1,2未供水時段的流速與時間關(guān)系曲線，見圖4.9(程序代碼見附錄11):圖4.11第1,2未供水時段的流速與時間關(guān)系曲線〔五次多項式擬合〕第1供水時段的流速那么用前后時刻的流速擬合得到。為使流速函數(shù)在t=9和t=11處連續(xù)，只取4個點，用3次多項式擬合得第1供水時段流速與時間關(guān)系曲線圖4.12，同圖4.1中相應(yīng)局部段較為吻合，程序代碼見附錄12：圖4.12第一供水時段流速與時間的三次多項式似合曲線對于混合時段，在第2供水時段之前取t=19.96,20.84兩點的流速，用第三末供水時段的3個記錄做差分得到兩個流速數(shù)據(jù)21.62,18.48，然后用這4個數(shù)據(jù)做3次多項式似合得到混合時段的流速與時間的關(guān)系曲線如圖4.13所示，同數(shù)據(jù)插值法得到的曲線圖4.1中相應(yīng)曲線段較為吻合，程序代碼見附錄13：圖4.11混合時段流速與時間的三次多項式擬合曲線3〕估計一天的用水總量分別對供水的兩個時段和不供水的兩個時段積分〔流量對時間積分〕并求和得到一天的用水總量約為526.8935〔此數(shù)據(jù)是用水總高度cm〕.表列出各段用水量，同數(shù)據(jù)插值法算得的表4.1數(shù)據(jù)相比，較為吻合。表4.3數(shù)據(jù)擬合算得各個時段的用水高度〔一〕時段第1未供水時段第2未供水時段第1供水時段混合時段全天用水高度145.67260.6646.6073.9635526.89354〕流量及用水總量的檢驗計算出各時刻的流量可以用水位記錄數(shù)據(jù)的數(shù)值微分來檢驗，各時段的用水高度可以用實際記錄的水位下降高度來檢驗。例如，算得第1未供水時段的用水高度是145.67，而實際記錄的水位下降高度為968-822=146，兩者是吻合的同樣地，算得第2未供水時段的用水高度是260.66，而實際記錄的水位下降高度為1082-822=260，兩者也是吻合的。從算法設(shè)計和分析可知，計算結(jié)果與各時段所用的擬合多項式的次數(shù)有關(guān)。表4.4給出的是對第1，2未供水時段分別用5，6次多項式擬合得到的用水高度。表4.4數(shù)據(jù)擬合法算得各個時段的用水高度〔二〕時段第1未供水時段第2未供水時段第1供水時段混合時段全天用水高度146.5150257.760546.131776.3076526.7148五.結(jié)果簡單分析對于數(shù)據(jù)差值法，由表4.1可以看出，使用三次樣條差值法得到的第1,2未供水時段的用水高度結(jié)果145.6870和258.6547與表中記錄的下降高度146和260相差不大，說明差值結(jié)果與原始數(shù)據(jù)比擬吻合。用三次樣條差值法估計全天的用水量約為539.4450x237.8x10=1282800.21〔升〕。對于數(shù)據(jù)擬合法，由表可得全天的用水總量約為526.7148x237.8x10=1252527.7944〔升〕，同數(shù)據(jù)差值法得到的結(jié)果也很接近。六．模型檢驗七．模型改良八.實驗小結(jié)本次實驗主要進行水塔水量的估算。第一種估計為數(shù)據(jù)差值方法，我們用了三種不同的差值方法進行估計。在求解的過程中，可以熟悉數(shù)據(jù)差值的理論和方法。第二種估計方法為數(shù)據(jù)擬合法，用多項式進行擬合，得到水塔水流量的估計。九．參考文獻附錄：下面程序是對上面問題的補充，均由matlab軟件編寫：附錄1：x=xlsread(x=xlsread('stgj.xls','sheet1','A2:A29');y=xlsread('stgj.xls','sheet1','B2:B29');plot(x,y,'b*')xlabel('x軸');ylabel('y軸');附錄2：x=xlsread(x=xlsread('stgjls.xls','sheet1','A2:A32');y=xlsread('stgjls.xls','sheet1','B2:B32');plot(x,y,'b*')xlabel('x軸');ylabel('y軸');附錄補充：拉格朗日函數(shù)文件ffunctiony=lglrcz(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end附錄3：t=[0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97];t=[0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97];v=[29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.27];t0=0:0.1:8.97;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,’spline’);sancytjf=0.1*trapz(scyt)plot(t,v,’*’,t0,lglr,’r’,t0,fdxx,’g’,t0,scyt,’b’)gtext(‘lglr’)gtext(‘fdxx’)gtext(‘scyt’)附錄4t=[10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84];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];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];t0=10.95:0.1:20.84;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,’spline’);sancytjf=0.1*trapz(scyt)plot(t,v,’*’,t0,lglr,’r’,t0,fdxx,’g’,t0,scyt,’b’)gtext(‘lglr’)gtext(‘fdxx’)gtext(‘scyt’)附錄5：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];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;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,’spline’);sancytjf=0.1*trapz(scyt)plot(t,v,’*’,t0,lglr,’r’,t0,fdxx,’g’,t0,scyt,’b’)gtext(‘lglr’)gtext(‘fdxx’)gtext(‘scyt’)附錄6：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];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:22.96;lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,’spline’);sancytjf=0.1*trapz(scyt)plot(t,v,’*’,t0,lglr,’r’,t0,fdxx,’g’,t0,scyt,’b’)gtext(‘lglr’)gtext(‘fdxx’)gtext(‘scyt’)附錄7：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,23.88,24.43,25.45,25.91];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,27.09,21.62,18.48,13.30];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,23.88,24.43,25.45,25.91];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,27.09,21.62,18.48,13.30];t0=0:0.1:25.91;fdxx=interp1(t,v,t0);fdxxjf=trapz(fdxx)scyt=interp1(t,v,t0,'spline');sancytjf=trapz(scyt)plot(t,v,'*',t0,fdxx,'g',t0,scyt,'b')gtext('fdxx')gtext('scyt')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,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66]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,13.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.21,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);%輸入項3表示擬合多項式的次數(shù)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,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66]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,13.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.21,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);%輸入項5表示擬合多項式的次數(shù)tp1=0:0.1:8.9;x1=polyval(c1,tp1);plot(tp1,x1)附錄10：c1=polyfit(t(1:10),h(1:10),3);%用三次多項式擬合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]);wgsysl1=100*trapz(tp1,x113);%計算第1未供水時段的總用水量x14=-polyval(a1,[7.93,8.97]);%為后面的程序準(zhǔn)備數(shù)據(jù)x23=-polyval(a2,tp2);x114=-polyval(a2,[10.95:0.01:20.84]);wgsysl2=100*trapz(tp2,x114);%計算第2未供水時段的總用水量x24=-polyval(a2,[10.95,12.03]);%為后面的程序準(zhǔn)備數(shù)據(jù)x25=-polyval(a2,[19.96,20.84]);%為后面的程序準(zhǔn)備數(shù)據(jù)subplot(1,2,1)plot(tp1,x13*10