B水庫排污問題數(shù)學(xué)112夏林龍

1、安徽工程大學(xué) 數(shù)學(xué)建模 課程設(shè)計論文水庫排污問題指導(dǎo)老師： 周金明 成 績： 完成日期： 2013年7月3日摘要本文針對水庫突發(fā)性事故排污問題，首先通過建立二維水質(zhì)污染物濃度模型，給出了單個水庫對干流造成大面積污染的可能性Q-（W1-W2）/(2/3*Q*1000);然后建立兩水庫排污模型，分析了在另一水庫有連續(xù)點源污染物排放及水流相互影響的情況下，兩水庫對干流造成大面積污染的可能性大小Q-(W1+W2+W3)/(q*1000);并進一步針對第三種情況的發(fā)生，給出在短時間內(nèi)控制污染的有效措施；且討論了若污染物具有揮發(fā)性，上述各情況造成干流發(fā)生大面積污染的可能性大小，為水庫事故性排污問題提供了有

2、價值的理論依據(jù)。一、問題重述問題的背景某條江流上有2條支流，每條支流上都興建了規(guī)模相當(dāng)?shù)乃畮臁Ｓ捎谡幵谟晁喟l(fā)季節(jié)，因此兩個水庫都以一定規(guī)模的流量進行泄洪。某天晚上10：00，在其中的一個水庫中發(fā)生了兩船相撞的事故，而其中的一條船裝載的p噸化學(xué)物質(zhì)（這里的化學(xué)物質(zhì)可以是具有揮發(fā)性的，也可能是急難揮發(fā)的）全部泄漏至水庫中。當(dāng)水上航運事故處置中心接獲事故報告，立即要求該水庫關(guān)閉水庫泄洪閘，以免化學(xué)物質(zhì)隨洪水流入干流，發(fā)生更大規(guī)模的污染。水庫閘門開始關(guān)閉時，已經(jīng)處在事故發(fā)生后的1個小時，而水庫閘門徹底關(guān)閉也需要1個小時的時間。 根據(jù)當(dāng)?shù)丨h(huán)境監(jiān)測的有關(guān)規(guī)定，干流大面積污染的危險警戒值設(shè)為：三小時內(nèi)q

3、噸該化學(xué)物質(zhì)發(fā)生泄漏。2、面臨的問題(1) 試建立合理的數(shù)學(xué)模型，討論由于此次事故的發(fā)生，干流發(fā)生大面積污染的可能性；(2) 如果在另外的一水庫中有一化工廠違規(guī)排放廢料。廢料中同樣含有該化學(xué)物質(zhì)。該工廠為躲避環(huán)境監(jiān)測站的監(jiān)控，均在晚上9：00-12：00違規(guī)進行周期性排放。在這種情形下，討論由于此次事故的發(fā)生，干流發(fā)生大面積污染的可能性；(3) 如果以上兩個水庫間有一條人工修建的水渠相連接，水渠中的水流流向不定，但保證兩水庫之間的水流能夠相互影響。那么上述結(jié)果是否會改變？請給出說明，若有改變，則給出修正的模型及結(jié)果；(4) 如果發(fā)生了大面積污染，那么針對第三種情況，試給出在短時間內(nèi)控制污染模型

4、。 二、問題的假設(shè)污染物為速溶物質(zhì)，因此藥品從船上流入水中的時間很少，可以忽略不計；污染物質(zhì)從水庫中一經(jīng)流出就進入干流；水庫和河流中的水流都是處于推流狀態(tài)；兩水庫事故發(fā)生條件相同，即兩水庫有相同的客觀條件；被污染的水庫關(guān)閉泄洪閘后不再有水流流入干流；不考慮生物等因素在水庫泄洪過程中的作用，污染物除了流出外不因腐爛沉積等手段從水中消失；外界因素不對水庫的體積變化產(chǎn)生影響，例如：雨水、地表徑流、底下徑流等；參與模型的變量是連續(xù)變化的，并且充分光滑；不考慮從不同的渠道流入與流出水庫之間的區(qū)別，只考慮攜帶污染物的水流入水庫和水庫中的水流出對水庫污染程度的影響，因此可以把水庫看成是單流入單流出的系統(tǒng)。三

5、、符號說明（1）：t時刻水庫水的流入速度；（2）：t時刻流入水庫的污染物的濃度；（3）：t時刻水庫水的流出的速度；（4）：t時刻流出水庫的污染物的濃度；（5）：t時刻水庫中污染物的濃度；（6）：t時刻水庫水的體積；（7）：計算體積元內(nèi)該污染物的增量；（8） ： 為時間；（9） ： 為從水庫中流出的水中的污染物的濃度；（10） ：為水庫的流出速度，即流量；（11）：為泄洪閘處到污染處的距離；（12） Q：為水庫的流量；（13） q：為排入河流的污水的流量；（14）：為河流中污染物的本底濃度；（15）:為水庫中的污染物的濃度； 注意：部分符號見論文中說明。四、模型的建立與求解由問題的分析中知道，流

6、入水庫的污染物能以很快的速度與水庫中的水均勻混合，也就是說水庫中的污染狀況在任何局部水體都是一樣的，污染程度與水體在水庫中的位置無關(guān)，因此我們可以建立下面模型。模型一問題一：此次事故的發(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ì)的

7、情況而帶來的影響，即可看作，另外由問題分析中知道：流出的污染物的濃度應(yīng)與水庫中污染物濃度相同，即這樣對于問題一我們可以得到求解公式： （3）進一步我們假設(shè)從水庫中流出的水的流量初始值（從t=0時算起）為，在關(guān)閉閘門的過程中，我們假定流量處于線性變化的趨勢。這一假設(shè)是基于流量與過流面積為線性關(guān)系上作出的，進一步可得： （4）從上面可看出為分?jǐn)?shù)函數(shù)，這主要是因為水庫閘門關(guān)閉是在事故發(fā)生一小時后作出的?，F(xiàn)在有了的變化的表達式，為了能求出的表達式。我們還要寫出的表達式。首先我們假設(shè)水庫的體積的初始值為（t=0時），值我們可以通過衛(wèi)星定位系統(tǒng)及所建立的模型求出（具體衛(wèi)星定位系統(tǒng)模型見附表）。而跟相關(guān)的還

8、有的值。我們假設(shè)為一定值，則隨隨時間變化的關(guān)系式為： 由于為分段函數(shù)可知：也響應(yīng)的為分段函數(shù)，具體函數(shù)表達式為： （5）把（4）式代入（3）式可以得到：當(dāng)秒時： （6）當(dāng)秒時： （7）對（5）式化簡有： （8）通過推導(dǎo)的出： （9）有已知條件可知：，故經(jīng)簡化后： 時 （10）對（7）簡化后得；（11）設(shè)，(r1-2r0)=b，得：當(dāng) （12）設(shè)：最后得到： （13）保持連續(xù)性，當(dāng)t=3600時，=，此時可得到相應(yīng)的值，但由于不確定因素很多，故確定不是很容易，這主要是缺少數(shù)據(jù)造成的。當(dāng)時得到（14）式 （14）同樣為保持連續(xù)性，要求當(dāng)t=3600時，=。最后： （15）那么時間內(nèi)流出水庫的污染物

9、的量便可表示為： （16）在（s）時流出的污染物的量為： （17）在（s）流出的量為： （18）流出的總量： （19）再用Q與2/3q進行比較，便得出是否會發(fā)生大面積污染。問題二：如果在另外的一水庫中有一化工廠均在晚上9：00-12：00違規(guī)進行周期性排放同樣含有該化學(xué)物質(zhì)的廢料。討論由于此次事故的發(fā)生，干流發(fā)生大面積污染的可能性。由題目分析中知道，水庫2中有一化工廠違規(guī)排放含有該污染物的廢料，而由于兩個水庫之間沒有聯(lián)系，故我們只需要單獨考慮水庫2的排污量，然后加上水庫1的污染物排放量，最后綜合考慮兩個水庫所排污染物的總量對干流的影響就可以了。下面我們將在水庫1模型的基礎(chǔ)上建立2水庫的模型。

10、（20）設(shè)（常數(shù)），即在9：0012：00這段時間內(nèi)污染物以一個恒定值流入水庫，考慮水庫水的流入速度為一定值，流出速度也為一個定值這樣水庫體積的表達式可寫成如下的公式： （21）代入上面的表達式可簡化為： （22） 用上式可求出，其中可通過求得。那么從9：0012：00這三個小時內(nèi)流出閘門的污染物的總量就可以求出來，污染物的量為： （23） （24）（為三小時內(nèi)從1水庫和2水庫流出的污染物總量）若則發(fā)生大面積污染；若則不會發(fā)生大面積污染；問題三：如果兩個水庫間有一條人工修建的水渠相連接，水渠中的水流流向不定，但保證兩水庫之間的水流能夠相互影響。那么問題二結(jié)果是否會改變？由題目分析可知，當(dāng)兩水庫

11、之間有一人工修建的水渠相互連通時，水渠中的水流必然是從高水位流向低水位，為了使模型簡化，我們有以下說明：1、由于兩水庫是連通的，因此在水庫1關(guān)閘前，兩水庫的液面必然是趨近于等高的，否則必然有水從一個水庫流向另一個水庫。2、在事故發(fā)生后一個小時內(nèi)，考慮兩水庫的規(guī)模相當(dāng)，且水庫1沒有關(guān)閉泄洪閘，此時兩水庫彼此不受影響。3、在事故發(fā)生一個小時后，考慮到水庫1要關(guān)閉泄洪閘，這勢必引起水庫1的水位上升，由說明1可知水庫1中的水必將通過水渠流向水庫2。從上面的分析可知，在有連通水渠的情況下，我們只需考慮從1水庫向2水庫的流入情況，而不必考慮從2水庫向1水庫的流入情況。下面我們就該問題給出進一步分析。在第一

12、個小時內(nèi)：（即9：0010：00）此時，水庫1中還未發(fā)生事故，只有水庫2中在排放廢料，用問題2的模型我們可求出流出的污染物的量為： （25）（由（20）式確定）在第二個小時內(nèi)（即10：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

13、）其中的表達式為： （）而： （）對水庫2： （30）其中 由基本假設(shè)知道，水庫1和水庫2規(guī)模相當(dāng)則：假定為一定值。聯(lián)立解出：由于所得表達式非常復(fù)雜，我們把表達式放在附錄（三）里表示。而由于過于復(fù)雜，我們這里就不給出解析解了。 （） （31）故有： （32）在12：00以后一個小時內(nèi)，1水庫完全關(guān)閉但1水庫的水將通過渠道流入2水庫內(nèi)。則有以下式：對水庫1有： （33）其中表達式為： 對水庫2有： （34）其中表達式為： 求得 （35） （36）若則干流不會發(fā)生大面積污染；若則干流會發(fā)生大面積污染；問題四：如果發(fā)生了大面積污染，那么針對第三種情況，試給出在短時間內(nèi)控制污染模型。由問題的分析中知道

14、，我們只能通過稀釋原理，建立污染物在干流中遷移遷移模型，使干流水體計算體積元內(nèi)該污染物的增量為負(fù)值時，從而使得在干流水體污染物的濃度低于危險警戒值時的濃度，才能在短時間內(nèi)達到控制污染。模型二為使模型清楚，我們先給出下面所用名詞解釋。源和漏：是對體積元內(nèi)污染物變化的一種描述，源是體積元內(nèi)污染物的增加速率，漏是體積元內(nèi)污染物的減少速率。這里的“漏”不意味著漏掉，而有更廣泛的意義，如污染物的降解、沉淀和揮發(fā)等都屬于“漏”。 由于污染物在干流中遷移符合遷移方程，由基本假設(shè)中我們知道河流中的水流處于推流狀態(tài)，也就是體積元中水分子以同一速度向下游運動如圖二和圖三。設(shè)C(x)和C(x+)分別為進入水片的水中

15、某中污染物的濃度，是計算體積元中所含的該污染物的質(zhì)量。按照質(zhì)量守衡原理，計算體積元里污染物質(zhì)量的增量為： （1） （2） （3） （4） + （37） （5）上式中（1）至（5）各項的意義如下（1）儲存量項； （38）（2）平流輸送項；（3）側(cè)向的源和漏，其中是單位時間內(nèi)單位長度上的源和漏，它通常是由側(cè)向分布流量帶入的，因此可把寫為：這里是該分布流量中污染物的濃度；（4）表面的源和漏，是單位時間內(nèi)單位面積上的源和漏（）；（5）體積元內(nèi)的源和漏，是單位時間內(nèi)單位體積內(nèi)的源和漏（）；圖 二圖 三把（38）式帶入（37）式，用除方程兩邊并令和，則得： （39）式（10）就是推流時的污染物遷移方程，以

16、后將廣泛地利用這個方程，在求解問題四時，我們只要、和的值為負(fù)數(shù)就可以了，即可以通過關(guān)閉人工水渠和部分關(guān)閉水庫2泄洪閘的手段實現(xiàn)短時間內(nèi)控制污染。五、模型的評價與改進5.1模型的優(yōu)點： 本模型在建立模型前，由于水體運動非常復(fù)雜，污染物的性質(zhì)也多種多樣，污染物在水中的遷移情況受各種因素限制，不容易全部考慮。本文綜合考慮上述因素，先通過合理的假設(shè)，首先從易溶急難揮發(fā)性污染物入手，考慮污染物和庫水混合均勻的情況，利用物質(zhì)平衡原理，針對問題中各種情況建立物質(zhì)平衡方程模型。在此基礎(chǔ)上考慮污染物的一維遷移模型。 在建立一維遷移擴展模型時，我們通過轉(zhuǎn)化思想，將污染物對干流形成大面積污染的可能轉(zhuǎn)換為求解在關(guān)閉泄

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

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

19、）附錄二：序號項目單位數(shù)據(jù)相應(yīng)時間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.634512平均出力萬千瓦49.2013最高負(fù)荷萬

20、千瓦109.62003-10-26 14最低負(fù)荷萬千瓦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 function ,it fulfils gen

21、etic 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 coding % you can choose f

22、loat;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=1 ; % length of bit of

23、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% inivitializationfor i=1:sizepop % p

24、roduce 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=individuals.chrom(bestindex,:);avg

25、fitness=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 individuals.chrom=Cross(pc

26、ross,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); en

27、d % 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=indivi

28、duals.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 it is openif ishandle(hfi

29、g) 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(適應(yīng)度曲線 終止代數(shù) num2str(maxgen);xlabel(進化代數(shù));ylabel(適應(yīng)度);legend(平均適應(yīng)度,最佳適應(yīng)度);disp(適應(yīng)度 變量);x=Decode(lenchrom,bound,bestchrom,fcode);% show in command windowvpa(bestfitness,10);vpa(x,

30、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 chromosome% lenchrom input

31、 : 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 coding pick=rand(1,sum(len

32、chrom); 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;endfunction ret=Cross(pcross,lench

33、rom,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 i

34、nput : 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)=index(2) pick=rand(1,2); in

35、dex=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(index(1)=chrom(index(1)

36、-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; while pick=0 pick=rand; end

37、 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,bitget(chrom(index(2),j);

38、 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 position of crossover pick=ra

39、nd; 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),pos)=pick*v1+(1-pick)*v2;

40、 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 % binary coding for i=length(le

41、nchrom):-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(lenchrom):-1:1 data(i)=b

42、itand(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 fulfils a mutation among ch

43、romosomes% 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

44、: 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 continue; end % mutation position pick=rand; while pick=0 pick=rand; end pos=ceil(pick*sum(lenchrom); v=bitget(chrom(index),pos); v=v; chrom1=bitset(chrom(index),pos,v); end ret=chrom; case float % multiple position mutation for i=1:sizepop % select child at random pick=r