1、安徽工程大學安徽工程大學 數學建模（選修課）課程論文數學建模（選修課）課程論文 題目：水庫排污問題題目：水庫排污問題 摘要：摘要： 本文主要需解決的問題是：本文主要需解決的問題是：討論由于此次事故的發生，干流發生大面積污染的可能性，以及如何避免該事故的發生。 針對問題一，我們只要考慮在事故發生到關閉水庫的兩個小時內，流出水庫的污染物的質量小于 q32噸。 針對問題二，我們只需單獨計算1,2號水庫的污染情況，然后再相加就可以了。 針對問題三，建立人工水渠就是在問題二的基礎上使水庫1和水庫2產生聯系，我們只需考慮從水庫1到水庫2的情況。 針對問題四，干流已經出現了大面積污染，在一些外在條件短時間內

2、無法干擾的情況下，我們可以采取稀釋的辦法，是單位體積內的化學物質的濃度少于危險警戒值時的濃度，這樣就可以短時間內控制污染。關鍵詞：控制 人工水渠 稀釋 姓名：殷競存 專業：數學與應用數學 班級：數學 112 學號：3110801234 指導老師：周金明指導老師：周金明 成績：成績： 一：問題重述一：問題重述1、 問題的背景某條江流上有 2 條支流，每條支流上都興建了規模相當的水庫。由于正處在雨水多發季節，因此兩個水庫都以一定規模的流量進行泄洪。某天晚上10：00，在其中的一個水庫中發生了兩船相撞的事故，而其中的一條船裝載的 p 噸化學物質（這里的化學物質可以是具有揮發性的，也可能是急難揮發的）

3、全部泄漏至水庫中。當水上航運事故處置中心接獲事故報告，立即要求該水庫關閉水庫泄洪閘，以免化學物質隨洪水流入干流，發生更大規模的污染。水庫閘門開始關閉時，已經處在事故發生后的 1 個小時，而水庫閘門徹底關閉也需要 1 個小時的時間。 根據當地環境監測的有關規定，干流大面積污染的危險警戒值設為：三小時內 q 噸該化學物質發生泄漏。2、面臨的問題(1) 試建立合理的數學模型，討論由于此次事故的發生，干流發生大面積污染的可能性；(2) 如果在另外的一水庫中有一化工廠違規排放廢料。廢料中同樣含有該化學物質。該工廠為躲避環境監測站的監控，均在晚上 9：00-12：00 違規進行周期性排放。在這種情形下，討

4、論由于此次事故的發生，干流發生大面積污染的可能性；(3) 如果以上兩個水庫間有一條人工修建的水渠相連接，水渠中的水流流向不定，但保證兩水庫之間的水流能夠相互影響。那么上述結果是否會改變？請給出說明，若有改變，則給出修正的模型及結果；(4) 如果發生了大面積污染，那么針對第三種情況，試給出在短時間內控制污染模型。 二：問題的假設（1）污染物為速溶物質，因此藥品從船上流入水中的時間很少，可以忽略不計；（2）污染物質從水庫中一經流出就進入干流；（3）水庫和河流中的水流都是處于推流狀態；（4）兩水庫事故發生條件相同，即兩水庫有相同的客觀條件；（5）被污染的水庫關閉泄洪閘后不再有水流流入干流；（6）不考

5、慮生物等因素在水庫泄洪過程中的作用，污染物除了流出外不因腐爛沉積等手段從水中消失；（7）外界因素不對水庫的體積變化產生影響，例如：雨水、地表徑流、底下徑流等；（8）參與模型的變量是連續變化的，并且充分光滑；（9）不考慮從不同的渠道流入與流出水庫之間的區別，只考慮攜帶污染物的水流入水庫和水庫中的水流出對水庫污染程度的影響，因此可以把水庫看成是單流入單流出的系統。 三：符號的約定三：符號的約定（1）：t 時刻水庫水的流入速度；)(1tr（2）：t 時刻流入水庫的污染物的濃度；)(1tp（3）：t 時刻水庫水的流出的速度；)(0tr（4）：t 時刻流出水庫的污染物的濃度；)(0tp（5）：t 時刻水

6、庫中污染物的濃度；)(tp（6）：t 時刻水庫水的體積；)(tv（7）：計算體積元內該污染物的增量；pm（8） ： 為時間；t（9） ： 為從水庫中流出的水中的污染物的濃度；C（10） ：為水庫的流出速度，即流量；Q（11）：為泄洪閘處到污染處的距離；X（12） Q：為水庫的流量；（13） q：為排入河流的污水的流量；（14） ：為河流中污染物的本底濃度；1C（15）:為水庫中的污染物的濃度；2C 四：模型的分析與建立四：模型的分析與建立由問題的分析中知道，流入水庫的污染物能以很快的速度與水庫中的水均勻混合，也就是說水庫中的污染狀況在任何局部水體都是一樣的，污染程度與水體在水庫中的位置無關，因

7、此我們可以建立下面模型。模型一問題一：此次事故的發生，干流發生大面積污染的可能性；根據物質平衡原理和題目假設可知：水庫 1 中污染物的改變量 = 流入的污染物的量 流出污染物的量于是對于充分小的，在時間（t，t+）內有：ttttrtptrtptvtpttvttp)()()()()()()()(0011兩邊同除，并使0 得：tt (1) )()()()()()()()(0011trtptrtpttvtpttvttp 現假設 f(t)=p(t)v(t)得：dtdtvdtdtpdtddtdttfttftptvtvtptf)()()()()()()()()(即原式可以寫為： (2) )()()()()

8、()(0011)()(trtptrtpdtdtvdtdtptptv在水庫 1 中發生撞船事故后，污染物處于非穩定排放即：，而由于水庫閘門的關閉也勢必會引0)(dtdtp起水庫中水的體積變化，故：。現不考慮流入0)(dtdtv水庫中的水所含有與泄漏污染物相同物質的情況而帶來的影響，即可看作，另外由問題分析中知道：0)(1tp流出的污染物的濃度應與水庫中污染物濃度相同，即這樣對于問題一我們可以得到求解公式：)()(0tptp (3) )()()()(0)()(trtpdtdtvdtdtptptv進一步我們假設從水庫中流出的水的流量初始值（從 t=0 時算起）為 ，在關閉閘門的過程中，我們0r假定流

9、量處于線性變化的趨勢。這一假設是基于流量與過流面積為線性關系上作出的，進一步可得： （4) )72003600.(.36002)36000.(.)(0000ttrrtrtr從上面可看出為分數函數，這主要是因為水庫)(0tr閘門關閉是在事故發生一小時后作出的。現在有了的變化的表達式，為了能求出的表達式。我們)(0tr)(tp還要寫出的表達式。首先我們假設水庫的體積的)(tv初始值為 （t=0 時） ， 值我們可以通過衛星定位系0v0v統及所建立的模型求出（具體衛星定位系統模型見附表） 。而跟相關的還有的值。我們假設為一)(tv)(1tr)(1tr定值 ，則隨隨時間變化的關系式為：1r)(tv t

10、trtrvtv)()(010由于為分段函數可知：也響應的為分段函數，)(0tr)(tv具體函數表達式為： (5) )72003600.(.)36002()36000.(.)(0010010tttrrtrvttrtrvtv把（4）式代入（3）式可以得到：當秒時：36000 t (6) 0)(01001)()()()(rtpdtdtrtrvrrtptp當秒時：72003600 t(7) )36002()()36002()18002()(00)(0010001trrtpdtdttrrtrvtrrrtptp對（5）式化簡有： (8) 1)(010)()(rtpdtdtrtrvtp通過推導的出： (9)

11、 cetptrtrvrrr)ln(1010010)(有已知條件可知：，故0/)0(vPp10100/rrrvvpc經簡化后： 00010/)()(011011vpvtrrvtprrrrrr(10)36000 t 對（7）簡化后得；(2001002001013600)2()3600(3600)2()()(lntvtrrvdttrtrtrrvdtrtp11)設，(r1-2r0)=b，得：cv 0ar36000當acb42)ln(21)3600()42arctan(42)()(ln2022ctbtaarbacbatbacrtp (12) )42arctan(42222bacbatbacab設：)42

12、arctan(4222bacbatbacz最后得到： (13) 1)72005.0()ln(5.012)(cetpzabcbtat保持連續性，當 t=3600 時，=，此時可)3600(1p)3600(p得到相應的 值，但由于不確定因素很多，故確定1C不是很容易，這主要是缺少數據造成的。1C當時得到（14）式acb42 1)72002()ln(2112)(Cetpzabcbtat （14）cacbbatacbbatacbz4242ln41222同樣為保持連續性，要求當 t=3600 時，=)3600(1p。)3600(p最后：)(72003600.4.)(72003600.4.)(36000.

13、)()(21)72002()ln(2121)72005 . 0()ln(5 . 0001022011011stacbcestacbcestvtrrvtpzabcbtatzabcbtatrrrrrr當當 （15）那么時間內流出水庫的污染物的量便可表示為：t ttrtpQ)()(0（16）在（s）時流出的污染物的量為：36000 t t drvpvtrrvQrrrrrr360000000101/)(011011（17）在（s）流出的量為：72003600 t 72003600002)36002()(dttrrtpQ（18）流出的總量： 21QQQ（19）再用 Q 與 2/3q 進行比較，便得出是否

14、會發生大面積污染。問題二：如果在另外的一水庫中有一化工廠均在晚上9：00-12：00 違規進行周期性排放同樣含有該化學物質的廢料。討論由于此次事故的發生，干流發生大面積污染的可能性。由題目分析中知道，水庫 2 中有一化工廠違規排放含有該污染物的廢料，而由于兩個水庫之間沒有聯系，故我們只需要單獨考慮水庫 2 的排污量，然后加上水庫 1 的污染物排放量，最后綜合考慮兩個水庫所排污染物的總量對干流的影響就可以了。下面我們將在水庫 1 模型的基礎上建立 2 水庫的模型。 )()()()()()(0011)()(trtptrtpdtdtvdtdtptptv（20）設（常數） ，即在 9：0012：00

15、這段時間內ktp)(1污染物以一個恒定值流入水庫，考慮水庫水的流入速度為一定值 ，流出速度也為一個定值1r0r這樣水庫體積的表達式可寫成如下的公式： trtrvtv010)()108000( t（21）代入上面的表達式可簡化為： 01)(010)()(rtprkdtdtrrvtp（22） 用上式可求出，其中可通2010211)()(Ctrrvktprrr2C過求得。那么從 9：0012：00 這三個小時內流kp)0(出閘門的污染物的總量就可以求出來，污染物的量為： 1080000)(dtrtpQ（23） QQQ（24）（ 為三小時內從 1 水庫和 2 水庫流出的污染物總Q量）若則發生大面積污染

16、；qQ 若則不會發生大面積污染；qQ 問題三：如果兩個水庫間有一條人工修建的水渠相連接，水渠中的水流流向不定，但保證兩水庫之間的水流能夠相互影響。那么問題二結果是否會改變？由題目分析可知，當兩水庫之間有一人工修建的水渠相互連通時，水渠中的水流必然是從高水位流向低水位，為了使模型簡化，我們有以下說明：1、由于兩水庫是連通的，因此在水庫 1 關閘前，兩水庫的液面必然是趨近于等高的，否則必然有水從一個水庫流向另一個水庫。2、在事故發生后一個小時內，考慮兩水庫的規模相當，且水庫 1 沒有關閉泄洪閘，此時兩水庫彼此不受影響。3、在事故發生一個小時后，考慮到水庫 1 要關閉泄洪閘，這勢必引起水庫 1 的水

17、位上升，由說明 1 可知水庫 1 中的水必將通過水渠流向水庫 2。從上面的分析可知，在有連通水渠的情況下，我們只需考慮從 1 水庫向 2 水庫的流入情況，而不必考慮從 2 水庫向 1 水庫的流入情況。下面我們就該問題給出進一步分析。在第一個小時內：（即 9：0010：00）此時，水庫 1 中還未發生事故，只有水庫 2 中在排放廢料，用問題 2 的模型我們可求出流出的污染物的量為： 3600001)(dtrtpQ（25）（由（20）式確定）)(tp在第二個小時內（即 10：0011：00）水庫 1 已發生裝船事故，而水庫 2 繼續排放廢料，但泄洪閘還未關閉，所以我們仍然獨立考慮。對水庫 1 我們

18、用問題一的模型求解得： 3600000001021/)(011011dtrvpvtrrvQrrrrrr（26）對水庫 2 我們仍有問題二的模型求解得到： （由（20）確定） 72003600022)(dtrtpQ)(tp（27 22212QQQ（28）在第三個小時內（即 11：0012：00） ；水庫 1 泄洪閘正在關閉，而水庫 2 繼續排放廢料。這時水庫 1 和水庫 2 就要結合在一起考慮了，如下圖所示：圖 一由上面分析可列出以下的方程式：對于水庫 1： （29）)()()()()()( 00)()(trtptrtpdtdtvdtdtptptv其中的表達式為：)(tv （）ttrttrrtr

19、vtv)()36002()( 00010172003600 t而： 36002)(000trrtr（）72003600 t對水庫 2： )()()()()()()()(0 011)()(trtptrtptrtpdtdtvdtdtptptv（30）其中 trttrtrvtv0 0102)()(由基本假設知道，水庫 1 和水庫 2 規模相當則：trttrtrvtv0 010)()(假定為一定值 。)( 0tr 0r聯立解出：由于所得表達式非常復雜，我)()(tptp與們把表達式放在附錄（三）里表示。而由于過)(tp)(tp于復雜，我們這里就不給出解析解了。 （）72003600031)()(dtt

20、rtpQ72003600 t dtrtpQ108007200032)(（31）36002)(010trrtr故有： 32313QQQ（32）在 12：00 以后一個小時內，1 水庫完全關閉但1 水庫的水將通過渠道流入 2 水庫內。則有以下式：對水庫 1 有： )()()()( 0)()(trtpdtdtvdtdtptptv（33）其中表達式為： )(tvtrtrvtv 010)(對水庫 2 有： （34）)()()()()()(0 0)()(trtptrtpdtdtvdtdtptptv其中表達式為： )(tvtrtrtrvtv01 00)(求得 dtrtpQ1440010800042)(（35

21、） 42321QQQQQ總（36）若則干流不會發生大面積污染；qQ 若則干流會發生大面積污染；qQ 問題四：如果發生了大面積污染，那么針對第三種情況，試給出在短時間內控制污染模型。由問題的分析中知道，我們只能通過稀釋原理，建立污染物在干流中遷移遷移模型，使干流水體計算體積元內該污染物的增量為負值時，從而使得在pm干流水體污染物的濃度低于危險警戒值時的濃度，才能在短時間內達到控制污染。模型二為使模型清楚，我們先給出下面所用名詞解釋。源和漏：是對體積元內污染物變化的一種描述，源是體積元內污染物的增加速率，漏是體積元內污染物的減少速率。這里的“漏”不意味著漏掉，而有更廣泛的意義，如污染物的降解、沉淀

22、和揮發等都屬于“漏” 。 由于污染物在干流中遷移符合遷移方程，由基本假設中我們知道河流中的水流處于推流狀態，也就是體積元中水分子以同一速度向下游運動如圖二和圖三。設 C(x)和 C(x+)分別為進入水片的水中某中污染物x的濃度，是計算體積元中所含的該污染物的質量。pm按照質量守衡原理，計算體積元里污染物質量的增量為： txbStxStxxQxxCxQxCmslp)()()()(（1） （2） （3） （4） txASv（37） （5）上式中（1）至（5）各項的意義如下（1）儲存量項； )(CAxmp（38）（2）平流輸送項；（3）側向的源和漏，其中 是單位時間內單位長lS度上的源和漏，它通常是

23、由側向分布流量 帶入的，q因此可把 寫為：lSilCqS這里 是該分布流量中污染物的濃度；iC（4）表面的源和漏， 是單位時間內單位面積sS上的源和漏（） ；smg2/（5）體積元內的源和漏， 是單位時間內單位vS體積內的源和漏（） ；smg3/圖 二圖 三把（38）式帶入（37）式，用除方程兩邊并令tx和，則得：0 x0t vslSASbSxQCtAC)()(（39）式（10）就是推流時的污染物遷移方程，以后將廣泛地利用這個方程，在求解問題四時，我們只要 、lS和 的值為負數就可以了，即可以通過關閉人工水sSvS渠和部分關閉水庫 2 泄洪閘的手段實現短時間內控制污染。 五：模型的評價與改進5

24、.1：模型的優點 本模型在建立模型前，由于水體運動非常復雜，污染物的性質也多種多樣，污染物在水中的遷移情況受各種因素限制，不容易全部考慮。本文綜合考慮上述因素，先通過合理的假設，首先從易溶急難揮發性污染物入手，考慮污染物和庫水混合均勻的情況，利用物質平衡原理，針對問題中各種情況建立物質平衡方程模型。在此基礎上考慮污染物的一維遷移模型。 在建立一維遷移擴展模型時，我們通過轉化思想，將污染物對干流形成大面積污染的可能轉換為求解在關閉泄洪閘的過程中，水庫中多大面積水域內發生事故時才會對干流形成大面積污染，使得問題求解方案明了，為了驗證模型的正確性，我們通過遺傳算法對事例進行求解，很好的驗證了模型。5

25、.2：模型的缺點模型在建立過程中我們所考慮的此污染物為易溶、不揮發的物質，而沒有考慮污染物為油性、揮發性和沉淀物質。另外，我們在考慮水體流動時沒有考慮水體上、中、下層的不同流動情況。在問題四的求解時，我們只考慮推流問題，而實際污染隨水流遷移是以水團形式的，這樣我們就需要考慮水體不同水層的流動問題。時間有限，未能對其他方面做深入的探討，使得模型十分完美。5.3：模型的改進從問題的分析及假設可以知道，我們的模型在建立過程中我們所考慮的只是一種污染物，即此污染物為易溶、不揮發的物質。實際上污染物可能為油性、揮發性和沉淀物質。對于油性物質我們應該考慮該物質的密度是否大于水的密度，若油性污染物小于水的密

26、度，因為油性物質不溶于水，該污染物應該是漂浮在水面上的，而且呈油膜狀，因此我們就應考慮水庫表層水的流動情況；若油性污染物的密度大于水的密度，則污染物將會沉如水底，那么我們就應該考慮深沉水體的水流問題；若油性污染物等于水的密度，那么污染物就會懸浮在水庫水體中，因此我們就應該考慮整個水體上、中、下層的流動情況（因為不同水層的水流的流動情況是不同的） 。對于揮發性的污染物，我們還需要考慮此污染物在水和空氣中的分配比例。在考慮問題四的時候我們考慮的是推流問題，實際污染隨水流遷移是以水團形式的，這樣我們就需要考慮水體不同水層的流動問題而。為使模型更符合實際，模型必須在上述方面在做深入探討。參考文件：【1

27、】應玖茜，魏權齡。非線性規劃與理論。北京：中華人民大學出版社，1994【2】岳天祥。資源環境數學模型手冊。北京：科學出版社，2003【3】陳春云，鄭彤。環境系統數學模型。北京：化學工業出版社，2003【4】周建華，黃燕 。MATLAB5.3。北京：北京大學出版社，2000【5】趙靜，但琦。數學建模與數學實驗。北京：高等教育出版社，2007附錄一：水庫實際庫容及水域面積計算公式：)(ihHsVsnS 式中，s 為 DEM 象元面積；H 為水庫水位高程； 為ih庫域 DEM 各像元高程值；n 為庫域 DEN 像元數。（曾永年，馬海州，沙占江等；龍羊峽庫區環境動態監測信息系統的建立與應用。遙感學報，

28、2000.4）附錄二：序號項目單位數據相應時間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大廠發電量億千瓦時44.634512平均出力萬千瓦49.201

29、3最高負荷萬千瓦109.62003-10-26 14最低負荷萬千瓦11.22004-7-25 15最大調峰萬千瓦69.82003-10-216最小調峰萬千瓦42.42004-9-12 17平均調峰萬千瓦32.6818調峰電量億千瓦時.0912 19棄水量億立方119.946米20平均耗水率立方米/千瓦時7.79212x2 發電量億千瓦時3.425422大、小廠總電量億千瓦時48.0599附錄一：水庫實際庫容及水域面積計算公式：)(ihHsVsnS 式中，s 為 DEM 象元面積；H 為水庫水位高程； 為ih庫域 DEM 各像元高程值；n 為庫域 DEN 像元數。（曾永年，馬海州，沙占江等；龍羊

30、峽庫區環境動態監測信息系統的建立與應用。遙感學報，2000.4）附錄二：序號項目單位數據相應時間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大廠發電量億

31、千瓦時44.634512平均出力萬千瓦49.2013最高負荷萬千瓦109.62003-10-26 14最低負荷萬千瓦11.22004-7-25 15最大調峰萬千瓦69.82003-10-216最小調峰萬千瓦42.42004-9-12 17平均調峰萬千瓦32.6818調峰電量億千瓦時.0912 19棄水量億立方米119.94620平均耗水率立方米/千瓦時7.79212x2 發電量億千瓦時3.425422大、小廠總電量億千瓦時48.0599附錄三)1%306030301arctan(:2%90036003600180036009002:1%)1%/ )2%152%301%)180036001800

32、1800ln(412%45exp()(2 00102 0 0020 0110210021 0021 0210101rrrtrrrrrrrrrrrvrrtrtrtrtrvCtp附錄四：遺傳算法程序function Genetic(AimFunc)% This is simple genetic algorithm(SGA)% In this function ,it fulfils genetic algorithm%-These can be modified as you like-maxgen=200; % maximum generationsizepop=100; % size of

33、population % AimFunc=strimFunc; % this is function of counting fitnessfselect=tournament; % method of select % you can choose tournament;roulettefcode=float; % method of coding % you can choose float;grey;binary pcross=0.6; % probablity of crossover,between 0 and 1fcross=float; % method of crossover

34、 % 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=1 ; % length of bit of every variblebound=0 100000; %-individuals=struct(fitness,zeros(1,sizepop),.%value,zeros(1,sizepop),. chr

35、om,); % structure of populationavgfitness=; % average fitness of population bestfitness=; % best fitness of populationbestchrom=; % chromosome of best fitness% inivitializationfor i=1:sizepop % produce new population at random individuals.chrom(i,:)=Code(lenchrom,fcode,bound); x=Decode(lenchrom,boun

36、d,individuals.chrom(i,:),fcode); individuals.fitness(i)=Aimfunc(x);end% find minimum value which is bestbestfitness bestindex=min(individuals.fitness);bestchrom=individuals.chrom(bestindex,:);avgfitness=sum(individuals.fitness)/sizepop;% record average and best fitness of every generationtrace=avgfi

37、tness bestfitness; % evolution beginfor i=1:maxgen % selection individuals=Select(individuals,sizepop,fselect); avgfitness=sum(individuals.fitness)/sizepop; % crossover individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,. sizepop,fcross,i maxgen); % mutation individuals.chrom=Mutation(pmutati

38、on,lenchrom,individuals.chrom,. sizepop,fmutation,i maxgen,bound); % calculate fitness for j=1:sizepop x=Decode(lenchrom,bound,individuals.chrom(j,:),fcode); individuals.fitness(j)=Aimfunc(x); end % substitute chromosome of worest fitness % find minimum value which is best newbestfitness,newbestinde

39、x=min(individuals.fitness); worestfitness,worestindex=max(individuals.fitness); % substitute chromosome of worest fitness if bestfitnessnewbestfitness bestfitness=newbestfitness; bestchrom=individuals.chrom(newbestindex,:); end individuals.chrom(worestindex,:)=bestchrom; individuals.fitness(worestin

40、dex)=bestfitness; avgfitness=sum(individuals.fitness)/sizepop; trace=trace;avgfitness bestfitness; end % draw fitness of every generationhfig=findobj(Tag,trace);% See if it 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,tr

41、ace(:,2),b-);title(適應度曲線 終止代數 num2str(maxgen);xlabel(進化代數);ylabel(適應度);legend(平均適應度,最佳適應度);disp(適應度 變量);x=Decode(lenchrom,bound,bestchrom,fcode);% show in command windowvpa(bestfitness,10);vpa(x,10);disp(bestfitness x);function ret=AimFunc(x)% 求最小值%ret=sum(x);% 求最小值z=quadl(f,0.0001,7200,x);%求最大值%ret

42、=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 chromosome% lenchrom input : length of chromosome% opts input : tag of coding method% bound input : boundary of varibles% ret outpu

43、t: 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 coding 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)

44、; 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;endfunction ret=Cross(pcross,lenchrom,chrom,sizepop,opts,pop)% In this function,it fulfils a crossover among chromosomes% pcorss input : pr

45、obability of crossover% 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 chromosomeswit

46、ch 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)=index(2) pick=rand(1,2); index=ceil(pick.*sizepop); end % probability of crossover pick=rand; if pickpcross continue; end % random p

47、osition 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(index(1)=chrom(index(1)-tail1+tail2; chrom(index(2)=chrom(index(2)-tail2+tail1; end ret=chrom; case uniform % uniform cross for

48、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; while pick=0 pick=rand; end if pickpcross continue; end % random position of crossover pick=rand; while pick=0 pick=rand; end mask=2

49、ceil(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,bitget(chrom(index(2),j); chrom2=bitset(chrom2,. j,bitget(chrom(index(1),j); end end chrom(index(1)=chrom1; chrom(index(2)=chrom2;

50、 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 position of crossover pick=rand; while pick=0 pick=rand; end if pickpcross continue; end % random position of crossover pick=rand; whi

51、le 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),pos)=pick*v1+(1-pick)*v2; end ret=chrom;endfunction ret=Decode(lenchrom,bound,code,opts)% In this function ,it decode chromosome %

52、 lenchrom input : length of chromosome% opts input : tag of coding method% bound input : boundary of varibles% ret output: value of variblesswitch opts case binary % binary coding for i=length(lenchrom):-1:1 data(i)=bitand(code,2lenchrom(i)-1); code=(code-data(i)/(2lenchrom(i); end ret=bound(:,1)+da

53、ta./(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(lenchrom):-1:1 data(i)=bitand(code,2lenchrom(i)-1); code=(code-data(i)/(2lenchrom(i); end ret=bound(:,1)+data./(2.lenchrom-1).*(b

54、ound(:,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 fulfils a mutation among chromosomes% pcorss input : probability of mutation% lenchrom input : length of a chromosome% chrom input :

55、 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