B谭金梅B数学112班水库排污问题正文

1、 水 库 排 污 问 题 姓 名： 谭 金 梅 班 级： 数 学 112 学 号： 3110801242 指导老师： 周 金 明成 绩： 完成日期： 2013年7月3日 摘 要数学建模是一种数学的思考方法，是运用数学的语言和方法，通过抽象、简化建立能近似刻画并解决实际问题的一种强有力的数学手段。数学建模就是用数学语言描述实际现象的过程。本文针对水库排污问题，首先通过建立二维水质污染物浓度模型，给出了单个水库对干流造成大面积污染的可能性；然后建立两水库排污模型，分析了在另一水库有连续点源污染物排放及水流相互影响的情况下，两水库对干流造成大面积污染的可能性；并进一步针对第三种情况的发生，给出在短时

2、间内控制污染的有效措施；且讨论了若污染物具有挥发性，上述各情况造成干流发生大面积污染的可能性大小，为水库事故性排污问题提供了有价值的理论依据。问题一，只要考虑在事故发生到关闭水库的两个小时内，流出水库的污染物的质量小于吨即可。问题二中，由于两个水库之间没有联系，只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响。问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都是从高水位流向低水位，因此，只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。问题四，通过稀释原理，建立污染物在干流中迁移迁移模

3、型，使干流水体计算体积元内该污染物的增量为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。关键字：水流速度、污染物浓度、流量、水库排污、MATLAB问题重述某条江流上有2条支流，每条支流上都兴建了规模相当的水库。由于正处在雨水多发季节，因此两个水库都以一定规模的流量进行泄洪。某天晚上10：00，在其中的一个水库中发生了两船相撞的事故，而其中的一条船装载的p吨化学物质（这里的化学物质可以是具有挥发性的，也可能是急难挥发的）全部泄漏至水库中。当水上航运事故处置中心接获事故报告，立即要求该水库关闭水库泄洪闸，以免化学物质随洪水流入干流，发生更大规模的污染。水

4、库闸门开始关闭时，已经处在事故发生后的1个小时，而水库闸门彻底关闭也需要1个小时的时间。 根据当地环境监测的有关规定，干流大面积污染的危险警戒值设为：三小时内q吨该化学物质发生泄漏。(1) 试建立合理的数学模型，讨论由于此次事故的发生，干流发生大面积污染的可能性；(2) 如果在另外的一水库中有一化工厂违规排放废料。废料中同样含有该化学物质。该工厂为躲避环境监测站的监控，均在晚上9：00-12：00违规进行周期性排放。在这种情形下，讨论由于此次事故的发生，干流发生大面积污染的可能性；(3) 如果以上两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那

5、么上述结果是否会改变？请给出说明，若有改变，则给出修正的模型及结果； (4) 如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模型。二、问题分析由题目知道，事故发生时两水库都正在泄洪，因此此时水库中的水流速水库较快。而泄露到水库中的化学物质不论是具有挥发性的，还是急难挥发的，它们对干流污染的情况总是类同的，因此我们总可以认为污染物是易溶急难挥发性物质。为使我们的模型简单，我们可以先假设事故发生在水库1中，污染物在水库中的分布是符合零维迁移模型的，此时流入水库的污染物能以很快的速度与水库中的水均匀混合，水库中任何部份水体的污染状况都是一样的，污染程度与水体在水库中的位置无关。而

6、实际上，污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，但在这样一个突发事件要求短时间内得到控制的问题中，我们总可以用污染物瞬时混合均匀状态模型来代替污染物一维迁移的过程。在解答问题一时，我们只要考虑在事故发生到关闭水库的两个小时内，流出水库的污染物的质量小于吨即可。问题二中提到的情况只是将干流的污染源从一个增加到两个，那么发生大面积污染的可能性就要增大。由于两个水库之间没有联系，我们只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都

7、是从高水位流向低水位，因此，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。 分析问题四，干流已经发生大面积污染，对比水体污染物处理的各种手段，不论是化学、物理还是生物手段都不可能在短时间内除去污染物，因此我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。三、模型假设污染物为速溶物质，因此药品从船上流入水中的时间很少，可以忽略不计；污染物质从水库中一经流出就进入干流；水库和河流中的水流都是处于推流状态；两水库事故发生条件相同，即两水库

8、有相同的客观条件；被污染的水库关闭泄洪闸后不再有水流流入干流；不考虑生物等因素在水库泄洪过程中的作用，污染物除了流出外不因腐烂沉积等手段从水中消失；外界因素不对水库的体积变化产生影响，例如：雨水、地表径流、底下径流等；参与模型的变量是连续变化的，并且充分光滑；不考虑从不同的渠道流入与流出水库之间的区别，只考虑携带污染物的水流入水库和水库中的水流出对水库污染程度的影响，因此可以把水库看成是单流入单流出的系统。符号说明 （1）：t时刻水库水的流入速度；（2）：t时刻流入水库的污染物的浓度；（3）：t时刻水库水的流出的速度；（4）：t时刻流出水库的污染物的浓度；（5）：t时刻水库中污染物的浓度；（6

9、）：t时刻水库水的体积；（7）：计算体积元内该污染物的增量；（8） ： 为时间；（9） ： 为从水库中流出的水中的污染物的浓度；（10） ：为水库的流出速度，即流量；（11）：为泄洪闸处到污染处的距离；（12） Q：为水库的流量；（13） q：为排入河流的污水的流量；（14）：为河流中污染物的本底浓度；（15）:为水库中的污染物的浓度； 注意：部分符号见论文中说明。五、模型的建立与求解5.1 模型一由问题的分析中知道，流入水库的污染物能以很快的速度与水库中的水均匀混合，也就是说水库中的污染状况在任何局部水体都是一样的，污染程度与水体在水库中的位置无关，因此我们可以建立下面模型。问题一、此次事故

10、的发生，干流发生大面积污染的可能性；根据物质平衡原理和题目假设可知：水库1中污染物的改变量 = 流入的污染物的量 流出污染物的量于是对于充分小的，在时间（t，t+）内有：两边同除，并使0得： （1）现假设f(t)=p(t)v(t)得：即原式可以写为： （2）在水库1中发生撞船事故后，污染物处于非稳定排放即：，而由于水库闸门的关闭也势必会引起水库中水的体积变化，故：。现不考虑流入水库中的水所含有与泄漏污染物相同物质的情况而带来的影响，即可看作，另外由问题分析中知道：流出的污染物的浓度应与水库中污染物浓度相同，即这样对于问题一我们可以得到求解公式： （3）进一步我们假设从水库中流出的水的流量初始值

11、（从t=0时算起）为，在关闭闸门的过程中，我们假定流量处于线性变化的趋势。这一假设是基于流量与过流面积为线性关系上作出的，进一步可得： （4）从上面可看出为分数函数，这主要是因为水库闸门关闭是在事故发生一小时后作出的。现在有了的变化的表达式，为了能求出的表达式。我们还要写出的表达式。首先我们假设水库的体积的初始值为（t=0时），值我们可以通过卫星定位系统及所建立的模型求出（具体卫星定位系统模型见附表）。而跟相关的还有的值。我们假设为一定值，则随随时间变化的关系式为： 由于为分段函数可知：也响应的为分段函数，具体函数表达式为： （5）把（4）式代入（3）式可以得到：当秒时： （6）当秒时： （7

12、）对（5）式化简有： （8）通过推导的出： （9）有已知条件可知：，故经简化后： 时 （10）对（7）简化后得；（11）设，(r1-2r0)=b，得：当 （12）设：最后得到： （13）保持连续性，当t=3600时，=，此时可得到相应的值，但由于不确定因素很多，故确定不是很容易，这主要是缺少数据造成的。当时得到（14）式 （14）同样为保持连续性，要求当t=3600时，=。最后： （15）那么时间内流出水库的污染物的量便可表示为： （16）在（s）时流出的污染物的量为： （17）在（s）流出的量为： （18）流出的总量： （19）再用Q与2/3q进行比较，便得出是否会发生大面积污染。问题二：如

13、果在另外的一水库中有一化工厂均在晚上9：00-12：00违规进行周期性排放同样含有该化学物质的废料。讨论由于此次事故的发生，干流发生大面积污染的可能性。由题目分析中知道，水库2中有一化工厂违规排放含有该污染物的废料，而由于两个水库之间没有联系，故我们只需要单独考虑水库2的排污量，然后加上水库1的污染物排放量，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。下面我们将在水库1模型的基础上建立2水库的模型。 （20）设（常数），即在9：0012：00这段时间内污染物以一个恒定值流入水库，考虑水库水的流入速度为一定值，流出速度也为一个定值这样水库体积的表达式可写成如下的公式： （21）代入

14、上面的表达式可简化为： （22） 用上式可求出，其中可通过求得。那么从9：0012：00这三个小时内流出闸门的污染物的总量就可以求出来，污染物的量为： （23） （24）（为三小时内从1水库和2水库流出的污染物总量）若则发生大面积污染；若则不会发生大面积污染；问题三：如果两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么问题二结果是否会改变？由题目分析可知，当两水库之间有一人工修建的水渠相互连通时，水渠中的水流必然是从高水位流向低水位，为了使模型简化，我们有以下说明：1、由于两水库是连通的，因此在水库1关闸前，两水库的液面必然是趋近于等高的，否

15、则必然有水从一个水库流向另一个水库。2、在事故发生后一个小时内，考虑两水库的规模相当，且水库1没有关闭泄洪闸，此时两水库彼此不受影响。3、在事故发生一个小时后，考虑到水库1要关闭泄洪闸，这势必引起水库1的水位上升，由说明1可知水库1中的水必将通过水渠流向水库2。从上面的分析可知，在有连通水渠的情况下，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。下面我们就该问题给出进一步分析。在第一个小时内：（即9：0010：00）此时，水库1中还未发生事故，只有水库2中在排放废料，用问题2的模型我们可求出流出的污染物的量为： （25）（由（20）式确定）在第二个小时内（即10

16、：0011：00）水库1已发生装船事故，而水库2继续排放废料，但泄洪闸还未关闭，所以我们仍然独立考虑。对水库1我们用问题一的模型求解得： （26）对水库2我们仍有问题二的模型求解得到： （由（20）确定） （27） （28）在第三个小时内（即11：0012：00）；水库1泄洪闸正在关闭，而水库2继续排放废料。这时水库1和水库2就要结合在一起考虑了，如下图所示：图 一由上面分析可列出以下的方程式：对于水库1： （29）其中的表达式为： （）而： （）对水库2： （30）其中 由基本假设知道，水库1和水库2规模相当则：假定为一定值。联立解出：由于所得表达式非常复杂，我们把表达式放在附录（三）里表示

17、。而由于过于复杂，我们这里就不给出解析解了。 （） （31）故有： （32）在12：00以后一个小时内，1水库完全关闭但1水库的水将通过渠道流入2水库内。则有以下式：对水库1有： （33）其中表达式为： 对水库2有： （34）其中表达式为： 求得 （35） （36）若则干流不会发生大面积污染；若则干流会发生大面积污染；问题四：如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模型。5.2 模型二由问题的分析中知道，我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间

18、内达到控制污染。为使模型清楚，我们先给出下面所用名词解释。源和漏：是对体积元内污染物变化的一种描述，源是体积元内污染物的增加速率，漏是体积元内污染物的减少速率。这里的“漏”不意味着漏掉，而有更广泛的意义，如污染物的降解、沉淀和挥发等都属于“漏”。 由于污染物在干流中迁移符合迁移方程，由基本假设中我们知道河流中的水流处于推流状态，也就是体积元中水分子以同一速度向下游运动如图二和图三。设C(x)和C(x+)分别为进入水片的水中某中污染物的浓度，是计算体积元中所含的该污染物的质量。按照质量守衡原理，计算体积元里污染物质量的增量为： （1） （2） （3） （4） + （37） （5）上式中（1）至（

19、5）各项的意义如下（1）储存量项； （38）（2）平流输送项；（3）侧向的源和漏，其中是单位时间内单位长度上的源和漏，它通常是由侧向分布流量带入的，因此可把写为：这里是该分布流量中污染物的浓度；（4）表面的源和漏，是单位时间内单位面积上的源和漏（）；（5）体积元内的源和漏，是单位时间内单位体积内的源和漏（）；图 二图 三把（38）式带入（37）式，用除方程两边并令和，则得： （39）式（10）就是推流时的污染物迁移方程，以后将广泛地利用这个方程，在求解问题四时，我们只要、和的值为负数就可以了，即可以通过关闭人工水渠和部分关闭水库2泄洪闸的手段实现短时间内控制污染。六、模型的扩展 从问题的分析中

20、可以知道，我们知道污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，为使模型比较符合实际，我们对上面模型进行扩展，建立污染物一维迁移扩展模型。扩展模型一从水库中流出的污染物的量可用下面的式子表示：由于在水库中污染物浓度在各个位置是不同的，而在整个水库中各个量又具有一定的连续性，故考虑时刻在闸门口处时间内通过的量为：，此式中我们考虑与，的关系，其中方向为水库流水的流向。这是因为对于事故泄漏排放的污染物，人们最关心的是其通过某一位置的时间，最大浓度等数值。因此，研究一维流场中的瞬时排放污染物的分布特征就非常必要了，我们用下式为基础来考虑分析泄漏时污染物在闸门口的浓度变化特征。在瞬时

21、排放条件下：，此时有下式： （40）考虑在泄洪闸处：为泄洪闸处到污染处的距离，且（守恒污染物） （41）其中： 考虑在一个极其短暂的时间内，在这个内污染物的浓度可近似看作是不变的。则在时刻有： 通过对上式积分可求出2个小时内流出水库的污染物的量为： （42）这里P，A，Q，均为已知，只是发生的地点未知。通过上式我们可以求出在两个小时内通过的量为吨时的，但由于式（42）的解析解很难求出，我们就把它写成下面的形式： （43）把（2）式的相应的量代入后，通过遗传算法可求出的最小值，而（43）式的最小值为零，由此可求出满足时所对应的，求出后，就可以得到发生大面积污染的可能的点所在的范围，如下图所示：A

22、2A1 图 四A1为发生大面积污染的点的集合。A2为不发生大面积污染的点的集合。因此对于问题一：发生大面积污染的可能为： ； （44）问题二：由题目分析知道，若水库2中有一工厂在9：0012：00违规进行排放，可知这样会增大发生大面积污染的可能性。因为另外一水库处于稳定流动状态，工厂排放的污染物为稳定的连续排放，环境中污染物的分布状况也是稳定的，综合考虑可有下式： （45）其中： ，。我们假定工厂距泄洪闸的距离可以确定为则： 污染物到达泄洪闸的时间为：则泄漏的量为： （46）由此我们可以得到在另外一水库当中有工厂进行周期性排放下流出的总量M为： （47）即我们把上式转化为的形式，这里我们可以看

23、出的表达式与问题（1）中有所不同，这是因为我们考虑的时间已从原来的2小时变成3小时了，对于问题（二）我们仍然用问题（一）的解法，即通过遗传算法(程序见附录四)求出发生大面积污染时的。通过这个我们可求出发生大面积污染的点的集合，进而求出发生大面积污染的可能性为，A是水库的面积,可知是个定值。 在这里我们可以给出扩展模型一求解下的事例； 其中：m/s =1513 部分数据见附录二。 通过MATLAB求解得到结果如下： 问题一发生大面积污染的点的集合类圆半径； 问题二发生大面积污染的点的集合类圆半径。 图形见图五图六：图五 问题一发生大面积污染的可能半径图六 问题二发生大面积污染的可能半径 扩展二从

24、问题的分析及假设可以知道，我们的模型在建立过程中我们所考虑的只是一种污染物，即此污染物为易溶、不挥发的物质。实际上污染物可能为油性、挥发性和沉淀物质。对于油性物质我们应该考虑该物质的密度是否大于水的密度，若油性污染物小于水的密度，因为油性物质不溶于水，该污染物应该是漂浮在水面上的，而且呈油膜状，因此我们就应考虑水库表层水的流动情况；若油性污染物的密度大于水的密度，则污染物将会沉如水底，那么我们就应该考虑深沉水体的水流问题；若油性污染物等于水的密度，那么污染物就会悬浮在水库水体中，因此我们就应该考虑整个水体上、中、下层的流动情况（因为不同水层的水流的流动情况是不同的）。对于挥发性的污染物，我们还

25、需要考虑此污染物在水和空气中的分配比例。在考虑问题四的时候我们考虑的是推流问题，实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题而。为使模型更符合实际，模型必须在上述方面在做深入探讨。七、模型的评价模型的优点： 本模型在建立模型前，由于水体运动非常复杂，污染物的性质也多种多样，污染物在水中的迁移情况受各种因素限制，不容易全部考虑。本文综合考虑上述因素，先通过合理的假设，首先从易溶急难挥发性污染物入手，考虑污染物和库水混合均匀的情况，利用物质平衡原理，针对问题中各种情况建立物质平衡方程模型。在此基础上考虑污染物的一维迁移模型。 在建立一维迁移扩展模型时，我们通过转化思想

26、，将污染物对干流形成大面积污染的可能转换为求解在关闭泄洪闸的过程中，水库中多大面积水域内发生事故时才会对干流形成大面积污染，使得问题求解方案明了，为了验证模型的正确性，我们通过遗传算法对事例进行求解，很好的验证了模型。模型的缺点：模型在建立过程中我们所考虑的此污染物为易溶、不挥发的物质，而没有考虑污染物为油性、挥发性和沉淀物质。另外，我们在考虑水体流动时没有考虑水体上、中、下层的不同流动情况。在问题四的求解时，我们只考虑推流问题，而实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题。时间有限，未能对其他方面做深入的探讨，使得模型十分完美。参考书目1 W.金士博联邦德国

27、水环境数学模型 中国建筑工业出版社 1987.102 岳天祥 资源环境数学模型手册 科学出版社 2003.103 余常昭 环境流体力学导论 清华大学出版社 1992.104 刘振航 数学建模 中国人民大学出版社 2004.055 陈春云 郑彤 环境系统数学模型 化学工业出版社 2003.046 周建华 黄燕 MATLAB5.3 北京大学出版社 2000.127 丁春利 精通MATLAB6 清华大学出版社 2002.06附录一：水库实际库容及水域面积计算公式：式中，s为DEM象元面积；H为水库水位高程；为库域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.6

29、34512平均出力万千瓦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附录三：附录四：遗传算法程序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(适应度曲线 终止代数 num2str(maxgen);xlabel(进化代数);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