B譚金梅B數(shù)學112班水庫排污問題正文

1、 水 庫 排 污 問 題 姓 名： 譚 金 梅 班 級： 數(shù) 學 112 學 號： 3110801242 指導老師： 周 金 明成 績： 完成日期： 2013年7月3日 摘 要數(shù)學建模是一種數(shù)學的思考方法，是運用數(shù)學的語言和方法，通過抽象、簡化建立能近似刻畫并解決實際問題的一種強有力的數(shù)學手段。數(shù)學建模就是用數(shù)學語言描述實際現(xiàn)象的過程。本文針對水庫排污問題，首先通過建立二維水質(zhì)污染物濃度模型，給出了單個水庫對干流造成大面積污染的可能性；然后建立兩水庫排污模型，分析了在另一水庫有連續(xù)點源污染物排放及水流相互影響的情況下，兩水庫對干流造成大面積污染的可能性；并進一步針對第三種情況的發(fā)生，給出在短時

2、間內(nèi)控制污染的有效措施；且討論了若污染物具有揮發(fā)性，上述各情況造成干流發(fā)生大面積污染的可能性大小，為水庫事故性排污問題提供了有價值的理論依據(jù)。問題一，只要考慮在事故發(fā)生到關(guān)閉水庫的兩個小時內(nèi)，流出水庫的污染物的質(zhì)量小于噸即可。問題二中，由于兩個水庫之間沒有聯(lián)系，只需要單獨考慮2庫的排污情況，然后加上1水庫的污染物排放情況，最后綜合考慮兩個水庫所排污染物的總量對干流的影響。問題三中建立人工水渠就是在問題二的基礎(chǔ)上使水庫1和水庫2發(fā)生聯(lián)系。由于水都是從高水位流向低水位，因此，只需考慮從1水庫向2水庫的流入情況，而不必考慮從2水庫向1水庫的流入情況。問題四，通過稀釋原理，建立污染物在干流中遷移遷移模

3、型，使干流水體計算體積元內(nèi)該污染物的增量為負值時，從而使得在干流水體污染物的濃度低于危險警戒值時的濃度，才能在短時間內(nèi)達到控制污染。關(guān)鍵字：水流速度、污染物濃度、流量、水庫排污、MATLAB問題重述某條江流上有2條支流，每條支流上都興建了規(guī)模相當?shù)乃畮?。由于正處在雨水多發(fā)季節(jié)，因此兩個水庫都以一定規(guī)模的流量進行泄洪。某天晚上10：00，在其中的一個水庫中發(fā)生了兩船相撞的事故，而其中的一條船裝載的p噸化學物質(zhì)（這里的化學物質(zhì)可以是具有揮發(fā)性的，也可能是急難揮發(fā)的）全部泄漏至水庫中。當水上航運事故處置中心接獲事故報告，立即要求該水庫關(guān)閉水庫泄洪閘，以免化學物質(zhì)隨洪水流入干流，發(fā)生更大規(guī)模的污染。水

4、庫閘門開始關(guān)閉時，已經(jīng)處在事故發(fā)生后的1個小時，而水庫閘門徹底關(guān)閉也需要1個小時的時間。 根據(jù)當?shù)丨h(huán)境監(jiān)測的有關(guān)規(guī)定，干流大面積污染的危險警戒值設(shè)為：三小時內(nèi)q噸該化學物質(zhì)發(fā)生泄漏。(1) 試建立合理的數(shù)學模型，討論由于此次事故的發(fā)生，干流發(fā)生大面積污染的可能性；(2) 如果在另外的一水庫中有一化工廠違規(guī)排放廢料。廢料中同樣含有該化學物質(zhì)。該工廠為躲避環(huán)境監(jiān)測站的監(jiān)控，均在晚上9：00-12：00違規(guī)進行周期性排放。在這種情形下，討論由于此次事故的發(fā)生，干流發(fā)生大面積污染的可能性；(3) 如果以上兩個水庫間有一條人工修建的水渠相連接，水渠中的水流流向不定，但保證兩水庫之間的水流能夠相互影響。那

5、么上述結(jié)果是否會改變？請給出說明，若有改變，則給出修正的模型及結(jié)果； (4) 如果發(fā)生了大面積污染，那么針對第三種情況，試給出在短時間內(nèi)控制污染模型。二、問題分析由題目知道，事故發(fā)生時兩水庫都正在泄洪，因此此時水庫中的水流速水庫較快。而泄露到水庫中的化學物質(zhì)不論是具有揮發(fā)性的，還是急難揮發(fā)的，它們對干流污染的情況總是類同的，因此我們總可以認為污染物是易溶急難揮發(fā)性物質(zhì)。為使我們的模型簡單，我們可以先假設(shè)事故發(fā)生在水庫1中，污染物在水庫中的分布是符合零維遷移模型的，此時流入水庫的污染物能以很快的速度與水庫中的水均勻混合，水庫中任何部份水體的污染狀況都是一樣的，污染程度與水體在水庫中的位置無關(guān)。而

6、實際上，污染物在水中達到分布均勻是有一個遷移過程，符合污染物一維遷移方程，但在這樣一個突發(fā)事件要求短時間內(nèi)得到控制的問題中，我們總可以用污染物瞬時混合均勻狀態(tài)模型來代替污染物一維遷移的過程。在解答問題一時，我們只要考慮在事故發(fā)生到關(guān)閉水庫的兩個小時內(nèi)，流出水庫的污染物的質(zhì)量小于噸即可。問題二中提到的情況只是將干流的污染源從一個增加到兩個，那么發(fā)生大面積污染的可能性就要增大。由于兩個水庫之間沒有聯(lián)系，我們只需要單獨考慮2庫的排污情況，然后加上1水庫的污染物排放情況，最后綜合考慮兩個水庫所排污染物的總量對干流的影響就可以了。問題三中建立人工水渠就是在問題二的基礎(chǔ)上使水庫1和水庫2發(fā)生聯(lián)系。由于水都

7、是從高水位流向低水位，因此，我們只需考慮從1水庫向2水庫的流入情況，而不必考慮從2水庫向1水庫的流入情況。 分析問題四，干流已經(jīng)發(fā)生大面積污染，對比水體污染物處理的各種手段，不論是化學、物理還是生物手段都不可能在短時間內(nèi)除去污染物，因此我們只能通過稀釋原理，建立污染物在干流中遷移遷移模型，使干流水體計算體積元內(nèi)該污染物的增量為負值時，從而使得在干流水體污染物的濃度低于危險警戒值時的濃度，才能在短時間內(nèi)達到控制污染。三、模型假設(shè)污染物為速溶物質(zhì)，因此藥品從船上流入水中的時間很少，可以忽略不計；污染物質(zhì)從水庫中一經(jīng)流出就進入干流；水庫和河流中的水流都是處于推流狀態(tài)；兩水庫事故發(fā)生條件相同，即兩水庫

8、有相同的客觀條件；被污染的水庫關(guān)閉泄洪閘后不再有水流流入干流；不考慮生物等因素在水庫泄洪過程中的作用，污染物除了流出外不因腐爛沉積等手段從水中消失；外界因素不對水庫的體積變化產(chǎn)生影響，例如：雨水、地表徑流、底下徑流等；參與模型的變量是連續(xù)變化的，并且充分光滑；不考慮從不同的渠道流入與流出水庫之間的區(qū)別，只考慮攜帶污染物的水流入水庫和水庫中的水流出對水庫污染程度的影響，因此可以把水庫看成是單流入單流出的系統(tǒng)。符號說明 （1）：t時刻水庫水的流入速度；（2）：t時刻流入水庫的污染物的濃度；（3）：t時刻水庫水的流出的速度；（4）：t時刻流出水庫的污染物的濃度；（5）：t時刻水庫中污染物的濃度；（6

9、）：t時刻水庫水的體積；（7）：計算體積元內(nèi)該污染物的增量；（8） ： 為時間；（9） ： 為從水庫中流出的水中的污染物的濃度；（10） ：為水庫的流出速度，即流量；（11）：為泄洪閘處到污染處的距離；（12） Q：為水庫的流量；（13） q：為排入河流的污水的流量；（14）：為河流中污染物的本底濃度；（15）:為水庫中的污染物的濃度； 注意：部分符號見論文中說明。五、模型的建立與求解5.1 模型一由問題的分析中知道，流入水庫的污染物能以很快的速度與水庫中的水均勻混合，也就是說水庫中的污染狀況在任何局部水體都是一樣的，污染程度與水體在水庫中的位置無關(guān)，因此我們可以建立下面模型。問題一、此次事故

10、的發(fā)生，干流發(fā)生大面積污染的可能性；根據(jù)物質(zhì)平衡原理和題目假設(shè)可知：水庫1中污染物的改變量 = 流入的污染物的量 流出污染物的量于是對于充分小的，在時間（t，t+）內(nèi)有：兩邊同除，并使0得： （1）現(xiàn)假設(shè)f(t)=p(t)v(t)得：即原式可以寫為： （2）在水庫1中發(fā)生撞船事故后，污染物處于非穩(wěn)定排放即：，而由于水庫閘門的關(guān)閉也勢必會引起水庫中水的體積變化，故：?，F(xiàn)不考慮流入水庫中的水所含有與泄漏污染物相同物質(zhì)的情況而帶來的影響，即可看作，另外由問題分析中知道：流出的污染物的濃度應與水庫中污染物濃度相同，即這樣對于問題一我們可以得到求解公式： （3）進一步我們假設(shè)從水庫中流出的水的流量初始值

11、（從t=0時算起）為，在關(guān)閉閘門的過程中，我們假定流量處于線性變化的趨勢。這一假設(shè)是基于流量與過流面積為線性關(guān)系上作出的，進一步可得： （4）從上面可看出為分數(shù)函數(shù)，這主要是因為水庫閘門關(guān)閉是在事故發(fā)生一小時后作出的?，F(xiàn)在有了的變化的表達式，為了能求出的表達式。我們還要寫出的表達式。首先我們假設(shè)水庫的體積的初始值為（t=0時），值我們可以通過衛(wèi)星定位系統(tǒng)及所建立的模型求出（具體衛(wèi)星定位系統(tǒng)模型見附表）。而跟相關(guān)的還有的值。我們假設(shè)為一定值，則隨隨時間變化的關(guān)系式為： 由于為分段函數(shù)可知：也響應的為分段函數(shù)，具體函數(shù)表達式為： （5）把（4）式代入（3）式可以得到：當秒時： （6）當秒時： （7

12、）對（5）式化簡有： （8）通過推導的出： （9）有已知條件可知：，故經(jīng)簡化后： 時 （10）對（7）簡化后得；（11）設(shè)，(r1-2r0)=b，得：當 （12）設(shè)：最后得到： （13）保持連續(xù)性，當t=3600時，=，此時可得到相應的值，但由于不確定因素很多，故確定不是很容易，這主要是缺少數(shù)據(jù)造成的。當時得到（14）式 （14）同樣為保持連續(xù)性，要求當t=3600時，=。最后： （15）那么時間內(nèi)流出水庫的污染物的量便可表示為： （16）在（s）時流出的污染物的量為： （17）在（s）流出的量為： （18）流出的總量： （19）再用Q與2/3q進行比較，便得出是否會發(fā)生大面積污染。問題二：如

13、果在另外的一水庫中有一化工廠均在晚上9：00-12：00違規(guī)進行周期性排放同樣含有該化學物質(zhì)的廢料。討論由于此次事故的發(fā)生，干流發(fā)生大面積污染的可能性。由題目分析中知道，水庫2中有一化工廠違規(guī)排放含有該污染物的廢料，而由于兩個水庫之間沒有聯(lián)系，故我們只需要單獨考慮水庫2的排污量，然后加上水庫1的污染物排放量，最后綜合考慮兩個水庫所排污染物的總量對干流的影響就可以了。下面我們將在水庫1模型的基礎(chǔ)上建立2水庫的模型。 （20）設(shè)（常數(shù)），即在9：0012：00這段時間內(nèi)污染物以一個恒定值流入水庫，考慮水庫水的流入速度為一定值，流出速度也為一個定值這樣水庫體積的表達式可寫成如下的公式： （21）代入

14、上面的表達式可簡化為： （22） 用上式可求出，其中可通過求得。那么從9：0012：00這三個小時內(nèi)流出閘門的污染物的總量就可以求出來，污染物的量為： （23） （24）（為三小時內(nèi)從1水庫和2水庫流出的污染物總量）若則發(fā)生大面積污染；若則不會發(fā)生大面積污染；問題三：如果兩個水庫間有一條人工修建的水渠相連接，水渠中的水流流向不定，但保證兩水庫之間的水流能夠相互影響。那么問題二結(jié)果是否會改變？由題目分析可知，當兩水庫之間有一人工修建的水渠相互連通時，水渠中的水流必然是從高水位流向低水位，為了使模型簡化，我們有以下說明：1、由于兩水庫是連通的，因此在水庫1關(guān)閘前，兩水庫的液面必然是趨近于等高的，否

15、則必然有水從一個水庫流向另一個水庫。2、在事故發(fā)生后一個小時內(nèi)，考慮兩水庫的規(guī)模相當，且水庫1沒有關(guān)閉泄洪閘，此時兩水庫彼此不受影響。3、在事故發(fā)生一個小時后，考慮到水庫1要關(guān)閉泄洪閘，這勢必引起水庫1的水位上升，由說明1可知水庫1中的水必將通過水渠流向水庫2。從上面的分析可知，在有連通水渠的情況下，我們只需考慮從1水庫向2水庫的流入情況，而不必考慮從2水庫向1水庫的流入情況。下面我們就該問題給出進一步分析。在第一個小時內(nèi)：（即9：0010：00）此時，水庫1中還未發(fā)生事故，只有水庫2中在排放廢料，用問題2的模型我們可求出流出的污染物的量為： （25）（由（20）式確定）在第二個小時內(nèi)（即10

16、：0011：00）水庫1已發(fā)生裝船事故，而水庫2繼續(xù)排放廢料，但泄洪閘還未關(guān)閉，所以我們?nèi)匀华毩⒖紤]。對水庫1我們用問題一的模型求解得： （26）對水庫2我們?nèi)杂袉栴}二的模型求解得到： （由（20）確定） （27） （28）在第三個小時內(nèi)（即11：0012：00）；水庫1泄洪閘正在關(guān)閉，而水庫2繼續(xù)排放廢料。這時水庫1和水庫2就要結(jié)合在一起考慮了，如下圖所示：圖 一由上面分析可列出以下的方程式：對于水庫1： （29）其中的表達式為： （）而： （）對水庫2： （30）其中 由基本假設(shè)知道，水庫1和水庫2規(guī)模相當則：假定為一定值。聯(lián)立解出：由于所得表達式非常復雜，我們把表達式放在附錄（三）里表示

17、。而由于過于復雜，我們這里就不給出解析解了。 （） （31）故有： （32）在12：00以后一個小時內(nèi)，1水庫完全關(guān)閉但1水庫的水將通過渠道流入2水庫內(nèi)。則有以下式：對水庫1有： （33）其中表達式為： 對水庫2有： （34）其中表達式為： 求得 （35） （36）若則干流不會發(fā)生大面積污染；若則干流會發(fā)生大面積污染；問題四：如果發(fā)生了大面積污染，那么針對第三種情況，試給出在短時間內(nèi)控制污染模型。5.2 模型二由問題的分析中知道，我們只能通過稀釋原理，建立污染物在干流中遷移遷移模型，使干流水體計算體積元內(nèi)該污染物的增量為負值時，從而使得在干流水體污染物的濃度低于危險警戒值時的濃度，才能在短時間

18、內(nèi)達到控制污染。為使模型清楚，我們先給出下面所用名詞解釋。源和漏：是對體積元內(nèi)污染物變化的一種描述，源是體積元內(nèi)污染物的增加速率，漏是體積元內(nèi)污染物的減少速率。這里的“漏”不意味著漏掉，而有更廣泛的意義，如污染物的降解、沉淀和揮發(fā)等都屬于“漏”。 由于污染物在干流中遷移符合遷移方程，由基本假設(shè)中我們知道河流中的水流處于推流狀態(tài)，也就是體積元中水分子以同一速度向下游運動如圖二和圖三。設(shè)C(x)和C(x+)分別為進入水片的水中某中污染物的濃度，是計算體積元中所含的該污染物的質(zhì)量。按照質(zhì)量守衡原理，計算體積元里污染物質(zhì)量的增量為： （1） （2） （3） （4） + （37） （5）上式中（1）至（

19、5）各項的意義如下（1）儲存量項； （38）（2）平流輸送項；（3）側(cè)向的源和漏，其中是單位時間內(nèi)單位長度上的源和漏，它通常是由側(cè)向分布流量帶入的，因此可把寫為：這里是該分布流量中污染物的濃度；（4）表面的源和漏，是單位時間內(nèi)單位面積上的源和漏（）；（5）體積元內(nèi)的源和漏，是單位時間內(nèi)單位體積內(nèi)的源和漏（）；圖 二圖 三把（38）式帶入（37）式，用除方程兩邊并令和，則得： （39）式（10）就是推流時的污染物遷移方程，以后將廣泛地利用這個方程，在求解問題四時，我們只要、和的值為負數(shù)就可以了，即可以通過關(guān)閉人工水渠和部分關(guān)閉水庫2泄洪閘的手段實現(xiàn)短時間內(nèi)控制污染。六、模型的擴展 從問題的分析中

20、可以知道，我們知道污染物在水中達到分布均勻是有一個遷移過程，符合污染物一維遷移方程，為使模型比較符合實際，我們對上面模型進行擴展，建立污染物一維遷移擴展模型。擴展模型一從水庫中流出的污染物的量可用下面的式子表示：由于在水庫中污染物濃度在各個位置是不同的，而在整個水庫中各個量又具有一定的連續(xù)性，故考慮時刻在閘門口處時間內(nèi)通過的量為：，此式中我們考慮與，的關(guān)系，其中方向為水庫流水的流向。這是因為對于事故泄漏排放的污染物，人們最關(guān)心的是其通過某一位置的時間，最大濃度等數(shù)值。因此，研究一維流場中的瞬時排放污染物的分布特征就非常必要了，我們用下式為基礎(chǔ)來考慮分析泄漏時污染物在閘門口的濃度變化特征。在瞬時

21、排放條件下：，此時有下式： （40）考慮在泄洪閘處：為泄洪閘處到污染處的距離，且（守恒污染物） （41）其中： 考慮在一個極其短暫的時間內(nèi)，在這個內(nèi)污染物的濃度可近似看作是不變的。則在時刻有： 通過對上式積分可求出2個小時內(nèi)流出水庫的污染物的量為： （42）這里P，A，Q，均為已知，只是發(fā)生的地點未知。通過上式我們可以求出在兩個小時內(nèi)通過的量為噸時的，但由于式（42）的解析解很難求出，我們就把它寫成下面的形式： （43）把（2）式的相應的量代入后，通過遺傳算法可求出的最小值，而（43）式的最小值為零，由此可求出滿足時所對應的，求出后，就可以得到發(fā)生大面積污染的可能的點所在的范圍，如下圖所示：A

22、2A1 圖 四A1為發(fā)生大面積污染的點的集合。A2為不發(fā)生大面積污染的點的集合。因此對于問題一：發(fā)生大面積污染的可能為： ； （44）問題二：由題目分析知道，若水庫2中有一工廠在9：0012：00違規(guī)進行排放，可知這樣會增大發(fā)生大面積污染的可能性。因為另外一水庫處于穩(wěn)定流動狀態(tài)，工廠排放的污染物為穩(wěn)定的連續(xù)排放，環(huán)境中污染物的分布狀況也是穩(wěn)定的，綜合考慮可有下式： （45）其中： ，。我們假定工廠距泄洪閘的距離可以確定為則： 污染物到達泄洪閘的時間為：則泄漏的量為： （46）由此我們可以得到在另外一水庫當中有工廠進行周期性排放下流出的總量M為： （47）即我們把上式轉(zhuǎn)化為的形式，這里我們可以看

23、出的表達式與問題（1）中有所不同，這是因為我們考慮的時間已從原來的2小時變成3小時了，對于問題（二）我們?nèi)匀挥脝栴}（一）的解法，即通過遺傳算法(程序見附錄四)求出發(fā)生大面積污染時的。通過這個我們可求出發(fā)生大面積污染的點的集合，進而求出發(fā)生大面積污染的可能性為，A是水庫的面積,可知是個定值。 在這里我們可以給出擴展模型一求解下的事例； 其中：m/s =1513 部分數(shù)據(jù)見附錄二。 通過MATLAB求解得到結(jié)果如下： 問題一發(fā)生大面積污染的點的集合類圓半徑； 問題二發(fā)生大面積污染的點的集合類圓半徑。 圖形見圖五圖六：圖五 問題一發(fā)生大面積污染的可能半徑圖六 問題二發(fā)生大面積污染的可能半徑 擴展二從

24、問題的分析及假設(shè)可以知道，我們的模型在建立過程中我們所考慮的只是一種污染物，即此污染物為易溶、不揮發(fā)的物質(zhì)。實際上污染物可能為油性、揮發(fā)性和沉淀物質(zhì)。對于油性物質(zhì)我們應該考慮該物質(zhì)的密度是否大于水的密度，若油性污染物小于水的密度，因為油性物質(zhì)不溶于水，該污染物應該是漂浮在水面上的，而且呈油膜狀，因此我們就應考慮水庫表層水的流動情況；若油性污染物的密度大于水的密度，則污染物將會沉如水底，那么我們就應該考慮深沉水體的水流問題；若油性污染物等于水的密度，那么污染物就會懸浮在水庫水體中，因此我們就應該考慮整個水體上、中、下層的流動情況（因為不同水層的水流的流動情況是不同的）。對于揮發(fā)性的污染物，我們還

25、需要考慮此污染物在水和空氣中的分配比例。在考慮問題四的時候我們考慮的是推流問題，實際污染隨水流遷移是以水團形式的，這樣我們就需要考慮水體不同水層的流動問題而。為使模型更符合實際，模型必須在上述方面在做深入探討。七、模型的評價模型的優(yōu)點： 本模型在建立模型前，由于水體運動非常復雜，污染物的性質(zhì)也多種多樣，污染物在水中的遷移情況受各種因素限制，不容易全部考慮。本文綜合考慮上述因素，先通過合理的假設(shè)，首先從易溶急難揮發(fā)性污染物入手，考慮污染物和庫水混合均勻的情況，利用物質(zhì)平衡原理，針對問題中各種情況建立物質(zhì)平衡方程模型。在此基礎(chǔ)上考慮污染物的一維遷移模型。 在建立一維遷移擴展模型時，我們通過轉(zhuǎn)化思想

26、，將污染物對干流形成大面積污染的可能轉(zhuǎn)換為求解在關(guān)閉泄洪閘的過程中，水庫中多大面積水域內(nèi)發(fā)生事故時才會對干流形成大面積污染，使得問題求解方案明了，為了驗證模型的正確性，我們通過遺傳算法對事例進行求解，很好的驗證了模型。模型的缺點：模型在建立過程中我們所考慮的此污染物為易溶、不揮發(fā)的物質(zhì)，而沒有考慮污染物為油性、揮發(fā)性和沉淀物質(zhì)。另外，我們在考慮水體流動時沒有考慮水體上、中、下層的不同流動情況。在問題四的求解時，我們只考慮推流問題，而實際污染隨水流遷移是以水團形式的，這樣我們就需要考慮水體不同水層的流動問題。時間有限，未能對其他方面做深入的探討，使得模型十分完美。參考書目1 W.金士博聯(lián)邦德國

27、水環(huán)境數(shù)學模型 中國建筑工業(yè)出版社 1987.102 岳天祥 資源環(huán)境數(shù)學模型手冊 科學出版社 2003.103 余常昭 環(huán)境流體力學導論 清華大學出版社 1992.104 劉振航 數(shù)學建模 中國人民大學出版社 2004.055 陳春云 鄭彤 環(huán)境系統(tǒng)數(shù)學模型 化學工業(yè)出版社 2003.046 周建華 黃燕 MATLAB5.3 北京大學出版社 2000.127 丁春利 精通MATLAB6 清華大學出版社 2002.06附錄一：水庫實際庫容及水域面積計算公式：式中，s為DEM象元面積；H為水庫水位高程；為庫域DEM各像元高程值；n為庫域DEN像元數(shù)。（曾永年，馬海州，沙占江等；龍羊峽庫區(qū)環(huán)境動態(tài)

28、監(jiān)測信息系統(tǒng)的建立與應用。遙感學報，2000.4）附錄二：序號項目單位數(shù)據(jù)相應時間1入庫水量億立方米469.0722入庫平均流量立方米/秒14363出庫水量億立方米494.287 4出庫平均流量立方米/秒15135最大入庫流量立方米/秒281002003-9-8 8:00 6最小入庫流量立方米/秒02004-2-20 20:00 7最大出庫流量立方米/秒120002003-9-3 2:00 8最小出庫流量立方米/秒4272004-7-10 8:00 9最高庫水位米156.982003-10-17 14:00 10最低庫水位米138.32004-7-16 8:00 11大廠發(fā)電量億千瓦時44.6

29、34512平均出力萬千瓦49.2013最高負荷萬千瓦109.62003-10-26 14最低負荷萬千瓦11.22004-7-25 15最大調(diào)峰萬千瓦69.82003-10-216最小調(diào)峰萬千瓦42.42004-9-12 17平均調(diào)峰萬千瓦32.6818調(diào)峰電量億千瓦時.0912 19棄水量億立方米119.94620平均耗水率立方米/千瓦時7.79212x2發(fā)電量億千瓦時3.425422大、小廠總電量億千瓦時48.0599附錄三：附錄四：遺傳算法程序function Genetic(AimFunc)% This is simple genetic algorithm(SGA)% In this

30、function ,it fulfils genetic algorithm%-These can be modified as you like-maxgen=200; % maximum generationsizepop=100; % size of population % AimFunc=strimFunc; % this is function of counting fitnessfselect=tournament; % method of select % you can choose tournament;roulettefcode=float; % method of c

31、oding % you can choose float;grey;binary pcross=0.6; % probablity of crossover,between 0 and 1fcross=float; % method of crossover % you can choose float;simple;uniformpmutation=0.2; % probability of mutation,between 0 and 1 fmutation=float; % method of mutation % you can choose float;simple;lenchrom

32、=1 ; % length of bit of every variblebound=0 100000; %-individuals=struct(fitness,zeros(1,sizepop),.%value,zeros(1,sizepop),. chrom,); % structure of populationavgfitness=; % average fitness of population bestfitness=; % best fitness of populationbestchrom=; % chromosome of best fitness% inivitializ

33、ationfor i=1:sizepop % produce new population at random individuals.chrom(i,:)=Code(lenchrom,fcode,bound); x=Decode(lenchrom,bound,individuals.chrom(i,:),fcode); individuals.fitness(i)=Aimfunc(x);end% find minimum value which is bestbestfitness bestindex=min(individuals.fitness);bestchrom=individual

34、s.chrom(bestindex,:);avgfitness=sum(individuals.fitness)/sizepop;% record average and best fitness of every generationtrace=avgfitness bestfitness; % evolution beginfor i=1:maxgen % selection individuals=Select(individuals,sizepop,fselect); avgfitness=sum(individuals.fitness)/sizepop; % crossover in

35、dividuals.chrom=Cross(pcross,lenchrom,individuals.chrom,. sizepop,fcross,i maxgen); % mutation individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,. sizepop,fmutation,i maxgen,bound); % calculate fitness for j=1:sizepop x=Decode(lenchrom,bound,individuals.chrom(j,:),fcode); individuals.f

36、itness(j)=Aimfunc(x); end % substitute chromosome of worest fitness % find minimum value which is best newbestfitness,newbestindex=min(individuals.fitness); worestfitness,worestindex=max(individuals.fitness); % substitute chromosome of worest fitness if bestfitnessnewbestfitness bestfitness=newbestf

37、itness; bestchrom=individuals.chrom(newbestindex,:); end individuals.chrom(worestindex,:)=bestchrom; individuals.fitness(worestindex)=bestfitness; avgfitness=sum(individuals.fitness)/sizepop; trace=trace;avgfitness bestfitness; end % draw fitness of every generationhfig=findobj(Tag,trace);% See if i

38、t is openif ishandle(hfig) figure(hfig);else hfig=figure(Tag,trace);endfigure(hfig);r c=size(trace);plot(1:r,trace(:,1),r-,1:r,trace(:,2),b-);title(適應度曲線 終止代數(shù) num2str(maxgen);xlabel(進化代數(shù));ylabel(適應度);legend(平均適應度,最佳適應度);disp(適應度 變量);x=Decode(lenchrom,bound,bestchrom,fcode);% show in command windowvp

39、a(bestfitness,10);vpa(x,10);disp(bestfitness x);function ret=AimFunc(x)% 求最小值%ret=sum(x);% 求最小值z=quadl(f,0.0001,7200,x);%求最大值%ret=1/sum(x);function y=f(t,x)y=7565./sqrt(4.8*pi*t).*exp(-(x-3*t).(2)./(4.8*t);function ret=Code(lenchrom,opts,bound)% In this function ,it converts a set varibles into a ch

40、romosome% lenchrom input : length of chromosome% opts input : tag of coding method% bound input : boundary of varibles% ret output: chromosomeswitch opts case binary % binary coding pick=rand(1,sum(lenchrom); bits=ceil(pick-0.5); temp=sum(lenchrom)-1:-1:0; ret=sum(bits.*(2.temp); case grey % grey co

41、ding pick=rand(1,sum(lenchrom); bits=ceil(pick-0.5); greybits=bits; for i=2:length(greybits) greybits(i)=bitxor(bits(i-1),bits(i); end temp=sum(lenchrom)-1:-1:0; ret=sum(greybits.*(2.temp); case float % float coding pick=rand(1,length(lenchrom); ret=bound(:,1)+(bound(:,2)-bound(:,1).*pick;endfunctio

42、n ret=Cross(pcross,lenchrom,chrom,sizepop,opts,pop)% In this function,it fulfils a crossover among chromosomes% pcorss input : probability of crossover% lenchrom input : length of a chromosome% chrom input : set of all chromosomes% sizepop input : size of population% opts input : tag for choosing me

43、thod of crossover% pop input : current serial number of generation and maximum gemeration% ret output : new set of chromosomeswitch opts case simple % cross at single position for i=1:sizepop % select two children at random pick=rand(1,2); index=ceil(pick.*sizepop); while prod(pick)=0 | index(1)=ind

44、ex(2) pick=rand(1,2); index=ceil(pick.*sizepop); end % probability of crossover pick=rand; if pickpcross continue; end % random position of crossover pick=rand; while pick=0 pick=rand; end pos=ceil(pick.*sum(lenchrom); tail1=bitand(chrom(index(1),2.pos-1); tail2=bitand(chrom(index(2),2.pos-1); chrom

45、(index(1)=chrom(index(1)-tail1+tail2; chrom(index(2)=chrom(index(2)-tail2+tail1; end ret=chrom; case uniform % uniform cross for i=1:sizepop % select two children at random pick=rand(1,2); while prod(pick)=0 pick=rand(1,2); end index=ceil(pick.*sizepop); % random position of crossover pick=rand; whi

46、le pick=0 pick=rand; end if pickpcross continue; end % random position of crossover pick=rand; while pick=0 pick=rand; end mask=2ceil(pick*sum(lenchrom); chrom1=chrom(index(1); chrom2=chrom(index(2); for j=1:sum(lenchrom) v=bitget(mask,j); % from lower to higher bit if v=1 chrom1=bitset(chrom1,. j,b

47、itget(chrom(index(2),j); chrom2=bitset(chrom2,. j,bitget(chrom(index(1),j); end end chrom(index(1)=chrom1; chrom(index(2)=chrom2; end ret=chrom; case float for i=1:sizepop % select two children at random pick=rand(1,2); while prod(pick)=0 pick=rand(1,2); end index=ceil(pick.*sizepop); % random posit

48、ion of crossover pick=rand; while pick=0 pick=rand; end if pickpcross continue; end % random position of crossover pick=rand; while pick=0 pick=rand; end pos=ceil(pick.*sum(lenchrom); pick=rand; v1=chrom(index(1),pos); v2=chrom(index(2),pos); chrom(index(1),pos)=pick*v2+(1-pick)*v1; chrom(index(2),p

49、os)=pick*v1+(1-pick)*v2; end ret=chrom;endfunction ret=Decode(lenchrom,bound,code,opts)% In this function ,it decode chromosome % lenchrom input : length of chromosome% opts input : tag of coding method% bound input : boundary of varibles% ret output: value of variblesswitch opts case binary % binar

50、y coding for i=length(lenchrom):-1:1 data(i)=bitand(code,2lenchrom(i)-1); code=(code-data(i)/(2lenchrom(i); end ret=bound(:,1)+data./(2.lenchrom-1).*(bound(:,2)-bound(:,1); case grey % grey coding for i=sum(lenchrom):-1:2 code=bitset(code,i-1,bitxor(bitget(code,i),bitget(code,i-1); end for i=length(

51、lenchrom):-1:1 data(i)=bitand(code,2lenchrom(i)-1); code=(code-data(i)/(2lenchrom(i); end ret=bound(:,1)+data./(2.lenchrom-1).*(bound(:,2)-bound(:,1); case float % float coding ret=code;endfunction y=f1(t)function ret=Mutation(pmutation,lenchrom,chrom,sizepop,opts,pop,bound)% In this function,it ful

52、fils a mutation among chromosomes% pcorss input : probability of mutation% lenchrom input : length of a chromosome% chrom input : set of all chromosomes% sizepop input : size of population% opts input : tag for choosing method of crossover% pop input : current serial number of generation and maximum gemeration% ret output : new set of chromosomeswitch opts case simple % mutation at single position for i=1:sizepop % select child at random pick=rand; while pick=0 pick=rand; end index=ceil(pick*sizepop); pick=rand; if pickpmutation