水塔流量估計(jì)

《數(shù)學(xué)模型與數(shù)學(xué)軟件綜合訓(xùn)練》論文MATLABiheLnn^itugeojrechnimlCompurhsg^7訓(xùn)練題目：水塔Cap^ripfiL 門jcMdllLHbT：虬Eiic.學(xué)生學(xué)號：07500119姓名：周才祥計(jì)通院信息與計(jì)算科學(xué)專業(yè)指導(dǎo)教師：黃燦云（理學(xué)院）2010年春季學(xué)期在生產(chǎn)實(shí)踐和科學(xué)研究中，常常遇到這樣的問題：由實(shí)驗(yàn)或測量得到的一批離散樣點(diǎn)，需要確定滿足特定要求的曲線或曲面(即變量之間的函數(shù)關(guān)系或預(yù)測樣點(diǎn)之外的數(shù)據(jù))。如果要求曲線(面)通過所給的所有數(shù)據(jù)點(diǎn)(即確定一個(gè)初等函數(shù)通過已知各數(shù)據(jù)，一般用多項(xiàng)式或分段多項(xiàng)式)，這就是數(shù)據(jù)插值。在數(shù)據(jù)較少的情況下，這樣做能夠取得好的效果。但是，如果數(shù)據(jù)較多，那么插值函數(shù)是一個(gè)次數(shù)很高的函數(shù)，比較復(fù)雜。如果不要求曲線(面)通過所有的數(shù)據(jù)點(diǎn)，而是要求它反映對象整體的變化趨勢，可得到更簡單實(shí)用的近似函數(shù)，這就是數(shù)據(jù)擬合。函數(shù)插值和曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者在數(shù)學(xué)方法上是完全不同的。針對水塔數(shù)據(jù)分析，利用數(shù)學(xué)軟件MATLAB進(jìn)行數(shù)據(jù)擬合。曲線擬合問題是指：已知平面上n個(gè)點(diǎn)(x，y),i=0,1,…，n，x互不相同，尋求函數(shù)TOC\o"1-5"\h\zi i iy=f(x)，使f(x)在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。線性最小二乘法是解決曲線擬合最常用的方法，其基本思路是，令f(x)—ar(x)+ar(x)H 卜ar(x)1 1 2 2 mm其中r(x)是事先選定的一組函數(shù)，系數(shù)a(k—0，1，?，m,m<n)待定。尋求a，使得殘差k k k平方和Q—￡(f(x)-y)2達(dá)到最小。這里的建模原理實(shí)質(zhì)上與實(shí)驗(yàn)七中的回歸分析是一致的。i ii二1摘要數(shù)學(xué)建模方法是處理科學(xué)理論的一種經(jīng)典方法，也是解決各類實(shí)際問題的常用方法。文章采用曲線擬合的方法，并利用數(shù)學(xué)軟件MATLAB對水塔流蚤進(jìn)行計(jì)算計(jì)算結(jié)果與實(shí)際記錄基本吻合。關(guān)鍵詞：建模，流量，擬合，MATLAB3 07500119周才祥10年春數(shù)學(xué)模型與數(shù)學(xué)軟件綜合訓(xùn)練目錄錯誤！未找到引用源。2 07500119周才祥10年春數(shù)學(xué)模型與數(shù)學(xué)軟件綜合訓(xùn)練1問題提出某居民區(qū)有一供居民用水的圓柱形水塔，一般可以通過測量其水位來估計(jì)水的流量。面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最底水位時(shí)，水泵自動啟動向水塔供水，到設(shè)定的最高水位的時(shí)候停止供水，這段時(shí)間無法測量水塔的水位和水泵的供水量。通常水泵每天供水一兩次，每次約2h.水塔是一個(gè)高為12.2m,直徑為17.4m的正圓柱。按照設(shè)計(jì)。水塔水位降至約&2m時(shí)，水泵自動啟動，水位升至約10.8m時(shí)水泵停止工作。下表是某一天的水位測量記錄（符號“//”表示水泵啟動），試估計(jì)任何時(shí)刻（包括水泵正供水時(shí)）從水塔流出的水流量，及一天的總用水量。表格1時(shí)刻（h）00.921.842.953.874.985.907.017.938.97水位（cm）968948931913898881869852839822時(shí)刻（h）9.9810.9210.9512.0312.9513.8814.9815.9016.8317.94水位（cm）////108210501021994965941918892時(shí)刻（h）19.0419.9620.8422.0122.9623.8824.9925.91水位（cm）866843822////1059103510182問題分析流量是單位時(shí)間流出水的體積。由于水塔為正圓柱形，橫截面積是常數(shù)，在水泵不工作的時(shí)段,流量很容易從水位對時(shí)間的變化率算出，問題是如何估計(jì)水泵供水時(shí)的流量。水泵供水時(shí)段的流量只能靠供水時(shí)段前后的流量擬合得到。作為用于擬合的原始數(shù)據(jù)，希望水泵不工作時(shí)段的流量越準(zhǔn)確越好。大體有兩種計(jì)算方法：一是直接對表：1中的水位用數(shù)值微分計(jì)算出個(gè)時(shí)刻的流量，用它們擬合其他時(shí)刻或連續(xù)時(shí)間的流量；二是先用表中的數(shù)據(jù)擬合水位一時(shí)間函數(shù)，求導(dǎo)數(shù)可的連續(xù)時(shí)間的流量。一般來說數(shù)值為份的凈度不高，何況測量記錄不等距，結(jié)束值微分計(jì)算麻煩。因此選用中二種方法進(jìn)行處理。有了任何時(shí)刻的流量，就可以酸楚一天的總水量。其實(shí)，水泵不工作時(shí)段的用水量可根據(jù)記錄直接得出，有表：1可知t=0到t=8.97h,水位下降了968?822cm,乘以水塔的截面積即得該時(shí)段的用水量。該數(shù)值可用來檢驗(yàn)擬合的結(jié)果。3模型假設(shè)3?1忽略水位對流速的影響流量值取決于水位差，與水位本身無關(guān)，按Torriceli定律從小孔流出的流體的速度正比于高度的平方根，水塔的高度和最低水位為10.8m、802m（設(shè)出水口的水位為0），因故可忽略水位對流速的影響。3?2供水時(shí)段的假設(shè)水泵第一次供水時(shí)段為-=9h到t2=11h，第二次供水時(shí)段為t3=20.8h到t4=23h這是根據(jù)最高和最低水位分別為10.8m和8.2m，以及表1的水位測量紀(jì)錄作出的假設(shè)，其中前3個(gè)時(shí)刻直接取自實(shí)測數(shù)據(jù)（精確到0.1h），最后1個(gè)時(shí)刻來自每次供水約2h的已知條件（從記錄上看，第2次供水時(shí)段應(yīng)在有記錄的22.96h之后即將結(jié)束）。3?3單位時(shí)間的供水量為常數(shù)水泵工作時(shí)單位時(shí)間的供水量基本為常數(shù)，這個(gè)廠書讀與單位時(shí)間的平均流量。3.4流量是對時(shí)間的連續(xù)函數(shù)流量是單位時(shí)間流出水的體積，這里假設(shè)流量是對時(shí)間的連續(xù)函數(shù)。3?5流量與水泵是否工作無關(guān)為簡化處理，不影響問題的解決，假設(shè)流量與水泵是工作無關(guān)。3?6流量定義的假設(shè)由于水塔截面積是常數(shù)，為簡單起見，在計(jì)算中將流量定義做單位時(shí)間流出的水的高度，即水位（水位時(shí)下降的）對時(shí)間變化率的絕對值，最后給出結(jié)果再乘以s即可。4流量估計(jì)4.1擬合水位一時(shí)間函數(shù)由表1記錄，一天有兩個(gè)供水時(shí)段（以下稱第一個(gè)供水時(shí)段和第二個(gè)供水時(shí)段），和三個(gè)水泵不工作時(shí)段（以下稱第1時(shí)段t=0到t=8.97，第2個(gè)時(shí)段t=10.95到t-20.84，和第3時(shí)段t=23以后）。對第1、2時(shí)段的測量數(shù)據(jù)直接分別做多項(xiàng)式函數(shù)擬合，得到水位函數(shù)。為使擬合曲線比較光滑，多項(xiàng)式次數(shù)不要太高，一把為3~6次。由于第3時(shí)段只有3個(gè)測量記錄，無法對這一時(shí)段的水位作出較好的擬合。

4.2確定流量一時(shí)間函數(shù)對于第1、2時(shí)段只需將水位函數(shù)求導(dǎo)即可，對于兩個(gè)供水時(shí)段的流量，則用供水時(shí)段前后（水泵不工作）的流量擬合得到，且將第2時(shí)段流量外推，將第3個(gè)時(shí)段流量包含在第2供水時(shí)段內(nèi)。4.3一天總用水量的估計(jì)總用水量等于兩個(gè)不工作時(shí)段和兩個(gè)工作時(shí)段用水量之和，它們都可由流量對時(shí)間的積分得到。5算法設(shè)計(jì)與編程下面利用數(shù)學(xué)軟件MATLAB對水塔流量進(jìn)行計(jì)算與曲線擬合。5?1擬合第1、2時(shí)段的水位，并到處流量5.1?1第1時(shí)段的流速:實(shí)現(xiàn)代碼實(shí)現(xiàn)代碼k1*k1*、、、、、t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,9.98,10.92,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.94,19.04,19.96,20.84,22.01,22.96,23.88,24.99,25.91];h=[968 94893191389888186985283982200 1082 1050 1021 9949659419188928668438220 0 1059 1035 1018];%水泵啟動是水位用“0”表示c1=polyfit(t(1:10),h(1:10),3);%用3次多項(xiàng)式擬合第1時(shí)段水位，c1輸出3次多項(xiàng)式的系數(shù)a1=polyder(c1);%輸出多項(xiàng)式(系數(shù)為c1)導(dǎo)數(shù)的系數(shù)tp1=0:0.1:9;x1=-polyval(a1,tp1);%x1輸出多項(xiàng)式(系數(shù)為a1)在tp1點(diǎn)的函數(shù)值(取負(fù)后變?yōu)檎?，即tp1時(shí)刻的流量，類似可計(jì)算第2時(shí)段各時(shí)刻的流量k1*、、、、、、、、rE3zigure1.1 II回IIS3丨FileEditVievuInsertToolsDesktopWindovuHelpD^HSlkl^^C?? 口因h□E§:5.1(時(shí)間■水位散點(diǎn)團(tuán))980 1 1 1 1 1 1 1 1 TOC\o"1-5"\h\z960- -++++++++++940- -+920- -+900- + -880-+-+860--+840- + -8201 i ill ill i T■圖圖.2第2時(shí)段的流速實(shí)現(xiàn)代碼實(shí)現(xiàn)代碼c2=polyfit(t(l:10),h(l:10),3);%用3次多項(xiàng)式擬合第1時(shí)段水位，c2輸出3次多項(xiàng)式的系數(shù)a2=polyder(c2);%輸出多項(xiàng)式(系數(shù)為cl)導(dǎo)數(shù)的系數(shù)tp2=11:0.1:20;x2=-polyval(a2,tp2);%x2輸出多項(xiàng)式(系數(shù)為a2)在tp2點(diǎn)的函數(shù)值(取負(fù)后變?yōu)檎?，即tp2時(shí)刻的流量，類似可計(jì)算第2時(shí)段各時(shí)刻的流量5?2擬合供水時(shí)段的流量5.2.1在第一供水時(shí)段(t=9?11)之前(即第1時(shí)段)和之后(即第2時(shí)段)各取幾點(diǎn)，其流量已經(jīng)得到，用它們擬合第1供水時(shí)段的流量。為使流量函數(shù)在t=9和t=11連續(xù)，我們?nèi)我馊?個(gè)點(diǎn)，擬合3次多項(xiàng)式，實(shí)現(xiàn)如下：實(shí)現(xiàn)代碼實(shí)現(xiàn)代碼xx1=-polyval(a1,[8xx2=-polyval(a1,[119]);%取第1時(shí)段在t=8,9的流量xx1=-polyval(a1,[8xx2=-polyval(a1,[1112]);%取第2時(shí)段在t=11,12的流量xx12=[xx1 xx2];c12=polyfit([89 1112],xx12,3);%擬合3次多項(xiàng)式tp12=9:0.1:11;x12=polyval(c12,tp12);%x12輸出第一供水時(shí)段各時(shí)刻的流量。5.2.2在第2供水時(shí)段之前取t=20，20.8兩點(diǎn)的流量，在該時(shí)段之后(第3時(shí)段)僅有3個(gè)水位個(gè)水位記錄，用差分得到流量，然后用這4個(gè)數(shù)值擬合第二供水時(shí)段的流量如下:實(shí)現(xiàn)代碼實(shí)現(xiàn)代碼dt3=diff(t(22:24));%最后3個(gè)時(shí)刻的兩兩之差dh3=diff(h(22:24));%最后3個(gè)水位的兩兩之差dth3=-dh3/dt3;%t(22)和t(23)的流量t3=[20,20.8,t(22),t(23)];xx3=[-polyval(a2,t3(1:4),dth3)];%t3各時(shí)刻的流量c3=polyfit(t3,xx3,3);%擬合3次多項(xiàng)式tp3=20.8:0.1:24;x3=polyval(c3,tp3)；%x輸出第二供水時(shí)段(外推至t=24)各時(shí)刻的流量5?3—天的總用水量的估計(jì)第1、2時(shí)段和第一、二供水時(shí)段流量的積分之和，就是的總用水量。雖然諸時(shí)段流量已表示為多項(xiàng)式函數(shù)，積分可解析的算出，使用數(shù)值積分計(jì)算如下：實(shí)現(xiàn)代碼實(shí)現(xiàn)代碼y1=0.1*trapz(x1);%第1時(shí)段用水量(仍然以高度表示)，0.1為幾分步長y2=0.1*trapz(x2);%第2時(shí)段用水量yl2=0.1*trapz(xl2);%第一供水時(shí)段用水量y3=0.1*trapz(x3);%第二供水時(shí)段用水量y=(y1+y2+y12+y3)*237.8*0.01%一天的總用水量、、、、、、、、、、、、、、、、、、、、5.4流量及總用水量的檢驗(yàn)計(jì)算出各時(shí)刻的流量可用水位記錄的數(shù)值微分來檢驗(yàn)。用水量y1可用第1時(shí)段測量記錄中的下降高度968-822=146來檢驗(yàn)，類似的，y2用1082-822=260檢驗(yàn)。供水時(shí)段流量的一種檢驗(yàn)辦法如下：供水時(shí)段的用水量加上水位升高值260是該時(shí)段的抽入的水量，除以時(shí)段長度得到水泵的功率（單位時(shí)間抽入的水量），而兩個(gè)供水時(shí)段水泵的功率大致相等。第一二時(shí)段水泵的功率可計(jì)算如下：、!>、!>、!>、!>、!>、!*、、、、、、、、、、實(shí)現(xiàn)代碼k1*k1*、、、、、P1=（y12+260）/2;%第一供水時(shí)段的水泵的功率（水量以高度表示）tp2=20.8:0.1:23;xp2=polyval(c3,tp2);%xp2輸出第二供水時(shí)段各時(shí)刻的流量p2=(0.1*trapz(x3)+260)/2;%第二供水時(shí)段水泵的功率(水量以高度表示)k1*、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、6計(jì)算結(jié)果下面是水塔一天的流量曲線圖：圖6.1當(dāng)取三次多項(xiàng)式擬合的流量曲線圖由圖6.1我們可以看到，流量曲線與原始記錄基本上相吻合，但在第1時(shí)段和第1泵水時(shí)段的交接處曲線不太光滑，這說明我們采用3次曲線通過4點(diǎn)的做法不夠好，應(yīng)該多取幾點(diǎn)進(jìn)行擬合。0點(diǎn)到10點(diǎn)很流量很低，10點(diǎn)到下午3點(diǎn)即中午時(shí)間段是用水高峰期。7分析與改進(jìn)由表2可看出，第1時(shí)段用水量與水位測量記錄中的下降高度146相差無幾，第2時(shí)段用水量與記錄中的下降高度260相差無幾,所以數(shù)據(jù)擬合、數(shù)值積分的精度是足夠的。對不同次數(shù)的擬合多項(xiàng)式第1、2供水時(shí)段用水量相差稍大，兩供水時(shí)段水泵的功率亦有差別，都說明供水時(shí)段用3次曲線通過4點(diǎn)的做法不夠好，應(yīng)多取幾點(diǎn)作擬合，但要注意讓流量曲線在不同時(shí)段相接處保持連續(xù)。由圖1可看出，流量曲線與原始記錄基本上相吻合,零點(diǎn)到10點(diǎn)鐘流量很低,10點(diǎn)到下午3點(diǎn)是用水高峰，全天流量平均在22cm/h左右。若按這個(gè)平均流量計(jì)算，一天總用水量應(yīng)為2224237.80.01=1255.6L，與表2中的結(jié)果非常接近。參考文獻(xiàn)［1］ 李大潛中國大學(xué)生數(shù)學(xué)建棋競賽［M］北京高等教育出版社1998［2］ 葉其孝大學(xué)生數(shù)學(xué)建模競賽輔導(dǎo)［M］長沙湖南孩育出版社1993［3］ 蕭樹鐵數(shù)學(xué)實(shí)臉［M］北京高等教育出版社1999［4］ 樂經(jīng)良數(shù)學(xué)實(shí)臉［M］北京高等教育出版社2005［5］ Justinmatlab數(shù)學(xué)實(shí)驗(yàn) /view/3f06c6c4bb4cf7ec4afed033.html2010-7-4附錄I部分源代碼t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,9.98,10.92,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.94,19.04,19.96,20.84,22.01,22.96,23.88,24.99,25.91];h=[968 94893191389888186985283982200 1082 1050 1021 9949659419188928668438220 0 1059 1035 1018];%水泵啟動是水位用“0”表示c1=polyfit(t(1:10),h(1:10),3);%用3次多項(xiàng)式擬合第1時(shí)段水位，c1輸出3次多項(xiàng)式的系數(shù)a1=polyder(c1);%輸出多項(xiàng)式(系數(shù)為c1)導(dǎo)數(shù)的系數(shù)tp1=0:0.1:9;x1=-polyval(a1,tp1);%x1輸出多項(xiàng)式(系數(shù)為a1)在tp1點(diǎn)的函數(shù)值(取負(fù)后變?yōu)檎?，即tp1時(shí)刻的流量，類似可計(jì)算第2時(shí)段各時(shí)刻的流量c2=polyfit(t(1:10),h(1:10),3);%用3次多項(xiàng)式擬合第1時(shí)段水位，c2輸出3次多項(xiàng)式的系數(shù)a2=polyder(c2);%輸出多項(xiàng)式(系數(shù)為c1)導(dǎo)數(shù)的系數(shù)tp2=