安徽工程大學安徽工程大學 數學建模（選修課）課程論文數學建模（選修課）課程論文 題目：水庫排污問題題目：水庫排污問題 摘要：摘要： 本文主要需解決的問題是：本文主要需解決的問題是：討論由于此次事故 的發生，干流發生大面積污染的可能性，以及如何避 免該事故的發生。 針對問題一，我們只要考慮在事故發生到關閉水 庫的兩個小時內，流出水庫的污染物的質量小于 q 3 2 噸。 針對問題二，我們只需單獨計算1,2號水庫的污染 情況，然后再相加就可以了。 針對問題三，建立人工水渠就是在問題二的基礎 上使水庫1和水庫2產生聯系，我們只需考慮從水庫1 到水庫2的情況。 針對問題四，干流已經出現了大面積污染，在一 些外在條件短時間內無法干擾的情況下，我們可以采 取稀釋的辦法，是單位體積內的化學物質的濃度少于 危險警戒值時的濃度，這樣就可以短時間內控制污染。 關鍵詞：控制 人工水渠 稀釋 姓名：殷競存 專業：數學與應用數學 班級：數學 112 學號：3110801234 指導老師：周金明指導老師：周金明 成績：成績： 完成日期：完成日期：2013.6.272013.6.27 一：問題重述一：問題重述 1、 問題的背景 某條江流上有 2 條支流，每條支流上都興建了規 模相當的水庫。由于正處在雨水多發季節，因此兩個 水庫都以一定規模的流量進行泄洪。某天晚上 10：00，在其中的一個水庫中發生了兩船相撞的事故， 而其中的一條船裝載的 p 噸化學物質（這里的化學物 質可以是具有揮發性的，也可能是急難揮發的）全部 泄漏至水庫中。當水上航運事故處置中心接獲事故報 告，立即要求該水庫關閉水庫泄洪閘，以免化學物質 隨洪水流入干流，發生更大規模的污染。水庫閘門開 始關閉時，已經處在事故發生后的 1 個小時，而水庫 閘門徹底關閉也需要 1 個小時的時間。 根據當地環 境監測的有關規定，干流大面積污染的危險警戒值設 為：三小時內 q 噸該化學物質發生泄漏。 2、面臨的問題 (1) 試建立合理的數學模型，討論由于此次事故 的發生，干流發生大面積污染的可能性； (2) 如果在另外的一水庫中有一化工廠違規排放 廢料。廢料中同樣含有該化學物質。該工廠為躲避環 境監測站的監控，均在晚上 9：00-12：00 違規進行 周期性排放。在這種情形下，討論由于此次事故的發 生，干流發生大面積污染的可能性； (3) 如果以上兩個水庫間有一條人工修建的水渠 相連接，水渠中的水流流向不定，但保證兩水庫之間 的水流能夠相互影響。那么上述結果是否會改變？請 給出說明，若有改變，則給出修正的模型及結果； (4) 如果發生了大面積污染，那么針對第三種情 況，試給出在短時間內控制污染模型。 二：問題的假設 （1）污染物為速溶物質，因此藥品從船上流入水 中的時間很少，可以忽略不計； （2）污染物質從水庫中一經流出就進入干流； （3）水庫和河流中的水流都是處于推流狀態； （4）兩水庫事故發生條件相同，即兩水庫有相同 的客觀條件； （5）被污染的水庫關閉泄洪閘后不再有水流流入 干流； （6）不考慮生物等因素在水庫泄洪過程中的作用， 污染物除了流出外不因腐爛沉積等手段從水中消 失； （7）外界因素不對水庫的體積變化產生影響，例 如：雨水、地表徑流、底下徑流等； （8）參與模型的變量是連續變化的，并且充分光 滑； （9）不考慮從不同的渠道流入與流出水庫之間的 區別，只考慮攜帶污染物的水流入水庫和水庫中 的水流出對水庫污染程度的影響，因此可以把水 庫看成是單流入單流出的系統。 三：符號的約定三：符號的約定 （1）：t 時刻水庫水的流入速度；)( 1 tr （2）：t 時刻流入水庫的污染物的濃度；)( 1 tp （3）：t 時刻水庫水的流出的速度；)( 0 tr （4）：t 時刻流出水庫的污染物的濃度；)( 0 tp （5）：t 時刻水庫中污染物的濃度；)(tp （6）：t 時刻水庫水的體積；)(tv （7）：計算體積元內該污染物的增量； p m （8） ： 為時間；t （9） ： 為從水庫中流出的水中的污染物的濃度；C （10） ：為水庫的流出速度，即流量；Q （11）：為泄洪閘處到污染處的距離； X （12） Q：為水庫的流量； （13） q：為排入河流的污水的流量； （14） ：為河流中污染物的本底濃度； 1 C （15）:為水庫中的污染物的濃度； 2 C 四：模型的分析與建立四：模型的分析與建立 由問題的分析中知道，流入水庫的污染物能以很 快的速度與水庫中的水均勻混合，也就是說水庫中的 污染狀況在任何局部水體都是一樣的，污染程度與水 體在水庫中的位置無關，因此我們可以建立下面模型。 模型一 問題一：此次事故的發生，干流發生大面積污染的可 能性； 根據物質平衡原理和題目假設可知： 水庫 1 中污染物的改變量 = 流入的污染物的量 流出污染物的量 于是對于充分小的，在時間（t，t+）內有：tt ttrtptrtptvtpttvttp)()()()()()()()( 0011 兩邊同除，并使0 得：tt (1) )()()()( )()()()( 0011 trtptrtp t tvtpttvttp 現假設 f(t)=p(t)v(t)得： dt d tv dt d tp dt d dt d t tfttf tptvtvtptf)()()()()( )()( )()( 即原式可以寫為： (2) )()()()()()( 0011 )()( trtptrtp dt d tv dt d tp tptv 在水庫 1 中發生撞船事故后，污染物處于非穩定 排放即：，而由于水庫閘門的關閉也勢必會引0 )( dt d tp 起水庫中水的體積變化，故：。現不考慮流入0 )( dt d tv 水庫中的水所含有與泄漏污染物相同物質的情況而帶 來的影響，即可看作，另外由問題分析中知道：0)( 1 tp 流出的污染物的濃度應與水庫中污染物濃度相同，即 這樣對于問題一我們可以得到求解公式：)()( 0 tptp (3) )()()()( 0 )()( trtp dt d tv dt d tp tptv 進一步我們假設從水庫中流出的水的流量初始值 （從 t=0 時算起）為 ，在關閉閘門的過程中，我們 0 r 假定流量處于線性變化的趨勢。這一假設是基于流量 與過流面積為線性關系上作出的，進一步可得： （4) )72003600.(. 3600 2 )36000.(. )( 0 0 0 0 t tr r tr tr 從上面可看出為分數函數，這主要是因為水庫)( 0 tr 閘門關閉是在事故發生一小時后作出的。現在有了 的變化的表達式，為了能求出的表達式。我們)( 0 tr)(tp 還要寫出的表達式。首先我們假設水庫的體積的)(tv 初始值為 （t=0 時） ， 值我們可以通過衛星定位系 0 v 0 v 統及所建立的模型求出（具體衛星定位系統模型見附 表） 。而跟相關的還有的值。我們假設為一)(tv)( 1 tr)( 1 tr 定值 ，則隨隨時間變化的關系式為： 1 r)(tv ttrtrvtv)()( 010 由于為分段函數可知：也響應的為分段函數，)( 0 tr)(tv 具體函數表達式為： (5) )72003600.(.) 3600 2( )36000.(. )( 0 010 010 tt tr rtrv ttrtrv tv 把（4）式代入（3）式可以得到： 當秒時：36000 t (6) 0 )( 01001 )()()()(rtp dt d trtrvrrtp tp 當秒時：72003600 t (7) ) 3600 2()() 3600 2() 1800 2()( 0 0 )( 0 010 0 01 tr rtp dt d t tr rtrvt r rrtp tp 對（5）式化簡有： (8) 1 )( 010 )()(rtp dt d trtrv tp 通過推導的出： (9) cetp trtrv rr r )ln( 1 010 01 0 )( 有已知條件可知：，故 0 /)0(vPp10 1 00 /rr r vvpc 經簡化后： 00010 /)()( 01 1 01 1 vpvtrrvtprr r rr r (10)36000 t 對（7）簡化后得； ( 20 010 0 20 010 1 3600 )2( ) 3600 ( 3600 )2( )()(ln t v trrv dttr t r trrv dt rtp 11) 設，(r1-2r0)=b，得：cv 0 a r 3600 0 當acb4 2 )ln( 2 1 ) 3600 () 4 2 arctan( 4 2 )()(ln 20 22 ctbta a r bac bat bac rtp (12) ) 4 2 arctan( 4 2 2 22 bac bat bac a b 設：) 4 2 arctan( 4 2 22 bac bat bac z 最后得到： (13) 1 )72005.0()ln(5.0 1 2 )(cetp zabcbtat 保持連續性，當 t=3600 時，=，此時可)3600( 1 p)3600(p 得到相應的 值，但由于不確定因素很多，故確定 1 C 不是很容易，這主要是缺少數據造成的。 1 C 當時得到（14）式acb4 2 1 )7200 2 ()ln( 2 1 1 2 )(Cetp za b cbtat （14） c acbbat acbbat acb z 42 42 ln 4 1 2 2 2 同樣為保持連續性，要求當 t=3600 時，=)3600( 1 p 。)3600(p 最后： )(72003600.4. )(72003600.4. )(36000.)( )( 2 1 )7200 2 ()ln( 2 1 2 1 )72005 . 0()ln(5 . 0 0010 2 2 01 1 01 1 stacbce stacbce stvtrrv tp za b cbtat zabcbtat rr r rr r 當 當 （15） 那么時間內流出水庫的污染物的量便可表示為：t ttrtpQ)()( 0 （16） 在（s）時流出的污染物的量為：36000 t t drvpvtrrvQrr r rr r 3600 0 0000101 /)( 01 1 01 1 （17） 在（s）流出的量為：72003600 t 7200 3600 0 02 ) 3600 2()(dt tr rtpQ （18） 流出的總量： 21 QQQ （19） 再用 Q 與 2/3q 進行比較，便得出是否會發生大面積 污染。 問題二：如果在另外的一水庫中有一化工廠均在晚上 9：00-12：00 違規進行周期性排放同樣含有該化學 物質的廢料。討論由于此次事故的發生，干流發生大 面積污染的可能性。 由題目分析中知道，水庫 2 中有一化工廠違規排 放含有該污染物的廢料，而由于兩個水庫之間沒有聯 系，故我們只需要單獨考慮水庫 2 的排污量，然后加 上水庫 1 的污染物排放量，最后綜合考慮兩個水庫所 排污染物的總量對干流的影響就可以了。下面我們將 在水庫 1 模型的基礎上建立 2 水庫的模型。 )()()()()()( 0011 )()( trtptrtp dt d tv dt d tp tptv （20） 設（常數） ，即在 9：0012：00 這段時間內ktp)( 1 污染物以一個恒定值流入水庫，考慮水庫水的流入速 度為一定值 ，流出速度也為一個定值 1 r 0 r 這樣水庫體積的表達式可寫成如下的公式： trtrvtv 010 )( )108000( t （21） 代入上面的表達式可簡化為： 01 )( 010 )()(rtprk dt d trrv tp （22） 用上式可求出，其中可通 2010 21 1 )()(Ctrrvktp rr r 2 C 過求得。那么從 9：0012：00 這三個小時內流kp)0( 出閘門的污染物的總量就可以求出來，污染物的量為： 10800 0 0 )(dtrtpQ （23） QQQ （24） （ 為三小時內從 1 水庫和 2 水庫流出的污染物總 Q 量） 若則發生大面積污染；qQ 若則不會發生大面積污染；qQ 問題三：如果兩個水庫間有一條人工修建的水渠相連 接，水渠中的水流流向不定，但保證兩水庫之間的水 流能夠相互影響。那么問題二結果是否會改變？由題 目分析可知，當兩水庫之間有一人工修建的水渠相互 連通時，水渠中的水流必然是從高水位流向低水位， 為了使模型簡化，我們有以下說明： 1、由于兩水庫是連通的，因此在水庫 1 關閘前， 兩水庫的液面必然是趨近于等高的，否則必然有水從 一個水庫流向另一個水庫。 2、在事故發生后一個小時內，考慮兩水庫的規模 相當，且水庫 1 沒有關閉泄洪閘，此時兩水庫彼此不 受影響。 3、在事故發生一個小時后，考慮到水庫 1 要關閉 泄洪閘，這勢必引起水庫 1 的水位上升，由說明 1 可 知水庫 1 中的水必將通過水渠流向水庫 2。 從上面的分析可知，在有連通水渠的情況下，我 們只需考慮從 1 水庫向 2 水庫的流入情況，而不必考 慮從 2 水庫向 1 水庫的流入情況。下面我們就該問題 給出進一步分析。 在第一個小時內：（即 9：0010：00） 此時，水庫 1 中還未發生事故，只有水庫 2 中在 排放廢料，用問題 2 的模型我們可求出流出的污染物 的量為： 3600 0 01 )(dtrtpQ （25） （由（20）式確定）)(tp 在第二個小時內（即 10：0011：00） 水庫 1 已發生裝船事故，而水庫 2 繼續排放廢料， 但泄洪閘還未關閉，所以我們仍然獨立考慮。 對水庫 1 我們用問題一的模型求解得： 3600 0 00001021 /)( 01 1 01 1 dtrvpvtrrvQrr r rr r （26） 對水庫 2 我們仍有問題二的模型求解得到： （由（20）確定） 7200 3600 022 )(dtrtpQ )(tp （27 22212 QQQ （28） 在第三個小時內（即 11：0012：00） ； 水庫 1 泄洪閘正在關閉，而水庫 2 繼續排放廢料。 這時水庫 1 和水庫 2 就要結合在一起考慮了，如下圖 所示： 圖 一 由上面分析可列出以下的方程式： 對于水庫 1： （29） )()()()()()( 00 )()( trtptrtp dt d tv dt d tp tptv 其中的表達式為：)(tv （）ttrt tr rtrvtv )() 3600 2()( 0 0 0101 72003600 t 而： 3600 2)( 0 00 tr rtr （）72003600 t 對水庫 2： )()()()()()()()( 0 0 1 1 )()( trtptrtptrtp dt d tv dt d tp tptv （30） 其中 trttrtrvtv 0 0 1 02 )()( 由基本假設知道，水庫 1 和水庫 2 規模相當則： trttrtrvtv 0 010 )()( 假定為一定值 。)( 0 tr 0 r 聯立解出：由于所得表達式非常復雜，我)()( tptp與 們把表達式放在附錄（三）里表示。而由于過)(tp)( tp 于復雜，我們這里就不給出解析解了。 （） 7200 3600 031 )()(dttrtpQ 72003600 t dtrtpQ 10800 7200 0 32 )( （31） 3600 2)( 0 10 tr rtr 故有： 32313 QQQ （32） 在 12：00 以后一個小時內，1 水庫完全關閉但 1 水庫的水將通過渠道流入 2 水庫內。則有以下式： 對水庫 1 有： )()()()( 0 )()( trtp dt d tv dt d tp tptv （33） 其中表達式為： )(tv trtrvtv 010 )( 對水庫 2 有： （34） )()()()()()( 0 0 )()( trtptrtp dt d tv dt d tp tptv 其中表達式為： )(tvtrtrtrvtv 0 1 00 )( 求得 dtrtpQ 14400 10800 0 42 )( （35） 42321 QQQQQ 總 （36） 若則干流不會發生大面積污染；qQ 若則干流會發生大面積污染；qQ 問題四：如果發生了大面積污染，那么針對第三種情 況，試給出在短時間內控制污染模型。 由問題的分析中知道，我們只能通過稀釋原理， 建立污染物在干流中遷移遷移模型，使干流水體計算 體積元內該污染物的增量為負值時，從而使得在 p m 干流水體污染物的濃度低于危險警戒值時的濃度，才 能在短時間內達到控制污染。 模型二 為使模型清楚，我們先給出下面所用名詞解釋。 源和漏：是對體積元內污染物變化的一種描述， 源是體積元內污染物的增加速率，漏是體積元內污染 物的減少速率。這里的“漏”不意味著漏掉，而有更 廣泛的意義，如污染物的降解、沉淀和揮發等都屬于 “漏” 。 由于污染物在干流中遷移符合遷移方程，由基本 假設中我們知道河流中的水流處于推流狀態，也就是 體積元中水分子以同一速度向下游運動如圖二和圖三。 設 C(x)和 C(x+)分別為進入水片的水中某中污染物x 的濃度，是計算體積元中所含的該污染物的質量。 p m 按照質量守衡原理，計算體積元里污染物質量的增量 為： txbStxStxxQxxCxQxCm slp )()()()( （1） （2） （3） （4） txASv （37） （5） 上式中（1）至（5）各項的意義如下 （1）儲存量項； )(CAxmp （38） （2）平流輸送項； （3）側向的源和漏，其中 是單位時間內單位長 l S 度上的源和漏，它通常是由側向分布流量 帶入的，q 因此可把 寫為： l S il CqS 這里 是該分布流量中污染物的濃度； i C （4）表面的源和漏， 是單位時間內單位面積 s S 上的源和漏（） ；smg 2 / （5）體積元內的源和漏， 是單位時間內單位 v S 體積內的源和漏（） ；smg 3 / 圖 二 圖 三 把（38）式帶入（37）式，用除方程兩邊并令tx 和，則得：0x0t vsl SASbS x QC t AC )()( （39） 式（10）就是推流時的污染物遷移方程，以后將廣泛 地利用這個方程，在求解問題四時，我們只要 、 l S 和 的值為負數就可以了，即可以通過關閉人工水 s S v S 渠和部分關閉水庫 2 泄洪閘的手段實現短時間內控制 污染。 五：模型的評價與改進 5.1：模型的優點 本模型在建立模型前，由于水體運動非常復雜， 污染物的性質也多種多樣，污染物在水中的遷移情況 受各種因素限制，不容易全部考慮。本文綜合考慮上 述因素，先通過合理的假設，首先從易溶急難揮發性 污染物入手，考慮污染物和庫水混合均勻的情況，利 用物質平衡原理，針對問題中各種情況建立物質平衡 方程模型。在此基礎上考慮污染物的一維遷移模型。 在建立一維遷移擴展模型時，我們通過轉化思想， 將污染物對干流形成大面積污染的可能轉換為求解在 關閉泄洪閘的過程中，水庫中多大面積水域內發生事 故時才會對干流形成大面積污染，使得問題求解方案 明了，為了驗證模型的正確性，我們通過遺傳算法對 事例進行求解，很好的驗證了模型。 5.2：模型的缺點 模型在建立過程中我們所考慮的此污染物為易溶、 不揮發的物質，而沒有考慮污染物為油性、揮發性和 沉淀物質。另外，我們在考慮水體流動時沒有考慮水 體上、中、下層的不同流動情況。在問題四的求解時， 我們只考慮推流問題，而實際污染隨水流遷移是以水 團形式的，這樣我們就需要考慮水體不同水層的流動 問題。時間有限，未能對其他方面做深入的探討，使 得模型十分完美。 5.3：模型的改進 從問題的分析及假設可以知道，我們的模型在建 立過程中我們所考慮的只是一種污染物，即此污染物 為易溶、不揮發的物質。實際上污染物可能為油性、 揮發性和沉淀物質。對于油性物質我們應該考慮該物 質的密度是否大于水的密度，若油性污染物小于水的 密度，因為油性物質不溶于水，該污染物應該是漂浮 在水面上的，而且呈油膜狀，因此我們就應考慮水庫 表層水的流動情況；若油性污染物的密度大于水的密 度，則污染物將會沉如水底，那么我們就應該考慮深 沉水體的水流問題；若油性污染物等于水的密度，那 么污染物就會懸浮在水庫水體中，因此我們就應該考 慮整個水體上、中、下層的流動情況（因為不同水層 的水流的流動情況是不同的） 。對于揮發性的污染物， 我們還需要考慮此污染物在水和空氣中的分配比例。 在考慮問題四的時候我們考慮的是推流問題，實際污 染隨水流遷移是以水團形式的，這樣我們就需要考慮 水體不同水層的流動問題而。為使模型更符合實際， 模型必須在上述方面在做深入探討。 參考文件： 【1】應玖茜，魏權齡。非線性規劃與理論。北京： 中華人民大學出版社，1994 【2】岳天祥。資源環境數學模型手冊。北京：科學 出版社，2003 【3】陳春云，鄭彤。環境系統數學模型。北京：化 學工業出版社，2003 【4】周建華，黃燕 。MATLAB5.3。北京：北京大學 出版社，2000 【5】趙靜，但琦。數學建模與數學實驗。北京：高 等教育出版社，2007 附錄一： 水庫實際庫容及水域面積計算公式： )( i hHsV snS 式中，s 為 DEM 象元面積；H 為水庫水位高程； 為 i h 庫域 DEM 各像元高程值；n 為庫域 DEN 像元數。 （曾永年，馬海州，沙占江等；龍羊峽庫區環境動態 監測信息系統的建立與應用。遙感學報，2000.4） 附錄二： 序號項目單位數據相應時間 1 入庫水量 億立方 米 469.072 2 入庫平均流量 立方米/ 秒 1436 3 出庫水量 億立方 米 494.287 4 出庫平均流量 立方米/ 秒 1513 5 最大入庫流量 立方米/ 秒 28100 2003-9-8 8:00 6 最小入庫流量 立方米/ 秒 0 2004-2-20 20:00 7 最大出庫流量 立方米/ 秒 12000 2003-9-3 2:00 8 最小出庫流量 立方米/ 秒 427 2004-7-10 8:00 9 最高庫水位米 156.98 2003-10-17 14:00 10 最低庫水位米 138.3 2004-7-16 8:00 11 大廠發電量 億千瓦 時 44.6345 12 平均出力萬千瓦 49.20 13 最高負荷萬千瓦 109.62003-10-26 14 最低負荷萬千瓦 11.22004-7-25 15 最大調峰萬千瓦 69.82003-10-2 16 最小調峰萬千瓦 42.42004-9-12 17 平均調峰萬千瓦 32.68 18 調峰電量 億千瓦 時 .0912 19 棄水量億立方 119.946 米 20 平均耗水率 立方米/ 千瓦時 7.79 21 2x2 發電量 億千瓦 時 3.4254 22 大、小廠總電 量 億千瓦 時 48.0599 附錄一： 水庫實際庫容及水域面積計算公式： )( i hHsV snS 式中，s 為 DEM 象元面積；H 為水庫水位高程； 為 i h 庫域 DEM 各像元高程值；n 為庫域 DEN 像元數。 （曾永年，馬海州，沙占江等；龍羊峽庫區環境動態 監測信息系統的建立與應用。遙感學報，2000.4） 附錄二： 序號項目單位數據相應時間 1 入庫水量 億立方 米 469.072 2 入庫平均流量 立方米/ 秒 1436 3 出庫水量 億立方 米 494.287 4 出庫平均流量 立方米/ 秒 1513 5 最大入庫流量 立方米/ 秒 28100 2003-9-8 8:00 6 最小入庫流量 立方米/ 秒 0 2004-2-20 20:00 7 最大出庫流量 立方米/ 秒 12000 2003-9-3 2:00 8 最小出庫流量 立方米/ 秒 427 2004-7-10 8:00 9 最高庫水位米 156.98 2003-10-17 14:00 10 最低庫水位米 138.3 2004-7-16 8:00 11 大廠發電量 億千瓦 時 44.6345 12 平均出力萬千瓦 49.20 13 最高負荷萬千瓦 109.62003-10-26 14 最低負荷萬千瓦 11.22004-7-25 15 最大調峰萬千瓦 69.82003-10-2 16 最小調峰萬千瓦 42.42004-9-12 17 平均調峰萬千瓦 32.68 18 調峰電量 億千瓦 時 .0912 19 棄水量 億立方 米 119.946 20 平均耗水率 立方米/ 千瓦時 7.79 21 2x2 發電量 億千瓦 時 3.4254 22 大、小廠總電 量 億千瓦 時 48.0599 附錄三 ) 1% 306030 30 1 arctan(:2% 90036003600180036009002:1% )1%/ )2%152%30 1%)1800360018001800ln( 4 1 2%45exp()( 2 0010 2 0 00 2 0 0110 2 100 2 1 00 2 1 0 2 10101 rrrtr rrrrrrrrrrv rr trtrtrtrvCtp 附錄四： 遺傳算法程序 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 generation sizepop=100; % size of population % AimFunc=strimFunc; % this is function of counting fitness fselect=tournament; % method of select % you can choose tournament;roulette fcode=float; % method of coding % you can choose float;grey;binary pcross=0.6; % probablity of crossover,between 0 and 1 fcross=float; % method of crossover % you can choose float;simple;uniform pmutation=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 varible bound=0 100000; %- individuals=struct(fitness,zeros(1,sizepop),.%value,zeros(1,sizepop),. chrom,); % structure of population avgfitness=; % average fitness of population bestfitness=; % best fitness of population bestchrom=; % chromosome of best fitness % inivitialization for 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 best bestfitness bestindex=min(individuals.fitness); bestchrom=individuals.chrom(bestindex,:); avgfitness=sum(individuals.fitness)/sizepop; % record average and best fitness of every generation trace=avgfitness bestfitness; % evolution begin for 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(pmutation,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,newbestindex=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(worestindex)=bestfitness; avgfitness=sum(individuals.fitness)/sizepop; trace=trace;avgfitness bestfitness; end % draw fitness of every generation hfig=findobj(Tag,trace); % See if it is open if ishandle(hfig) figure(hfig); else hfig=figure(Tag,trace); end figure(hfig); r c=size(trace); plot(1:r,trace(:,1),r-,1:r,trace(:,2),b-); title(適應度曲線 終止代數 num2str(maxgen); xlabel(進化代數);ylabel(適應度); legend(平均適應度,最佳適應度); disp(適應度 變量); x=Decode(lenchrom,bound,bestchrom,fcode); % show in command window vpa(bestfitness,10); vpa(x,10); disp(bestfitness x); function ret=AimFunc(x) % 求最小值 %ret=sum(x); % 求最小值 z=quadl(f,0.0001,7200,x); ret=abs(z-49.771090910000000000000000000000); %求最大值 %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 chromosome % lenchrom input : length of chromosome % opts input : tag of coding method % bound input : boundary of varibles % ret output: chromosome switch 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); 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; end function 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 method of crossover % pop input : current serial number of generation and maximum gemeration % ret output : new set of chromosome switch 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 position of crossover pick=rand; while pick=0 pick=rand; end pos=ceil(pick.*sum(lenchrom); tail