1、第1 页安徽工程大学 数学建模 课程设计论文水库排污问题的数学模型姓 名：刘智浩 学 号：3110801208 班 级：数学 112 指导老师：周金明 成 绩： 完成日期： 2013 年 7 月 3 日第2 页摘 要 众所周知，水库污染是一个很严重的问题，讨论干流发生大面积污染的可能性，则必须知道在两小时内流入干流的化学物质的质量，问题的关键在于对流入干流的化学物质的质量的计算，通过实际流入干流的化学物质的质量与警戒值的比，则可讨论干流发生大面积污染的可能。所以本文针对水库突发性事故排污问题，根据物质平衡原理，首先通过建立二维水质污染物浓度模型，给出了单个水库对干流造成大面积污染的可能性；然后
2、建立两水库排污模型，分析了在另一水库有连续点源污染物排放及水流相互影响的情况下，两水库对干流造成大面积污染的可能性大小并进一步针对第三种情况的发生，给出在短时间内控制污染的有效措施；且讨论了若污染物具有挥发性，上述各情况造成干流发生大面积污染的可能性大小，为水库事故性排污问题提供了有价值的理论依据。通过建立这两个模型大大可以减少污染对干流的影响，并且能够节省大量的时间。关键词：水库排污；污染物浓度；流量；水流速度第3 页一、问题的提出1、问题的背景某条江流上有 2 条支流，每条支流上都兴建了规模相当的水库。由于正处在雨水多发季节，因此两个水库都以一定规模的流量进行泄洪。某天晚上 10：00，在
3、其中的一个水库中发生了两船相撞的事故，而其中的一条船装载的 p 吨化学物质（这里的化学物质可以是具有挥发性的，也可能是急难挥发的）全部泄漏至水库中。当水上航运事故处置中心接获事故报告，立即要求该水库关闭水库泄洪闸，以免化学物质随洪水流入干流，发生更大规模的污染。水库闸门开始关闭时，已经处在事故发生后的 1 个小时，而水库闸门彻底关闭也需要 1 个小时的时间。 根据当地环境监测的有关规定，干流大面积污染的危险警戒值设为：三小时内 q 吨该化学物质发生泄漏。2、面临的问题(1) 试建立合理的数学模型，讨论由于此次事故的发生，干流发生大面积污染的可能性；(2) 如果在另外的一水库中有一化工厂违规排放
4、废料。废料中同样含有该化学物质。该工厂为躲避环境监测站的监控，均在晚上 9：00-12：00 违规进行周期性排放。在这种情形下，讨论由于此次事故的发生，干流发生大面积污染的可能性；(3) 如果以上两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么上述结果是否会改变？请给出说明，若有改变，则给出修正的模型及结果；(4) 如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模型。二、问题的基本假设（1） 污染物为速溶物质，因此药品从船上流入水中的时间很少，可以忽略不计；（2） 污染物质从水库中一经流出就进入干流；（3） 水库和河流中的
5、水流都是处于推流状态；（4） 两水库事故发生条件相同，即两水库有相同的客观条件；（5） 被污染的水库关闭泄洪闸后不再有水流流入干流；（6） 不考虑生物等因素在水库泄洪过程中的作用，污染物除了流出外不因腐烂沉积等手段从水中消失；（7） 外界因素不对水库的体积变化产生影响，例如：雨水、地表径流、底下径流等；第4 页（8） 参与模型的变量是连续变化的，并且充分光滑；（9） 不考虑从不同的渠道流入与流出水库之间的区别，只考虑携带污染物的水流入水库和水库中的水流出对水库污染程度的影响，因此可以把水库看成是单流入单流出的系统。三、符号说明（1） ：t 时刻水库水的流入速度；)(r（2） ：t 时刻流入水库
6、的污染物的浓度；1p（3） ：t 时刻水库水的流出的速度；)(0r（4） ：t 时刻流出水库的污染物的浓度；（5） ：t 时刻水库中污染物的浓度；)(p（6） ：t 时刻水库水的体积；v（7） ：计算体积元内该污染物的增量；pm（8） ： 为时间；t（9） ： 为从水库中流出的水中的污染物的浓度；C（10） ：为水库的流出速度，即流量；Q（11） ：为泄洪闸处到污染处的距离；X（12） Q：为水库的流量；（13） q：为排入河流的污水的流量；（14） ：为河流中污染物的本底浓度；1C（15） :为水库中的污染物的浓度；2注意：部分符号见论文中说明。四、模型的建立与求解由问题的分析中知道，流入水
7、库的污染物能以很快的速度与水库中的水均匀混合，也就是说水库中的污染状况在任何局部水体都是一样的，污染程度与水体在水库中的位置无关，因此我们可以建立下面模型。模型一问题一：此次事故的发生，干流发生大面积污染的可能性；根据物质平衡原理和题目假设可知：第5 页水库 1 中污染物的改变量 = 流入的污染物的量 流出污染物的量于是对于充分小的 ，在时间（t，t+ ）内有：ttrtptrpvtvp )()()()()( 01两边同除 ，并使 0 得：tt（1）)()(01tttt现假设 f(t)=p(t)v(t)得： dttvdttpdttdttff pvpf )()()()()( 即原式可以写为：（2）
8、)()()()( 01)()( trtrtvdtpp在水库 1 中发生撞船事故后，污染物处于非稳定排放即： ，而由于水库闸)(dtp门的关闭也势必会引起水库中水的体积变化，故： 。现不考虑流入水库中的水0)(tdv所含有与泄漏污染物相同物质的情况而带来的影响，即可看作 ，另外由问题)(1tp分析中知道：流出的污染物的浓度应与水库中污染物浓度相同，即 这样对)(0t于问题一我们可以得到求解公式：（3）)()()( 0)()( trtpdttvdtpp进一步我们假设从水库中流出的水的流量初始值（从 t=0 时算起）为 ，在关闭闸0r门的过程中，我们假定流量处于线性变化的趋势。这一假设是基于流量与过
9、流面积为线性关系上作出的，进一步可得：（4）)72036.(3602.)(0 ttrtr从上面可看出 为分数函数，这主要是因为水库闸门关闭是在事故发生一小时)(0t后作出的。现在有了 的变化的表达式，为了能求出 的表达式。我们还要写出r )(tp的表达式。首先我们假设水库的体积的初始值为 （t=0 时） ， 值我们可以通过)(tv 0v0v卫星定位系统及所建立的模型求出（具体卫星定位系统模型见附表） 。而跟 相关的)(t第6 页还有 的值。我们假设 为一定值 ，则 随随时间变化的关系式为：)(1tr)(1tr1r)(tvv00由于 为分段函数可知： 也响应的为分段函数，具体函数表达式为：)(0
10、tr)(t（5） )72036.()3602(.)(10 ttrtrvt把（4）式代入（3）式可以得到：当 秒时：60t（6）0)(0101 )()() rtpdttrtvrtp p当 秒时：7203t（7）)3602()3602()18() )(101 trtttrtrvtrrtp p对（5）式化简有：（8）1)(010 )( rtpdttrtvp通过推导的出：（9）cetptrtvr)ln(1 0100)(有已知条件可知： ，故0/vP100/rvpc经简化后：时 （10）00010 /)()( 101 vvtrvtp rr 36t对（7）简化后得；（11） 20010200101 36)
11、2()36()2()(ln tvtrvdtrtrvdt设 ，(r 1-2r0)=b， 得：cv0 a360当 ab42第7 页)ln(21)360()42arctn(42)(ln 2ctbtarbtbrtp （12）)rt(22ctca设： )4rtn(422btbcz最后得到：（13）1)7205.0()ln(5.012)( cetp zabcbtat 保持连续性，当 t=3600 时， = ，此时可得到相应的 值，但由于不确361p C定因素很多，故确定 不是很容易，这主要是缺少数据造成的。1C当 时得到（14）式acb421)720()ln(211 )( Cetpzabcta（14）ca
12、batcbz 42l4 22同样为保持连续性，要求当 t=3600 时， = 。)360(1p)(最后： )(72036.4 )(0)()( 21)720()ln(215.l5.0012 11 stacbce tvtrvtpzabctat rr当当（15）那么 时间内流出水库的污染物的量便可表示为：t（16）trtpQ)(0在 （s）时流出的污染物的量为：360t（17）tdrvptrvrr 360 00011 /)(11在 （s）流出的量为：72t第8 页（18）720360)36()dtrtpQ流出的总量： （19）21再用 Q 与 2/3q 进行比较，便得出是否会发生大面积污染。问题二：
13、如果在另外的一水库中有一化工厂均在晚上 9：00-12：00 违规进行周期性排放同样含有该化学物质的废料。讨论由于此次事故的发生，干流发生大面积污染的可能性。由题目分析中知道，水库 2 中有一化工厂违规排放含有该污染物的废料，而由于两个水库之间没有联系，故我们只需要单独考虑水库 2 的排污量，然后加上水库 1 的污染物排放量，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。下面我们将在水库 1 模型的基础上建立 2 水库的模型。（20）)()()()( 001)()( trtptrtpdttvdttp 设 （常数） ，即在 9：0012：00 这段时间内污染物以一个恒定值流入水库，k
14、t1考虑水库水的流入速度为一定值 ，流出速度也为一个定值1r0r这样水库体积的表达式可写成如下的公式：（21）ttvt010)( )18(t代入上面的表达式可简化为：（22） 01)(010 )()( rtprkdttrp用上式可求出 ，其中 可通过 求得。那20101)()( Cvktpr2kp)(么从 9：0012 ：00 这三个小时内流出闸门的污染物的总量就可以求出来，污染物的量为：（23）1080)(dtrtpQ（24）Q（ 为三小时内从 1 水库和 2 水库流出的污染物总量）Q若 则发生大面积污染；q若 则不会发生大面积污染；问题三：如果两个水库间有一条人工修建的水渠相连接，水渠中的
15、水流流向不定，但保证两水库之间的水流能够相互影响。那么问题二结果是否会改变？由题目分析可知，第9 页当两水库之间有一人工修建的水渠相互连通时，水渠中的水流必然是从高水位流向低水位，为了使模型简化，我们有以下说明：1、由于两水库是连通的，因此在水库 1 关闸前，两水库的液面必然是趋近于等高的，否则必然有水从一个水库流向另一个水库。2、在事故发生后一个小时内，考虑两水库的规模相当，且水库 1 没有关闭泄洪闸，此时两水库彼此不受影响。3、在事故发生一个小时后，考虑到水库 1 要关闭泄洪闸，这势必引起水库 1 的水位上升，由说明 1 可知水库 1 中的水必将通过水渠流向水库 2。从上面的分析可知，在有
16、连通水渠的情况下，我们只需考虑从 1 水库向 2 水库的流入情况，而不必考虑从 2 水库向 1 水库的流入情况。下面我们就该问题给出进一步分析。在第一个小时内：（即 9：0010：00）此时，水库 1 中还未发生事故，只有水库 2 中在排放废料，用问题 2 的模型我们可求出流出的污染物的量为：（25）36001)(dtrtpQ（ 由（20）式确定）)(tp在第二个小时内（即 10：0011：00）水库 1 已发生装船事故，而水库 2 继续排放废料，但泄洪闸还未关闭，所以我们仍然独立考虑。对水库 1 我们用问题一的模型求解得：（26） 360 00012 /)(11 dtrvptrvQrr对水库
17、 2 我们仍有问题二的模型求解得到：（ 由（20）确定） （27）720360)(dtrtp)(t（28）212Q在第三个小时内（即 11：0012：00） ；水库 1 泄洪闸正在关闭，而水库 2 继续排放废料。这时水库 1 和水库 2 就要结合在一起考虑了，如下图所示：第10 页图 一由上面分析可列出以下的方程式：对于水库 1： （29）)()()()( 00)()( trtprtpdtvdtpp 其中 的表达式为：)(tv（ ）trtrtvt )(3602(101 72036t而： （ ）)(00ttr t对水库 2： （30）)()()()()( 001)()( trptrtprtpdt
18、vdtpp 其中 tttv00102由基本假设知道，水库 1 和水库 2 规模相当则： trtrtt 0010)()(假定 为一定值 。)(0tr0r联立解出： ：由于所得表达式非常复杂，我们把 表达式放在附录（三）里表)(tp与 )(tp示。而 由于过于复杂，我们这里就不给出解析解了。)(t第11 页（ ）7203601)(dtrtpQ72036t（31）tt087203)(6)(10trtr故有： （32）3213Q在 12：00 以后一个小时内，1 水库完全关闭但 1 水库的水将通过渠道流入 2 水库内。则有以下式：对水库 1 有： （33）)()()( 0)()( trtpdtvdtp
19、p其中 表达式为： )(tv rtv010对水库 2 有： （34）)()()()( 00)()( trptrtpdttdtptv 其中 表达式为： )(tv ttrt010求得 （35）tpQ14802)(（36）42321Q总若 则干流不会发生大面积污染；q若 则干流会发生大面积污染；Q问题四：如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模型。由问题的分析中知道，我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量 为负值时，从而使得在干流水体污pm染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。模型二为使模型清楚，
20、我们先给出下面所用名词解释。源和漏：是对体积元内污染物变化的一种描述，源是体积元内污染物的增加速率，第12 页漏是体积元内污染物的减少速率。这里的“漏”不意味着漏掉，而有更广泛的意义，如污染物的降解、沉淀和挥发等都属于“漏” 。由于污染物在干流中迁移符合迁移方程，由基本假设中我们知道河流中的水流处于推流状态，也就是体积元中水分子以同一速度向下游运动如图二和图三。设 C(x)和C(x+ )分别为进入水片的水中某中污染物的浓度， 是计算体积元中所含的该污染x pm物的质量。按照质量守衡原理，计算体积元里污染物质量的增量为：txbStxStxQxCxQm slp )()()(（1） （2） （3）
21、（4）+ （37）txASv（5）上式中（1）至（5）各项的意义如下（1）储存量项；（38）)(CAxmp（2）平流输送项；（3）侧向的源和漏，其中 是单位时间内单位长度上的源和漏，它通常是由侧向lS分布流量 带入的，因此可把 写为：ql ilCq这里 是该分布流量中污染物的浓度；iC（4）表面的源和漏， 是单位时间内单位面积上的源和漏（ ） ；sS smg2/（5）体积元内的源和漏， 是单位时间内单位体积内的源和漏（ ） ；v 3图 二第13 页图 三把（38）式带入（37）式，用 除方程两边并令 和 ，则得：tx0xt（39）vsl SAbSQCtA)()(式（10）就是推流时的污染物迁移
22、方程，以后将广泛地利用这个方程，在求解问题四时，我们只要 、 和 的值为负数就可以了，即可以通过关闭人工水渠和部分关闭lSsv水库 2 泄洪闸的手段实现短时间内控制污染。五、模型的评价模型的优点：本模型在建立模型前，由于水体运动非常复杂，污染物的性质也多种多样，污染物在水中的迁移情况受各种因素限制，不容易全部考虑。本文综合考虑上述因素，先通过合理的假设，首先从易溶急难挥发性污染物入手，考虑污染物和库水混合均匀的情况，利用物质平衡原理，针对问题中各种情况建立物质平衡方程模型。在此基础上考虑污染物的一维迁移模型。在建立一维迁移扩展模型时，我们通过转化思想，将污染物对干流形成大面积污染的可能转换为求
23、解在关闭泄洪闸的过程中，水库中多大面积水域内发生事故时才会对干流形成大面积污染，使得问题求解方案明了，为了验证模型的正确性，我们通过遗传算法对事例进行求解，很好的验证了模型。模型的缺点：模型在建立过程中我们所考虑的此污染物为易溶、不挥发的物质，而没有考虑污染物为油性、挥发性和沉淀物质。另外，我们在考虑水体流动时没有考虑水体上、中、下层的不同流动情况。在问题四的求解时，我们只考虑推流问题，而实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题。时间有限，未能对其他方面做深入的探讨，使得模型十分完美。参考书目1 W.金士博 联邦德国 水环境数学模型 中国建筑工业出版社 第14
24、 页1987.102 岳天祥 资源环境数学模型手册 科学出版社 2003.103 余常昭 环境流体力学导论 清华大学出版社 1992.104 刘振航 数学建模 中国人民大学出版社 2004.055 陈春云 郑彤 环境系统数学模型 化学工业出版社 2003.046 周建华 黄燕 MATLAB5.3 北京大学出版社 2000.127 丁春利 精通 MATLAB6 清华大学出版社 2002.06附录一：水库实际库容及水域面积计算公式： )(ihHsVnS式中，s 为 DEM 象元面积；H 为水库水位高程； 为库域 DEM 各像元高程值；n 为库i域 DEN 像元数。（曾永年，马海州，沙占江等；龙羊峡
25、库区环境动态监测信息系统的建立与应用。遥感学报，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:0
26、0 10 最低库水位 米 138.3 2004-7-16 8:00 11 大厂发电量 亿千瓦时 44.6345 12 平均出力 万千瓦 49.20 13 最高负荷 万千瓦 109.6 2003-10-26 14 最低负荷 万千瓦 11.2 2004-7-25 第15 页15 最大调峰 万千瓦 69.8 2003-10-216 最小调峰 万千瓦 42.4 2004-9-12 17 平均调峰 万千瓦 32.68 18 调峰电量 亿千瓦时 .0912 19 弃水量 亿立方米 119.946 20 平均耗水率 立方米/千瓦时 7.79 21 2x2 发电量 亿千瓦时 3.4254 22 大、小厂总电量
27、 亿千瓦时 48.0599附录三： )1%306301arctn(:2% 9036089: )1/2530 1%)180ln(4exp()( 2 2201020 202011 rrt rrrvr trtrtrvCt 附录四：遗传算法程序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 generationsi
28、zepop=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 float;grey;binary pcross=0.6; % probablity of crossover,between 0 and 1fcross=float; %
29、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;第16 页lenchrom=1 ; % length of bit of every variblebound=0 100000;%-individuals=struct(fitness,zeros(1,sizepop),.%value,z
30、eros(1,sizepop),.chrom,); % structure of populationavgfitness=; % average fitness of population bestfitness=; % best fitness of populationbestchrom=; % chromosome of best fitness% inivitializationfor i=1:sizepop% produce new population at randomindividuals.chrom(i,:)=Code(lenchrom,fcode,bound);x=Dec
31、ode(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,:);avgfitness=sum(individuals.fitness)/sizepop;% record average and best fitness of every gener
32、ationtrace=avgfitness bestfitness; % evolution beginfor i=1:maxgen% selectionindividuals=Select(individuals,sizepop,fselect);avgfitness=sum(individuals.fitness)/sizepop;% crossoverindividuals.chrom=Cross(pcross,lenchrom,individuals.chrom,.sizepop,fcross,i maxgen);% mutationindividuals.chrom=Mutation
33、(pmutation,lenchrom,individuals.chrom,.sizepop,fmutation,i maxgen,bound);% calculate fitness for j=1:sizepopx=Decode(lenchrom,bound,individuals.chrom(j,:),fcode);individuals.fitness(j)=Aimfunc(x); end% substitute chromosome of worest fitness% find minimum value which is best第17 页newbestfitness,newbe
34、stindex=min(individuals.fitness);worestfitness,worestindex=max(individuals.fitness);% substitute chromosome of worest fitnessif bestfitnessnewbestfitnessbestfitness=newbestfitness;bestchrom=individuals.chrom(newbestindex,:);endindividuals.chrom(worestindex,:)=bestchrom;individuals.fitness(worestinde
35、x)=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);elsehfig=figure(Tag,trace);endfigure(hfig);r c=size(trace);plot(1:r,trace(:,1),r-,1:r,trace(:,2
36、),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);第18 页% 求最小值z=quadl(f,0.0001,7200,x);ret=abs(z-
37、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% b
38、ound input : boundary of varibles% ret output: chromosomeswitch optscase binary % binary codingpick=rand(1,sum(lenchrom);bits=ceil(pick-0.5);temp=sum(lenchrom)-1:-1:0;ret=sum(bits.*(2.temp);case grey % grey codingpick=rand(1,sum(lenchrom);bits=ceil(pick-0.5);greybits=bits;for i=2:length(greybits)gre
39、ybits(i)=bitxor(bits(i-1),bits(i);endtemp=sum(lenchrom)-1:-1:0;ret=sum(greybits.*(2.temp);case float % float codingpick=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 chr
40、omosomes% pcorss input : probability of crossover% lenchrom input : length of a chromosome第19 页% 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 ou
41、tput : new set of chromosomeswitch optscase simple % cross at single positionfor i=1:sizepop% select two children at randompick=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 crossoverpick=rand;if pickpcrossconti
42、nue;end% random position of crossoverpick=rand;while pick=0pick=rand;endpos=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;endret=chrom;case uniform % uniform cross
43、for i=1:sizepop% select two children at randompick=rand(1,2);while prod(pick)=0pick=rand(1,2);endindex=ceil(pick.*sizepop);% random position of crossoverpick=rand;while pick=0pick=rand;第20 页endif pickpcrosscontinue;end% random position of crossoverpick=rand;while pick=0pick=rand;endmask=2ceil(pick*s
44、um(lenchrom);chrom1=chrom(index(1);chrom2=chrom(index(2);for j=1:sum(lenchrom)v=bitget(mask,j); % from lower to higher bitif v=1chrom1=bitset(chrom1,.j,bitget(chrom(index(2),j);chrom2=bitset(chrom2,.j,bitget(chrom(index(1),j);endendchrom(index(1)=chrom1;chrom(index(2)=chrom2;endret=chrom;case floatf
45、or i=1:sizepop% select two children at randompick=rand(1,2);while prod(pick)=0pick=rand(1,2);endindex=ceil(pick.*sizepop);% random position of crossoverpick=rand;while pick=0pick=rand;endif pickpcrosscontinue;end% random position of crossoverpick=rand;while pick=0第21 页pick=rand;endpos=ceil(pick.*sum
46、(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;endret=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 var
