黃河小浪底調(diào)水調(diào)沙問(wèn)題

黃河小浪底調(diào)水調(diào)沙問(wèn)題摘要：本文利用插值擬合的方法通過(guò)Matlab工具模擬出了排沙量與時(shí)間、排沙量與水流量的函數(shù)關(guān)系，并且求出了總排沙量為1.704億噸。整個(gè)模型簡(jiǎn)單且方便計(jì)算，其中排沙量與水流量的函數(shù)關(guān)系為分段函數(shù)。關(guān)鍵詞：調(diào)水調(diào)沙Matlab插值擬合一、問(wèn)題重述2004年6月至7月黃河進(jìn)行了第三次調(diào)水調(diào)沙試驗(yàn)，特別是首次由小浪底、三門(mén)峽和萬(wàn)家寨三大水庫(kù)聯(lián)合調(diào)度，采用接力式防洪預(yù)泄放水，形成人造洪峰進(jìn)行調(diào)沙試驗(yàn)獲得成功．整個(gè)試驗(yàn)期為20多天，小浪底從6月19日開(kāi)始預(yù)泄放水，直到7月13日恢復(fù)正常供水結(jié)束．小浪底水利工程按設(shè)計(jì)攔沙量為75.5億立方米，在這之前，小浪底共積泥沙達(dá)14.15億噸．這次調(diào)水調(diào)試驗(yàn)一個(gè)重要目的就是由小浪底上游的三門(mén)峽和萬(wàn)家寨水庫(kù)泄洪，在小浪底形成人造洪峰，沖刷小浪底庫(kù)區(qū)沉積的泥沙．在小浪底水庫(kù)開(kāi)閘泄洪以后，從6月27日開(kāi)始三門(mén)峽水庫(kù)和萬(wàn)家寨水庫(kù)陸續(xù)開(kāi)閘放水，人造洪峰于29日先后到達(dá)小浪底，7表1:試驗(yàn)觀測(cè)數(shù)據(jù)單位：水流為立方米/秒，含沙量為公斤/立方米日期6.296.307.17.27.37.4時(shí)間8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10時(shí)間8:0020:008:0020:008:0020:008:0020:008:0020:008:0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085現(xiàn)在，根據(jù)試驗(yàn)數(shù)據(jù)建立數(shù)學(xué)模型研究下面的問(wèn)題：(1)給出估算任意時(shí)刻的排沙量及總排沙量的方法；(2)確定排沙量與水流量的變化關(guān)系。二、模型假設(shè)假設(shè)所給數(shù)據(jù)客觀準(zhǔn)確的反響了現(xiàn)實(shí)情況假設(shè)所給數(shù)據(jù)遵循一定規(guī)律變化，即是連續(xù)的假設(shè)模型中不需要考慮一些外在因素假設(shè)可將時(shí)間化為等分的時(shí)間點(diǎn)進(jìn)行計(jì)算三、符號(hào)說(shuō)明t:時(shí)間或時(shí)間點(diǎn)v:水流量S:含沙量V:排沙量四、問(wèn)題分析假設(shè)水流量和含沙量都是連續(xù)的，那么某一時(shí)刻的排沙量V=v(t)S(t)，其中v(t)為t時(shí)刻的水流量，而S(t)為t時(shí)刻的含沙量。通過(guò)觀察數(shù)據(jù)，這些數(shù)據(jù)是每個(gè)12小時(shí)采集一次，所以我們可以將時(shí)間設(shè)為時(shí)間點(diǎn)t，依次為1，2，3，……，24，單位時(shí)間為12h。為了找到排沙量與時(shí)間的關(guān)系，我們就要先找到水流量和含沙量與時(shí)間的關(guān)系，一但找到水流量和含沙量與時(shí)間的關(guān)系，那么所要求的問(wèn)題也就不難解決了。五、模型的建立與求解通過(guò)分析，我們假設(shè)水流量和含沙量都是連續(xù)的，那么我們開(kāi)始對(duì)問(wèn)題“(1)給出估算任意時(shí)刻的排沙量及總排沙量的方法〞進(jìn)行求解。我們通過(guò)Matlab工具將所知道的數(shù)據(jù)顯示為直觀的圖像，如下所示，具體程序見(jiàn)附錄的%tuxing.m。通過(guò)觀察圖像，我們可以看出其變化并不光滑，而且也沒(méi)有特定的表現(xiàn)出服從某種分布的趨勢(shì)。但是為了得到具體的計(jì)算函數(shù)，我們就必須對(duì)數(shù)據(jù)進(jìn)行擬合，所以通過(guò)Matlab先利用spline方法對(duì)數(shù)據(jù)進(jìn)行插值，從而提高精確度，使圖像變得光滑，然后利用多項(xiàng)式進(jìn)行擬合，當(dāng)多項(xiàng)式次數(shù)越高擬合也越準(zhǔn)確，但是由于數(shù)據(jù)受到的影響較多，所以這里的數(shù)據(jù)也不是準(zhǔn)確值，因此我們可以只取三次進(jìn)行擬合，也方便了后續(xù)的計(jì)算。于是我們分別對(duì)含沙量和水流量進(jìn)行插值擬合，便可以得到下面圖像和結(jié)果，具體程序見(jiàn)附錄%hansha.m和%liuliang.m。所得到的擬合函數(shù)為：y=0.014*x^{3}-1.3*x^{2}+21*x+16即含沙量與時(shí)間的關(guān)系式為：S=0.014*t^3-1.3*t^2+21*t+16所得到的擬合函數(shù)為：y=0.13*x^{3}-14*x^{2}+2.4e+002*x+1.5e+003即水流量與時(shí)間的關(guān)系式為：v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003因?yàn)槟骋粫r(shí)刻的排沙量V=v(t)S(t)，所以我們可以將所擬合出來(lái)的多項(xiàng)式帶入上式，通過(guò)Matlab進(jìn)行計(jì)算可以得到下面答案，程序見(jiàn)附錄%jisuan.m。ans=91/50000*t^6-73/200*t^5+2429/100*t^4-14573/25*t^3+2866*t^2+35340*t+24000即排沙量與時(shí)間的關(guān)系為：V=0.0018*t^6-0.365*t^5+24.29*t^4-582.92*t^3+2866*t^2+35340*t+24000由于這里的多項(xiàng)式次數(shù)過(guò)高，不便于計(jì)算和傳播，所以我們可以對(duì)其再進(jìn)行一次擬合，有下面結(jié)果，程序見(jiàn)附錄%paisha.m。所以擬合后的函數(shù)為V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004，通過(guò)圖像可以看出排沙量與時(shí)間服從正態(tài)分布，所以也可以化成的形式e的指數(shù)形式進(jìn)行擬合，這里就不再重復(fù)計(jì)算。我們得到了擬合函數(shù)，下面就可以計(jì)算出這幾天的總排沙量02412×60×60×95*an即總含沙量為1.704億噸。下面我們對(duì)問(wèn)題“(2)確定排沙量與水流量的變化關(guān)系。〞進(jìn)行分析計(jì)算。以下所有相關(guān)程序見(jiàn)附錄%paishui.m，下面就不重復(fù)說(shuō)明。我們先利用Matlab將排沙量和水流量的相關(guān)數(shù)據(jù)反映到圖像中。通過(guò)觀察可以看出，其關(guān)系是分段的，所以我們按時(shí)間進(jìn)行分段擬合，擬合原理同問(wèn)題〔1〕相同，于是可以得到分段前后的擬合多項(xiàng)式。y=-7.5e-005*x^{3}+0.43*x^{2}-5.2e+002*x+3.6e+004y=2.3e-005*x^{3}-0.066*x^{2}+1.9e+002*x-1.9e+005綜上，我們可以得到排沙量與水流量的關(guān)系式為-7.5e-5*v^3+0.43*v^2-5.2e+2*v+3.6e+40<=t<9V=2.3e-5*v^3-0.066*v^2+1.9e+2*v-1.9e+59<=t<=24六、模型評(píng)估本模型的優(yōu)點(diǎn)是：建模簡(jiǎn)單，方便計(jì)算，適用度廣。但也有最大的缺點(diǎn)為：精確度較低。為了減少誤差，我們可以通過(guò)增大模型中擬合多項(xiàng)式的次數(shù)。當(dāng)然在日后的模型改良中可以參加誤差評(píng)估系統(tǒng)，來(lái)對(duì)模型進(jìn)行完善。附錄%tuxing.mT=1:24;S=[32 60 75 85 90 98 100 102 108 112 115 116118 120 118 105 80 60 50 30 26 20 8 5];W=[1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 26502600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900];subplot(2,1,1);plot(T,S);holdon;plot(T,S,'.');title('時(shí)間與含沙量關(guān)系');xlabel('時(shí)間t/12h');ylabel('含沙量/公斤每立方米');subplot(2,1,2);plot(T,W);holdon;plot(T,W,'.');title('時(shí)間與水流量關(guān)系');xlabel('時(shí)間t/12h');ylabel('水流量/立方米每秒');%hansha.mT=1:24;S=[32 60 75 85 90 98 100 102 108 112 115 116118 120 118 105 80 60 50 30 26 20 8 5];x=1:0.1:24;y=interp1(T,S,x,'spline');plot(T,S,'.',x,y);title('時(shí)間與含沙量關(guān)系擬合圖');xlabel('時(shí)間t/12h');ylabel('含沙量/公斤每立方米');%liuliang.mT=1:24;W=[1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 26502600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900];x=1:0.1:24;y=interp1(T,W,x,'spline');plot(T,W,'.',x,y);title('時(shí)間與水流量關(guān)系擬合圖');xlabel('時(shí)間t/12h');ylabel('水流量/立方米每秒');%jisuan.msymst;S=0.014*t^3-1.3*t^2+21*t+16;v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003;V=v*S;simple(V);symst;V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004;int(12*60*60*V,t,0,24)%paisha.mt=1:24;V=0.0018*t.^6-0.365*t.^5+24.29*t.^4-582.92*t.^3+2866*t.^2+35340*t+24000;plot(t,V);title('時(shí)間與排沙量關(guān)系圖')%paishui.mt=1:24;v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;V=95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;plot(v,V,'.');title('整理圖')figure;t=1:9;v=0.13*t.^3-14*t.