B水库排污问题

安徽工程大学 数学建模 课程设计论文

水库排污问题

班 级： 数学112 姓 名：学 号： 3110801204 指导老师： 周金明 成 绩： 完成日期： 2013年7月3日

摘 要

本文针对水库突发性事故排污问题，首先通过问题一建立二维水质污染物浓度模

2

型，只要考虑在事故发生到关闭水库的两个小时内，流出水库的污染物的质量小于q

3吨，给出了单个水库对干流造成大面积污染的可能性 Q

W1W2

；然后根据问题二、

2

q10003

三建立两水库排污模型，问题二由于两个水库之间没有联系，我们只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了，问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都是从高水位流向低水位，因此，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。分析了在另一水库有连续点源污染物排放及水流相互影响的情况下，两水库对干流造成大面积污染的可能性大小

Q

W1'W2'W3'

q1000

。

进一步针对第三种情况的发生，给出在短时间内控制污染的有效措施；且讨论了若污染物具有挥发性，上述各情况造成干流发生大面积污染的可能性大小，为水库事故性排污问题提供了有价值的理论依据。环境容量是环境科学的一个基本理论问题，也是环境管理中的一个重要的实际应用问题。在实践中，环境容量是环境目标管理的基本依据，是环境规划的主要约束条件，也是污染物总量控制的关键技术支持。从环境管理、监测与监督的角度出发，水环境容量是指水体在设计水文条件和规定的环境目标下所能容纳的最大污染物量。分析库区的污染状况，提出了水库环境容量的计算原则、设计水文条件和水质保护目标。根据岸边环境容量计算结果，提出了三峡水库污染混合区的控制标准。同时，推荐了三峡水库水环境容量的综合方案。

关键词：水库排污；污染物浓度；流量；水流速度；环境容量

某条江流上有2条支流，每条支流上都兴建了规模相当的水库。由于正处在雨水多发季节，因此两个水库都以一定规模的流量进行泄洪。某天晚上10：00，在其中的一个水库中发生了两船相撞的事故，而其中的一条船装载的p吨化学物质（这里的化学物质可以是具有挥发性的，也可能是急难挥发的）全部泄漏至水库中。当水上航运事故处置中心接获事故报告，立即要求该水库关闭水库泄洪闸，以免化学物质随洪水流入干流，发生更大规模的污染。水库闸门开始关闭时，已经处在事故发生后的1个小时，而水库闸门彻底关闭也需要1个小时的时间。 根据当地环境监测的有关规定，干流大面积污染的危险警戒值设为：三小时内q吨该化学物质发生泄漏。

阅读完以上资料后，结合查询资料，请做出合理假设并完成以下要求：

(1) 试建立合理的数学模型，讨论由于此次事故的发生，干流发生大面积污染的可能性；

(2) 如果在另外的一水库中有一化工厂违规排放废料。废料中同样含有该化学物质。该工厂为躲避环境监测站的监控，均在晚上9：00-12：00违规进行周期性排放。在这种情形下，讨论由于此次事故的发生，干流发生大面积污染的可能性；

(3) 如果以上两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么上述结果是否会改变？请给出说明，若有改变，则给出修正的模型及结果；

(4) 如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模型。

由题目的背景知道，事故发生时两水库都正在泄洪，因此此时水库中的水流速水库较快。而泄露到水库中的化学物质不论是具有挥发性的，还是急难挥发的，它们对干流污染的情况总是类同的，因此我们总可以认为污染物是易溶急难挥发性物质。为使我们的模型简单，我们可以先假设事故发生在水库1中，污染物在水库中的分布是符合零维迁移模型的，此时流入水库的污染物能以很快的速度与水库中的水均匀混合，水库中任何部份水体的污染状况都是一样的，污染程度与水体在水库中的位置无关。而实际上，污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，但在这样一个突发事件要求短时间内得到控制的问题中，我们总可以用污染物瞬时混合均匀状态模型来代替污染物一维迁移的过程。

在解答问题一时，我们只要考虑在事故发生到关闭水库的两个小时内，流出水库的

2

污染物的质量小于q吨即可。

3

问题二中提到的情况只是将干流的污染源从一个增加到两个，那么发生大面积污染的可能性就要增大。由于两个水库之间没有联系，我们只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。

问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都是从高水位流向低水位，因此，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。

分析问题四，干流已经发生大面积污染，对比水体污染物处理的各种手段，不论是化学、物理还是生物手段都不可能在短时间内除去污染物，因此我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量mp为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。

则我们可假设：

（1） 污染物为速溶物质，因此药品从船上流入水中的时间很少，可以忽略不计； （2） 污染物质从水库中一经流出就进入干流； （3） 水库和河流中的水流都是处于推流状态；

（4） 两水库事故发生条件相同，即两水库有相同的客观条件； （5） 被污染的水库关闭泄洪闸后不再有水流流入干流； （6） 不考虑生物等因素在水库泄洪过程中的作用，污染物除了流出外不因腐烂沉积等

手段从水中消失；

（7） 外界因素不对水库的体积变化产生影响，例如：雨水、地表径流、底下径流等； （8） 参与模型的变量是连续变化的，并且充分光滑； （9） 不考虑从不同的渠道流入与流出水库之间的区别，只考虑携带污染物的水流入水

库和水库中的水流出对水库污染程度的影响，因此可以把水库看成是单流入单流出的系统。

三、符号的约定

注意：部分符号见论文中说明。

四、模型的建立与求解

由问题的分析中知道，流入水库的污染物能以很快的速度与水库中的水均匀混合，也就是说水库中的污染状况在任何局部水体都是一样的，污染程度与水体在水库中的位置无关，因此我们可以建立下面模型。

4.1 模型一

问题一：此次事故的发生，干流发生大面积污染的可能性； 根据物质平衡原理和题目假设可知：

水库1中污染物的改变量 = 流入的污染物的量 － 流出污染物的量

于是对于充分小的t，在时间（t，t+t）内有：

p(tt)v(tt)p(t)v(t)[p1(t)r1(t)p0(t)r0(t)]t

两边同除t，并使t－>0得：

p(tt)v(tt)p(t)v(t)

p1(t)r1(t)p0(t)r0(t) （1）

t

现假设f(t)=p(t)×v(t)得：

df(t)d(p(t)v(t))dv(t)dp(t)f(tt)f(t)p(t)v(t)

tdtdtdtdt

即原式可以写为：

p(t)

dv(t)dt

v(t)

dp(t)dt

p1(t)r1(t)p0(t)r0(t) （2）

dp(t)dt

在水库1中发生撞船事故后，污染物处于非稳定排放即：

dv(t)dt

0，而由于水库闸

门的关闭也势必会引起水库中水的体积变化，故：0。现不考虑流入水库中的水

所含有与泄漏污染物相同物质的情况而带来的影响，即可看作p1(t)0，另外由问题分析中知道：流出的污染物的浓度应与水库中污染物浓度相同，即p(t)p0(t)这样对于问题一我们可以得到求解公式：

p(t)

dv(t)dt

v(t)

dp(t)dt

p(t)r0(t) （3）

进一步我们假设从水库中流出的水的流量初始值（从t=0时算起）为r0，在关闭闸

门的过程中，我们假定流量处于线性变化的趋势。这一假设是基于流量与过流面积为线

性关系上作出的，进一步可得：

r0...................(0t3600)

r0t r0(t) （4）

2r...........(3600t720)003600

从上面可看出r0(t)为分数函数，这主要是因为水库闸门关闭是在事故发生一小时后

作出的。现在有了r0(t)的变化的表达式，为了能求出p(t)的表达式。我们还要写出v(t)的表达式。首先我们假设水库的体积的初始值为v0（t=0时），v0值我们可以通过卫星定位系统及所建立的模型求出（具体卫星定位系统模型见附表）。而跟v(t)相关的还有r1(t)的值。我们假设r1(t)为一定值r1，则v(t)随随时间变化的关系式为：

v(t)v0r1tr0(t)t

由于r0(t)为分段函数可知：v(t)也响应的为分段函数，具体函数表达式为：

.........(0t3600)v0r1tr0t..........

r0t v(t) （5）

vrt(2r)t...........(3600t7200)0103600

把（4）式代入（3）式可以得到：

当0t3600秒时：

p(t)(r1r0)(v0r1tr0t)当3600t7200秒时：

p(t)(r12r0

dp(t)r0rtrt

t)[v0r1t(2r00)t]p(t)(2r00) （7） 18003600dt3600

dp(t)dt

p(t)r0 （6）

对（5）式化简有：

(v0r1tr0t)通过推导的出：

dp(t)dt

p(t)r1 （8）

p(t)e

r0[

1

ln(v0r1tr0t)]r1r0

c （9）

r1

r0r1

有已知条件可知：p(0)P/v0，故cp/v0v0经简化后：

r1

r1r0

p(t)[v0(r1r0)t]v0

r1

r1r0

p/v0

0t360时0 （10）

对（7）简化后得；

lnp(t)(r1)

dt

r

v0(r12r0)t0t2

3600

r0

a得： 3600

(

r0

)3600

tdt

v

v0(r12r0)t0t2

3600

（11）

设v0c，(r1-2r0)=b，当b24ac

lnp(t)(r)[

24acb2

2atb4acb2

)](

r01

)[ln(at2btc) 36002a

24acb

b22atb)] （12）

222a4acb4acb2atb4acb

2

设：z

2

)

最后得到：

p1(t)e

0.5lnat(2btc)(0.5b720a)0z

c1 （13）

保持连续性，当t=3600时，p1(3600)=p(3600)，此时可得到相应的C1值，但由于不确定因素很多，故确定C1不是很容易，这主要是缺少数据造成的。

当b24ac时得到（14）式

p1(t)e

1

1b

ln(at2btc)(7200a)z'22

C1

2atbb24ac

c （14） z'ln

2b24ac2atbb4ac

同样为保持连续性，要求当t=3600时，p1(3600)=p(3600)。 最后：

r1r1

v0(r1r0)tr1r0v0r1r0...................................0t3600(s)0.5ln(at2btc)(0.5b7200a)z2

p(t)ec..........当b4ac........3600t7200(s) （15）1

1b2'ln(atbtc)(7200a)z

222ec...............当b4ac.........3600t7200(s)1

那么t时间内流出水库的污染物的量便可表示为：

Qp(t)r0(t)t （16） 在0t3600（s）时流出的污染物的量为：

r1

r1r1r0

Q1[v0(r1r0)t]v0r1r0p/v0r0dt （17）



在3600t7200（s）流出的量为：

3600

Q23600p(t)(2r0

7200

r0t

)dt （18） 3600

流出的总量： QQ1Q2 （19） 再用Q与2/3q进行比较，便得出是否会发生大面积污染。

问题二：如果在另外的一水库中有一化工厂均在晚上9：00-12：00违规进行周期性排放同样含有该化学物质的废料。讨论由于此次事故的发生，干流发生大面积污染的可能性。 由题目分析中知道，水库2中有一化工厂违规排放含有该污染物的废料，而由于两个水库之间没有联系，故我们只需要单独考虑水库2的排污量，然后加上水库1的污染物排放量，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。下面我们将在水库1模型的基础上建立2水库的模型。 p(t)

dv(t)dt

v(t)

dp(t)dt

p1(t)r1(t)p0(t)r0(t) （20）

设p1(t)k（常数），即在9：00~12：00这段时间内污染物以一个恒定值流入水库，考虑水库水的流入速度为一定值r1，流出速度也为一个定值r0

这样水库体积的表达式可写成如下的公式：

v(t)v0r1tr0t

dp(t)dt

(0t1080) 0 （21）

代入上面的表达式可简化为： [v0(r1r0)t]

kr1p(t)r0 （22）

r1

r1r2

用上式可求出p(t)k{v0(r1r0)t}C2，其中C2可通过p(0)k求得。那

么从9：00~12：00这三个小时内流出闸门的污染物的总量就可以求出来，污染物的量为：

Q

10800

p(t)r0dt （23）

Q'QQ （24）

（Q'为三小时内从1水库和2水库流出的污染物总量） 若Q'q则发生大面积污染； 若Q'q则不会发生大面积污染；

问题三：如果两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么问题二结果是否会改变？由题目分析可知，当两水库之间有一人工修建的水渠相互连通时，水渠中的水流必然是从高水位流向低水位，为了使模型简化，我们有以下说明： 1、由于两水库是连通的，因此在水库1关闸前，两水库的液面必然是趋近于等高的，否则必然有水从一个水库流向另一个水库。 2、在事故发生后一个小时内，考虑两水库的规模相当，且水库1没有关闭泄洪闸，此时两水库彼此不受影响。 3、在事故发生一个小时后，考虑到水库1要关闭泄洪闸，这势必引起水库1的水位上升，由说明1可知水库1中的水必将通过水渠流向水库2。 从上面的分析可知，在有连通水渠的情况下，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。下面我们就该问题给出进一步分析。 在第一个小时内：（即9：00~10：00） 此时，水库1中还未发生事故，只有水库2中在排放废料，用问题2的模型我们可求出流出的污染物的量为：

Q1

3600

p(t)r0dt （25）

（p(t)由（20）式确定）

在第二个小时内（即10：00~11：00） 水库1已发生装船事故，而水库2继续排放废料，但泄洪闸还未关闭，所以我们仍然独立考虑。 对水库1我们用问题一的模型求解得：

r1

r1r1r0

Q21[v0(r1r0)t]v0r1r0p/v0r0dt （26）



对水库2我们仍有问题二的模型求解得到：

3600

Q22

7200

3600

p(t)r0dt （p(t)由（20）确定） （27）

Q2Q21Q22 （28）

在第三个小时内（即11：00~12：00）； 水库1泄洪闸正在关闭，而水库2继续排放废料。这时水库1和水库2就要结合在一起考虑了，如下图所示：

图 一

由上面分析可列出以下的方程式： 对于水库1：

p(t)

其中v(t)的表达式为：

v1(t)v0r1t(2r0

dv(t)dt

v(t)

dp(t)dt

p(t)r0(t)p(t)r0''(t) （29）

r0t

)tr0''(t)t （3600t7200） 3600

而： r0(t)2r0对水库2： p(t)

'

r0t

（3600t7200） 3600

dv(t)dt

v(t)

dp'(t)dt

p1'(t)r1'(t)p(t)r0''(t)p'(t)r0'(t) （30）

其中

'

v2(t)v0r1'tr0''(t)tr0't

由基本假设知道，水库1和水库2规模相当则：

v(t)v0r1tr0''(t)tr0't

假定r0''(t)为一定值r0''。

'

联立解出：p(t)与p(t)：由于所得表达式非常复杂，我们把p(t)表达式放在附录（三）里表

示。而p'(t)由于过于复杂，我们这里就不给出解析解了。

Q31

7200

3600

p(t)r0(t)dt （3600t7200）

Q32

10800

7200

p'(t)r0'dt （31）

r0t

3600

r0(t)2r1

故有：

Q3Q31Q32 （32）

在12：00以后一个小时内，1水库完全关闭但1水库的水将通过渠道流入2水库内。则有以下式：

对水库1有：

p(t)其中v(t)表达式为： 对水库2有： p(t)其中v(t)表达式为：

'

dv(t)dt

v(t)

dp(t)dt

p(t)r0''(t) （33）

v(t)v0r1tr0''t

dv(t)dt

dp'(t)dt

p(t)r0''(t)p'(t)r0'(t) （34）

v(t)

v(t)v0r0''tr1'tr0't

14400

''

Qp(t)r求得 420dt （35） 10800

Q总Q1Q2Q3Q42 （36）

若Qq则干流不会发生大面积污染； 若Qq则干流会发生大面积污染；

问题四：如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模

型。

由问题的分析中知道，我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量mp为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。

4.2 模型二

为使模型清楚，我们先给出下面所用名词解释。

源和漏：是对体积元内污染物变化的一种描述，源是体积元内污染物的增加速率，漏是体积元内污染物的减少速率。这里的“漏”不意味着漏掉，而有更广泛的意义，如污染物的降解、沉淀和挥发等都属于“漏”。 由于污染物在干流中迁移符合迁移方程，由基本假设中我们知道河流中的水流处于推流状态，也就是体积元中水分子以同一速度向下游运动如图二和图三。设C(x)和C(x+x)分别为进入水片的水中某中污染物的浓度，mp是计算体积元中所含的该污染物的质量。按照质量守衡原理，计算体积元里污染物质量的增量为：

m[C(x)Q(x)C(xx)Q(xx)]tSxtSbxt

（1） （2） （3） （4）

+SvAxt （37） （5）

上式中（1）至（5）各项的意义如下 （1）——储存量项； （2）——平流输送项；

（3）——侧向的源和漏，其中Sl是单位时间内单位长度上的源和漏，它通常是由侧向分布流量q带入的，因此可把Sl写为：

mp(xAC) （38）

SlqCi

这里Ci是该分布流量中污染物的浓度；

（4）——表面的源和漏，Ss是单位时间内单位面积上的源和漏（g/m2s）； （5）——体积元内的源和漏，Sv是单位时间内单位体积内的源和漏（g/m3s）；

图 二

图 三

把（38）式带入（37）式，用xt除方程两边并令x0和t0，则得：

(AC)(QC)

SlbSsASv （39）

tx

式（10）就是推流时的污染物迁移方程，以后将广泛地利用这个方程，在求解问题四时，我们只要Sl、Ss和Sv的值为负数就可以了，即可以通过关闭人工水渠和部分关闭水库2泄洪闸的手段实现短时间内控制污染。

五、模型的评价与改进

5.1 模型的优点

本模型在建立模型前，由于水体运动非常复杂，污染物的性质也多种多样，污染物在水中的迁移情况受各种因素限制，不容易全部考虑。本文综合考虑上述因素，先通过合理的假设，首先从易溶急难挥发性污染物入手，考虑污染物和库水混合均匀的情况，利用物质平衡原理，针对问题中各种情况建立物质平衡方程模型。在此基础上考虑污染物的一维迁移模型。

在建立一维迁移扩展模型时，我们通过转化思想，将污染物对干流形成大面积污染的可能转换为求解在关闭泄洪闸的过程中，水库中多大面积水域内发生事故时才会对干流形成大面积污染，使得问题求解方案明了，为了验证模型的正确性，我们通过遗传算法对事例进行求解，很好的验证了模型。

5.2

模型的缺点

模型在建立过程中我们所考虑的此污染物为易溶、不挥发的物质，而没有考虑污染物为油性、挥发性和沉淀物质。另外，我们在考虑水体流动时没有考虑水体上、中、下层的不同流动情况。在问题四的求解时，我们只考虑推流问题，而实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题。时间有限，未能对其他方面做深入的探讨，使得模型十分完美。

5.3 模型的扩展

从问题的分析中可以知道，我们知道污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，为使模型比较符合实际，我们对上面模型进行扩展，建立污染物一维迁移扩展模型。 扩展模型一

从水库中流出的污染物的量可用下面的式子表示：

MQCt

由于在水库中污染物浓度在各个位置是不同的，而在整个水库中各个量又具有一定的连续性，故考虑t时刻在闸门口处t时间内通过的量为：mQCt，此式中我们考虑C与x，t的关系，其中x方向为水库流水的流向。这是因为对于事故泄漏排放的污染物，人们最关心的是其通过某一位置的时间，最大浓度等数值。因此，研究一维流场中的瞬时排放污染物的分布特征就非常必要了，我们用下式为基础来考虑分析泄漏时污染物在闸门口的浓度变化特征。

在瞬时排放条件下：CP/Q，QAUx，此时有下式：

(XUxt)2

exp[]exp(kt) （40） C(X,t)

4DtA4Dxtx

P

考虑在泄洪闸处X：X为泄洪闸处到污染处的距离，且k0（守恒污染物）

C(X0,t)

P

AUx02t0

2Dx0tUx0

2

X

(t0)2

Ux

exp] （41） 2

2t0

其中： t0

考虑在一个极其短暂的时间t内，在这个t内污染物的浓度可近似看作是不变的。

mQC(X0,t)t

则在t时刻有：

Q

P

AUx02t0

X(t0)2

Ux

exp]t2

2t0

通过对上式积分可求出2个小时内流出水库的污染物的量为：

M1

7200

M1

dm

Q

P

AUx02t0

X(t0)2

Ux

exp]dt （42） 2

2t0

这里P，A，Q，Ux0，Dx0均为已知，只是发生的地点未知。

2

通过上式我们可以求出在两个小时内通过的量为q吨时的X0，但由于式（42）的解析

3

解很难求出，我们就把它写成下面的形式：

7200

abs(

Q

P

AUx02t0

X02

(t)

Ux2

exp]dtq) （43） 2

32t0

把（2）式的相应的量代入后，通过遗传算法可求出的最小值，而（43）式的最小值为零，由此可求出满足0时所对应的X0，求出X0后，就可以得到发生大面积污染的可能的点所在的范围，如下图所示：

图 四

A1为发生大面积污染的点的集合。 A2为不发生大面积污染的点的集合。

因此对于问题一：发生大面积污染的可能为：

Z

A1

； （44）

A1A2

问题二：

由题目分析知道，若水库2中有一工厂在9：00—12：00违规进行排放，可知这样会增大发生大面积污染的可能性。因为另外一水库处于稳定流动状态，工厂排放的污染物为稳定的连续排放，环境中污染物的分布状况也是稳定的，综合考虑可有下式：

X

CC0erfc()

4Dmt

（45）

X

C0(1erf())

4Dmt

2QC1qC2

erf(z)其中：C0 ，

Qq我们假定工厂距泄洪闸的距离可以确定为X0 则： CC0erf(X0

'



z

e

t2

dt。

'

4Dmt

)

X

污染物到达泄洪闸的时间为：t0

Ux

'

'

则泄漏的量为：

M2'

t10800t'

QC0erfc(

X0

'

4Dmt

)dt （46）

由此我们可以得到在另外一水库当中有工厂进行周期性排放下流出的总量M为： MM1M2 （47） 即

7200

M

Q

P

AUx02t0

(t

exp[

X02

)Ux

2

2t0

]dt'

t

10800t'

QC0erfc(

X0

'

4Dmt

)dt

我们把上式转化为abs(Mq)的形式，这里我们可以看出的表达式与问题（1）中有所不同，这是因为我们考虑的时间已从原来的2小时变成3小时了，对于问题（二）我们仍然用问题（一）的解法，即通过遗传算法(程序见附录四)求出发生大面积污染时的X0。通过这个X0我们可求出发生大面积污染的点的集合A1'，进而求出发生大面积

A1'

污染的可能性为Z，A是水库的面积,可知是个定值。

A

'

在这里我们可以给出扩展模型一求解下的事例； 其中：ux30m/s

Q

=1513m3/s

Dx1.2m2/s

q0.7m3/s

C11.0103kg/m3C20.5kg/m3 X30000m部分数据见附录二。

通过MATLAB求解得到结果如下：

问题一发生大面积污染的点的集合类圆半径R21882m； 问题二发生大面积污染的点的集合类圆半径R21870m。

扩展二

从问题的分析及假设可以知道，我们的模型在建立过程中我们所考虑的只是一种污染物，即此污染物为易溶、不挥发的物质。实际上污染物可能为油性、挥发性和沉淀物

质。对于油性物质我们应该考虑该物质的密度是否大于水的密度，若油性污染物小于水的密度，因为油性物质不溶于水，该污染物应该是漂浮在水面上的，而且呈油膜状，因此我们就应考虑水库表层水的流动情况；若油性污染物的密度大于水的密度，则污染物将会沉如水底，那么我们就应该考虑深沉水体的水流问题；若油性污染物等于水的密度，那么污染物就会悬浮在水库水体中，因此我们就应该考虑整个水体上、中、下层的流动情况（因为不同水层的水流的流动情况是不同的）。对于挥发性的污染物，我们还需要考虑此污染物在水和空气中的分配比例。在考虑问题四的时候我们考虑的是推流问题，实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题而。为使模型更符合实际，模型必须在上述方面在做深入探讨。

参考文献

[1] W.金士博[联邦德国] 《水环境数学模型》 中国建筑工业出版社 1987.10 [2] 岳天祥 《资源环境数学模型手册》 科学出版社 2003.10 [3] 余常昭 《环境流体力学导论》 清华大学出版社 1992.10 [4] 刘振航 《数学建模》 中国人民大学出版社 2004.05 [5] 陈春云 郑彤 《环境系统数学模型》 化学工业出版社 2003.04 [6] 周建华 黄燕 《MATLAB5.3》 北京大学出版社 2000.12 [7] 丁春利 《精通MATLAB6》 清华大学出版社 2002.06 [8] 林芳荣等编译，面污染源管理与控制手册，广州：科学普及出版社广州分社， 15－43，1987，1，6。

[9] 方子云主编，水资源保护工作工作手册，河海大学出版社，1998，7。

附录

附录一：

水库实际库容及水域面积计算公式：

Vs(Hhi)

Ssn

式中，s为DEM象元面积；H为水库水位高程；hi为库域DEM各像元高程值；n为库域DEN像元数。

（曾永年，马海州，沙占江等；龙羊峡库区环境动态监测信息系统的建立与应用。遥感学报，2000.4⑵）

附录三：

1''

p(t)C1exp((45%2ln(1800v01800r1t3600r0tr1t21800r0t)%12

430r0%215%2r0)/%1)

%1:2v0r0900r13600r0r11800r1r03600r03600r0r0900r01''

r0t30r160r030r0

%2:)2

%1

2

''

2

''

''2

''

12

1

附录四： 遗传算法程序

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 bestfitness>newbestfitness

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.[***********][1**********]0);

%求最大值

%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 pick>pcross

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))-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

if pick>pcross

continue;

end

% random position of crossover

pick=rand;

while pick==0

pick=rand;

end

mask=2^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;

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 pick>pcross

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;

end

ret=chrom;

end

function 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 varibles

switch opts

case 'binary' % binary coding

for i=length(lenchrom):-1:1

data(i)=bitand(code,2^lenchrom(i)-1);

code=(code-data(i))/(2^lenchrom(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)=bitand(code,2^lenchrom(i)-1);

code=(code-data(i))/(2^lenchrom(i));

end

ret=bound(:,1)'+data./(2.^lenchrom-1).*(bound(:,2)-bound(:,1))';

case 'float' % float coding

ret=code;

end

function y=f1(t)

y=0.[***********][1**********]632*erfc((30000-30*t)./sqrt(4.8*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 : 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' % 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 pick>pmutation

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=rand;

while pick==0

pick=rand;

end

index=ceil(pick*sizepop);

pick=rand;

if pick>pmutation

continue;

end

% mutation position

pick=rand;

while pick==0

pick=rand;

end

pos=ceil(pick*sum(lenchrom));

v=chrom(i,pos);

v1=v-bound(pos,1);

v2=bound(pos,2)-v;

pick=rand;

if pick>0.5

delta=v2*(1-pick^((1-pop(1)/pop(2))^2)); chrom(i,pos)=v+delta;

else

delta=v1*(1-pick^((1-pop(1)/pop(2))^2)); chrom(i,pos)=v-delta;

end

end

ret=chrom;

end

function ret=select(individuals,sizepop,opts)

% In this function,it fulfils a selection among chromosomes

% individuals input : set of individuals

% sizepop input : size of population

% opts input : tag for choosing method of selection

% ret output : new set of individuals

switch opts

case 'roulette' % roulette wheel model

sumfitness=sum(1./individuals.fitness);

sumf=(1./individuals.fitness)./sumfitness;

index=[];

for i=1:sizepop

pick=rand;

while pick==0

pick=rand;

end

for i=1:sizepop

pick=pick-sumf(i);

if pick

index=[index i];

break;

end

end

end

individuals.chrom=individuals.chrom(index,:);

individuals.fitness=individuals.fitness(index);

ret=individuals;

case 'tournament' % tourment model allindex=[];

for i=1:sizepop

pick=rand(1,2);

while prod(pick)==0

pick=rand(1,2);

end

index=ceil(pick.*sizepop);

if individuals.fitness(index(1))

else

allindex=[allindex index(2)];

end

end

individuals.chrom=individuals.chrom(allindex,:);

individuals.fitness=individuals.fitness(allindex);

ret=individuals;

end

z=zeros(1);

z=quadl('f1',1000,9800);

disp

安徽工程大学 数学建模 课程设计论文

水库排污问题

班 级： 数学112 姓 名：学 号： 3110801204 指导老师： 周金明 成 绩： 完成日期： 2013年7月3日

摘 要

本文针对水库突发性事故排污问题，首先通过问题一建立二维水质污染物浓度模

2

型，只要考虑在事故发生到关闭水库的两个小时内，流出水库的污染物的质量小于q

3吨，给出了单个水库对干流造成大面积污染的可能性 Q

W1W2

；然后根据问题二、

2

q10003

三建立两水库排污模型，问题二由于两个水库之间没有联系，我们只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了，问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都是从高水位流向低水位，因此，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。分析了在另一水库有连续点源污染物排放及水流相互影响的情况下，两水库对干流造成大面积污染的可能性大小

Q

W1'W2'W3'

q1000

。

进一步针对第三种情况的发生，给出在短时间内控制污染的有效措施；且讨论了若污染物具有挥发性，上述各情况造成干流发生大面积污染的可能性大小，为水库事故性排污问题提供了有价值的理论依据。环境容量是环境科学的一个基本理论问题，也是环境管理中的一个重要的实际应用问题。在实践中，环境容量是环境目标管理的基本依据，是环境规划的主要约束条件，也是污染物总量控制的关键技术支持。从环境管理、监测与监督的角度出发，水环境容量是指水体在设计水文条件和规定的环境目标下所能容纳的最大污染物量。分析库区的污染状况，提出了水库环境容量的计算原则、设计水文条件和水质保护目标。根据岸边环境容量计算结果，提出了三峡水库污染混合区的控制标准。同时，推荐了三峡水库水环境容量的综合方案。

关键词：水库排污；污染物浓度；流量；水流速度；环境容量

某条江流上有2条支流，每条支流上都兴建了规模相当的水库。由于正处在雨水多发季节，因此两个水库都以一定规模的流量进行泄洪。某天晚上10：00，在其中的一个水库中发生了两船相撞的事故，而其中的一条船装载的p吨化学物质（这里的化学物质可以是具有挥发性的，也可能是急难挥发的）全部泄漏至水库中。当水上航运事故处置中心接获事故报告，立即要求该水库关闭水库泄洪闸，以免化学物质随洪水流入干流，发生更大规模的污染。水库闸门开始关闭时，已经处在事故发生后的1个小时，而水库闸门彻底关闭也需要1个小时的时间。 根据当地环境监测的有关规定，干流大面积污染的危险警戒值设为：三小时内q吨该化学物质发生泄漏。

阅读完以上资料后，结合查询资料，请做出合理假设并完成以下要求：

(1) 试建立合理的数学模型，讨论由于此次事故的发生，干流发生大面积污染的可能性；

(2) 如果在另外的一水库中有一化工厂违规排放废料。废料中同样含有该化学物质。该工厂为躲避环境监测站的监控，均在晚上9：00-12：00违规进行周期性排放。在这种情形下，讨论由于此次事故的发生，干流发生大面积污染的可能性；

(3) 如果以上两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么上述结果是否会改变？请给出说明，若有改变，则给出修正的模型及结果；

(4) 如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模型。

由题目的背景知道，事故发生时两水库都正在泄洪，因此此时水库中的水流速水库较快。而泄露到水库中的化学物质不论是具有挥发性的，还是急难挥发的，它们对干流污染的情况总是类同的，因此我们总可以认为污染物是易溶急难挥发性物质。为使我们的模型简单，我们可以先假设事故发生在水库1中，污染物在水库中的分布是符合零维迁移模型的，此时流入水库的污染物能以很快的速度与水库中的水均匀混合，水库中任何部份水体的污染状况都是一样的，污染程度与水体在水库中的位置无关。而实际上，污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，但在这样一个突发事件要求短时间内得到控制的问题中，我们总可以用污染物瞬时混合均匀状态模型来代替污染物一维迁移的过程。

在解答问题一时，我们只要考虑在事故发生到关闭水库的两个小时内，流出水库的

2

污染物的质量小于q吨即可。

3

问题二中提到的情况只是将干流的污染源从一个增加到两个，那么发生大面积污染的可能性就要增大。由于两个水库之间没有联系，我们只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。

问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都是从高水位流向低水位，因此，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。

分析问题四，干流已经发生大面积污染，对比水体污染物处理的各种手段，不论是化学、物理还是生物手段都不可能在短时间内除去污染物，因此我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量mp为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。

则我们可假设：

（1） 污染物为速溶物质，因此药品从船上流入水中的时间很少，可以忽略不计； （2） 污染物质从水库中一经流出就进入干流； （3） 水库和河流中的水流都是处于推流状态；

（4） 两水库事故发生条件相同，即两水库有相同的客观条件； （5） 被污染的水库关闭泄洪闸后不再有水流流入干流； （6） 不考虑生物等因素在水库泄洪过程中的作用，污染物除了流出外不因腐烂沉积等

手段从水中消失；

（7） 外界因素不对水库的体积变化产生影响，例如：雨水、地表径流、底下径流等； （8） 参与模型的变量是连续变化的，并且充分光滑； （9） 不考虑从不同的渠道流入与流出水库之间的区别，只考虑携带污染物的水流入水

库和水库中的水流出对水库污染程度的影响，因此可以把水库看成是单流入单流出的系统。

三、符号的约定

注意：部分符号见论文中说明。

四、模型的建立与求解

由问题的分析中知道，流入水库的污染物能以很快的速度与水库中的水均匀混合，也就是说水库中的污染状况在任何局部水体都是一样的，污染程度与水体在水库中的位置无关，因此我们可以建立下面模型。

4.1 模型一

问题一：此次事故的发生，干流发生大面积污染的可能性； 根据物质平衡原理和题目假设可知：

水库1中污染物的改变量 = 流入的污染物的量 － 流出污染物的量

于是对于充分小的t，在时间（t，t+t）内有：

p(tt)v(tt)p(t)v(t)[p1(t)r1(t)p0(t)r0(t)]t

两边同除t，并使t－>0得：

p(tt)v(tt)p(t)v(t)

p1(t)r1(t)p0(t)r0(t) （1）

t

现假设f(t)=p(t)×v(t)得：

df(t)d(p(t)v(t))dv(t)dp(t)f(tt)f(t)p(t)v(t)

tdtdtdtdt

即原式可以写为：

p(t)

dv(t)dt

v(t)

dp(t)dt

p1(t)r1(t)p0(t)r0(t) （2）

dp(t)dt

在水库1中发生撞船事故后，污染物处于非稳定排放即：

dv(t)dt

0，而由于水库闸

门的关闭也势必会引起水库中水的体积变化，故：0。现不考虑流入水库中的水

所含有与泄漏污染物相同物质的情况而带来的影响，即可看作p1(t)0，另外由问题分析中知道：流出的污染物的浓度应与水库中污染物浓度相同，即p(t)p0(t)这样对于问题一我们可以得到求解公式：

p(t)

dv(t)dt

v(t)

dp(t)dt

p(t)r0(t) （3）

进一步我们假设从水库中流出的水的流量初始值（从t=0时算起）为r0，在关闭闸

门的过程中，我们假定流量处于线性变化的趋势。这一假设是基于流量与过流面积为线

性关系上作出的，进一步可得：

r0...................(0t3600)

r0t r0(t) （4）

2r...........(3600t720)003600

从上面可看出r0(t)为分数函数，这主要是因为水库闸门关闭是在事故发生一小时后

作出的。现在有了r0(t)的变化的表达式，为了能求出p(t)的表达式。我们还要写出v(t)的表达式。首先我们假设水库的体积的初始值为v0（t=0时），v0值我们可以通过卫星定位系统及所建立的模型求出（具体卫星定位系统模型见附表）。而跟v(t)相关的还有r1(t)的值。我们假设r1(t)为一定值r1，则v(t)随随时间变化的关系式为：

v(t)v0r1tr0(t)t

由于r0(t)为分段函数可知：v(t)也响应的为分段函数，具体函数表达式为：

.........(0t3600)v0r1tr0t..........

r0t v(t) （5）

vrt(2r)t...........(3600t7200)0103600

把（4）式代入（3）式可以得到：

当0t3600秒时：

p(t)(r1r0)(v0r1tr0t)当3600t7200秒时：

p(t)(r12r0

dp(t)r0rtrt

t)[v0r1t(2r00)t]p(t)(2r00) （7） 18003600dt3600

dp(t)dt

p(t)r0 （6）

对（5）式化简有：

(v0r1tr0t)通过推导的出：

dp(t)dt

p(t)r1 （8）

p(t)e

r0[

1

ln(v0r1tr0t)]r1r0

c （9）

r1

r0r1

有已知条件可知：p(0)P/v0，故cp/v0v0经简化后：

r1

r1r0

p(t)[v0(r1r0)t]v0

r1

r1r0

p/v0

0t360时0 （10）

对（7）简化后得；

lnp(t)(r1)

dt

r

v0(r12r0)t0t2

3600

r0

a得： 3600

(

r0

)3600

tdt

v

v0(r12r0)t0t2

3600

（11）

设v0c，(r1-2r0)=b，当b24ac

lnp(t)(r)[

24acb2

2atb4acb2

)](

r01

)[ln(at2btc) 36002a

24acb

b22atb)] （12）

222a4acb4acb2atb4acb

2

设：z

2

)

最后得到：

p1(t)e

0.5lnat(2btc)(0.5b720a)0z

c1 （13）

保持连续性，当t=3600时，p1(3600)=p(3600)，此时可得到相应的C1值，但由于不确定因素很多，故确定C1不是很容易，这主要是缺少数据造成的。

当b24ac时得到（14）式

p1(t)e

1

1b

ln(at2btc)(7200a)z'22

C1

2atbb24ac

c （14） z'ln

2b24ac2atbb4ac

同样为保持连续性，要求当t=3600时，p1(3600)=p(3600)。 最后：

r1r1

v0(r1r0)tr1r0v0r1r0...................................0t3600(s)0.5ln(at2btc)(0.5b7200a)z2

p(t)ec..........当b4ac........3600t7200(s) （15）1

1b2'ln(atbtc)(7200a)z

222ec...............当b4ac.........3600t7200(s)1

那么t时间内流出水库的污染物的量便可表示为：

Qp(t)r0(t)t （16） 在0t3600（s）时流出的污染物的量为：

r1

r1r1r0

Q1[v0(r1r0)t]v0r1r0p/v0r0dt （17）



在3600t7200（s）流出的量为：

3600

Q23600p(t)(2r0

7200

r0t

)dt （18） 3600

流出的总量： QQ1Q2 （19） 再用Q与2/3q进行比较，便得出是否会发生大面积污染。

问题二：如果在另外的一水库中有一化工厂均在晚上9：00-12：00违规进行周期性排放同样含有该化学物质的废料。讨论由于此次事故的发生，干流发生大面积污染的可能性。 由题目分析中知道，水库2中有一化工厂违规排放含有该污染物的废料，而由于两个水库之间没有联系，故我们只需要单独考虑水库2的排污量，然后加上水库1的污染物排放量，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。下面我们将在水库1模型的基础上建立2水库的模型。 p(t)

dv(t)dt

v(t)

dp(t)dt

p1(t)r1(t)p0(t)r0(t) （20）

设p1(t)k（常数），即在9：00~12：00这段时间内污染物以一个恒定值流入水库，考虑水库水的流入速度为一定值r1，流出速度也为一个定值r0

这样水库体积的表达式可写成如下的公式：

v(t)v0r1tr0t

dp(t)dt

(0t1080) 0 （21）

代入上面的表达式可简化为： [v0(r1r0)t]

kr1p(t)r0 （22）

r1

r1r2

用上式可求出p(t)k{v0(r1r0)t}C2，其中C2可通过p(0)k求得。那

么从9：00~12：00这三个小时内流出闸门的污染物的总量就可以求出来，污染物的量为：

Q

10800

p(t)r0dt （23）

Q'QQ （24）

（Q'为三小时内从1水库和2水库流出的污染物总量） 若Q'q则发生大面积污染； 若Q'q则不会发生大面积污染；

问题三：如果两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么问题二结果是否会改变？由题目分析可知，当两水库之间有一人工修建的水渠相互连通时，水渠中的水流必然是从高水位流向低水位，为了使模型简化，我们有以下说明： 1、由于两水库是连通的，因此在水库1关闸前，两水库的液面必然是趋近于等高的，否则必然有水从一个水库流向另一个水库。 2、在事故发生后一个小时内，考虑两水库的规模相当，且水库1没有关闭泄洪闸，此时两水库彼此不受影响。 3、在事故发生一个小时后，考虑到水库1要关闭泄洪闸，这势必引起水库1的水位上升，由说明1可知水库1中的水必将通过水渠流向水库2。 从上面的分析可知，在有连通水渠的情况下，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。下面我们就该问题给出进一步分析。 在第一个小时内：（即9：00~10：00） 此时，水库1中还未发生事故，只有水库2中在排放废料，用问题2的模型我们可求出流出的污染物的量为：

Q1

3600

p(t)r0dt （25）

（p(t)由（20）式确定）

在第二个小时内（即10：00~11：00） 水库1已发生装船事故，而水库2继续排放废料，但泄洪闸还未关闭，所以我们仍然独立考虑。 对水库1我们用问题一的模型求解得：

r1

r1r1r0

Q21[v0(r1r0)t]v0r1r0p/v0r0dt （26）



对水库2我们仍有问题二的模型求解得到：

3600

Q22

7200

3600

p(t)r0dt （p(t)由（20）确定） （27）

Q2Q21Q22 （28）

在第三个小时内（即11：00~12：00）； 水库1泄洪闸正在关闭，而水库2继续排放废料。这时水库1和水库2就要结合在一起考虑了，如下图所示：

图 一

由上面分析可列出以下的方程式： 对于水库1：

p(t)

其中v(t)的表达式为：

v1(t)v0r1t(2r0

dv(t)dt

v(t)

dp(t)dt

p(t)r0(t)p(t)r0''(t) （29）

r0t

)tr0''(t)t （3600t7200） 3600

而： r0(t)2r0对水库2： p(t)

'

r0t

（3600t7200） 3600

dv(t)dt

v(t)

dp'(t)dt

p1'(t)r1'(t)p(t)r0''(t)p'(t)r0'(t) （30）

其中

'

v2(t)v0r1'tr0''(t)tr0't

由基本假设知道，水库1和水库2规模相当则：

v(t)v0r1tr0''(t)tr0't

假定r0''(t)为一定值r0''。

'

联立解出：p(t)与p(t)：由于所得表达式非常复杂，我们把p(t)表达式放在附录（三）里表

示。而p'(t)由于过于复杂，我们这里就不给出解析解了。

Q31

7200

3600

p(t)r0(t)dt （3600t7200）

Q32

10800

7200

p'(t)r0'dt （31）

r0t

3600

r0(t)2r1

故有：

Q3Q31Q32 （32）

在12：00以后一个小时内，1水库完全关闭但1水库的水将通过渠道流入2水库内。则有以下式：

对水库1有：

p(t)其中v(t)表达式为： 对水库2有： p(t)其中v(t)表达式为：

'

dv(t)dt

v(t)

dp(t)dt

p(t)r0''(t) （33）

v(t)v0r1tr0''t

dv(t)dt

dp'(t)dt

p(t)r0''(t)p'(t)r0'(t) （34）

v(t)

v(t)v0r0''tr1'tr0't

14400

''

Qp(t)r求得 420dt （35） 10800

Q总Q1Q2Q3Q42 （36）

若Qq则干流不会发生大面积污染； 若Qq则干流会发生大面积污染；

问题四：如果发生了大面积污染，那么针对第三种情况，试给出在短时间内控制污染模

型。

由问题的分析中知道，我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量mp为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。

4.2 模型二

为使模型清楚，我们先给出下面所用名词解释。

源和漏：是对体积元内污染物变化的一种描述，源是体积元内污染物的增加速率，漏是体积元内污染物的减少速率。这里的“漏”不意味着漏掉，而有更广泛的意义，如污染物的降解、沉淀和挥发等都属于“漏”。 由于污染物在干流中迁移符合迁移方程，由基本假设中我们知道河流中的水流处于推流状态，也就是体积元中水分子以同一速度向下游运动如图二和图三。设C(x)和C(x+x)分别为进入水片的水中某中污染物的浓度，mp是计算体积元中所含的该污染物的质量。按照质量守衡原理，计算体积元里污染物质量的增量为：

m[C(x)Q(x)C(xx)Q(xx)]tSxtSbxt

（1） （2） （3） （4）

+SvAxt （37） （5）

上式中（1）至（5）各项的意义如下 （1）——储存量项； （2）——平流输送项；

（3）——侧向的源和漏，其中Sl是单位时间内单位长度上的源和漏，它通常是由侧向分布流量q带入的，因此可把Sl写为：

mp(xAC) （38）

SlqCi

这里Ci是该分布流量中污染物的浓度；

（4）——表面的源和漏，Ss是单位时间内单位面积上的源和漏（g/m2s）； （5）——体积元内的源和漏，Sv是单位时间内单位体积内的源和漏（g/m3s）；

图 二

图 三

把（38）式带入（37）式，用xt除方程两边并令x0和t0，则得：

(AC)(QC)

SlbSsASv （39）

tx

式（10）就是推流时的污染物迁移方程，以后将广泛地利用这个方程，在求解问题四时，我们只要Sl、Ss和Sv的值为负数就可以了，即可以通过关闭人工水渠和部分关闭水库2泄洪闸的手段实现短时间内控制污染。

五、模型的评价与改进

5.1 模型的优点

本模型在建立模型前，由于水体运动非常复杂，污染物的性质也多种多样，污染物在水中的迁移情况受各种因素限制，不容易全部考虑。本文综合考虑上述因素，先通过合理的假设，首先从易溶急难挥发性污染物入手，考虑污染物和库水混合均匀的情况，利用物质平衡原理，针对问题中各种情况建立物质平衡方程模型。在此基础上考虑污染物的一维迁移模型。

在建立一维迁移扩展模型时，我们通过转化思想，将污染物对干流形成大面积污染的可能转换为求解在关闭泄洪闸的过程中，水库中多大面积水域内发生事故时才会对干流形成大面积污染，使得问题求解方案明了，为了验证模型的正确性，我们通过遗传算法对事例进行求解，很好的验证了模型。

5.2

模型的缺点

模型在建立过程中我们所考虑的此污染物为易溶、不挥发的物质，而没有考虑污染物为油性、挥发性和沉淀物质。另外，我们在考虑水体流动时没有考虑水体上、中、下层的不同流动情况。在问题四的求解时，我们只考虑推流问题，而实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题。时间有限，未能对其他方面做深入的探讨，使得模型十分完美。

5.3 模型的扩展

从问题的分析中可以知道，我们知道污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，为使模型比较符合实际，我们对上面模型进行扩展，建立污染物一维迁移扩展模型。 扩展模型一

从水库中流出的污染物的量可用下面的式子表示：

MQCt

由于在水库中污染物浓度在各个位置是不同的，而在整个水库中各个量又具有一定的连续性，故考虑t时刻在闸门口处t时间内通过的量为：mQCt，此式中我们考虑C与x，t的关系，其中x方向为水库流水的流向。这是因为对于事故泄漏排放的污染物，人们最关心的是其通过某一位置的时间，最大浓度等数值。因此，研究一维流场中的瞬时排放污染物的分布特征就非常必要了，我们用下式为基础来考虑分析泄漏时污染物在闸门口的浓度变化特征。

在瞬时排放条件下：CP/Q，QAUx，此时有下式：

(XUxt)2

exp[]exp(kt) （40） C(X,t)

4DtA4Dxtx

P

考虑在泄洪闸处X：X为泄洪闸处到污染处的距离，且k0（守恒污染物）

C(X0,t)

P

AUx02t0

2Dx0tUx0

2

X

(t0)2

Ux

exp] （41） 2

2t0

其中： t0

考虑在一个极其短暂的时间t内，在这个t内污染物的浓度可近似看作是不变的。

mQC(X0,t)t

则在t时刻有：

Q

P

AUx02t0

X(t0)2

Ux

exp]t2

2t0

通过对上式积分可求出2个小时内流出水库的污染物的量为：

M1

7200

M1

dm

Q

P

AUx02t0

X(t0)2

Ux

exp]dt （42） 2

2t0

这里P，A，Q，Ux0，Dx0均为已知，只是发生的地点未知。

2

通过上式我们可以求出在两个小时内通过的量为q吨时的X0，但由于式（42）的解析

3

解很难求出，我们就把它写成下面的形式：

7200

abs(

Q

P

AUx02t0

X02

(t)

Ux2

exp]dtq) （43） 2

32t0

把（2）式的相应的量代入后，通过遗传算法可求出的最小值，而（43）式的最小值为零，由此可求出满足0时所对应的X0，求出X0后，就可以得到发生大面积污染的可能的点所在的范围，如下图所示：

图 四

A1为发生大面积污染的点的集合。 A2为不发生大面积污染的点的集合。

因此对于问题一：发生大面积污染的可能为：

Z

A1

； （44）

A1A2

问题二：

由题目分析知道，若水库2中有一工厂在9：00—12：00违规进行排放，可知这样会增大发生大面积污染的可能性。因为另外一水库处于稳定流动状态，工厂排放的污染物为稳定的连续排放，环境中污染物的分布状况也是稳定的，综合考虑可有下式：

X

CC0erfc()

4Dmt

（45）

X

C0(1erf())

4Dmt

2QC1qC2

erf(z)其中：C0 ，

Qq我们假定工厂距泄洪闸的距离可以确定为X0 则： CC0erf(X0

'



z

e

t2

dt。

'

4Dmt

)

X

污染物到达泄洪闸的时间为：t0

Ux

'

'

则泄漏的量为：

M2'

t10800t'

QC0erfc(

X0

'

4Dmt

)dt （46）

由此我们可以得到在另外一水库当中有工厂进行周期性排放下流出的总量M为： MM1M2 （47） 即

7200

M

Q

P

AUx02t0

(t

exp[

X02

)Ux

2

2t0

]dt'

t

10800t'

QC0erfc(

X0

'

4Dmt

)dt

我们把上式转化为abs(Mq)的形式，这里我们可以看出的表达式与问题（1）中有所不同，这是因为我们考虑的时间已从原来的2小时变成3小时了，对于问题（二）我们仍然用问题（一）的解法，即通过遗传算法(程序见附录四)求出发生大面积污染时的X0。通过这个X0我们可求出发生大面积污染的点的集合A1'，进而求出发生大面积

A1'

污染的可能性为Z，A是水库的面积,可知是个定值。

A

'

在这里我们可以给出扩展模型一求解下的事例； 其中：ux30m/s

Q

=1513m3/s

Dx1.2m2/s

q0.7m3/s

C11.0103kg/m3C20.5kg/m3 X30000m部分数据见附录二。

通过MATLAB求解得到结果如下：

问题一发生大面积污染的点的集合类圆半径R21882m； 问题二发生大面积污染的点的集合类圆半径R21870m。

扩展二

从问题的分析及假设可以知道，我们的模型在建立过程中我们所考虑的只是一种污染物，即此污染物为易溶、不挥发的物质。实际上污染物可能为油性、挥发性和沉淀物

质。对于油性物质我们应该考虑该物质的密度是否大于水的密度，若油性污染物小于水的密度，因为油性物质不溶于水，该污染物应该是漂浮在水面上的，而且呈油膜状，因此我们就应考虑水库表层水的流动情况；若油性污染物的密度大于水的密度，则污染物将会沉如水底，那么我们就应该考虑深沉水体的水流问题；若油性污染物等于水的密度，那么污染物就会悬浮在水库水体中，因此我们就应该考虑整个水体上、中、下层的流动情况（因为不同水层的水流的流动情况是不同的）。对于挥发性的污染物，我们还需要考虑此污染物在水和空气中的分配比例。在考虑问题四的时候我们考虑的是推流问题，实际污染随水流迁移是以水团形式的，这样我们就需要考虑水体不同水层的流动问题而。为使模型更符合实际，模型必须在上述方面在做深入探讨。

参考文献

[1] W.金士博[联邦德国] 《水环境数学模型》 中国建筑工业出版社 1987.10 [2] 岳天祥 《资源环境数学模型手册》 科学出版社 2003.10 [3] 余常昭 《环境流体力学导论》 清华大学出版社 1992.10 [4] 刘振航 《数学建模》 中国人民大学出版社 2004.05 [5] 陈春云 郑彤 《环境系统数学模型》 化学工业出版社 2003.04 [6] 周建华 黄燕 《MATLAB5.3》 北京大学出版社 2000.12 [7] 丁春利 《精通MATLAB6》 清华大学出版社 2002.06 [8] 林芳荣等编译，面污染源管理与控制手册，广州：科学普及出版社广州分社， 15－43，1987，1，6。

[9] 方子云主编，水资源保护工作工作手册，河海大学出版社，1998，7。

附录

附录一：

水库实际库容及水域面积计算公式：

Vs(Hhi)

Ssn

式中，s为DEM象元面积；H为水库水位高程；hi为库域DEM各像元高程值；n为库域DEN像元数。

（曾永年，马海州，沙占江等；龙羊峡库区环境动态监测信息系统的建立与应用。遥感学报，2000.4⑵）

附录三：

1''

p(t)C1exp((45%2ln(1800v01800r1t3600r0tr1t21800r0t)%12

430r0%215%2r0)/%1)

%1:2v0r0900r13600r0r11800r1r03600r03600r0r0900r01''

r0t30r160r030r0

%2:)2

%1

2

''

2

''

''2

''

12

1

附录四： 遗传算法程序

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 bestfitness>newbestfitness

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.[***********][1**********]0);

%求最大值

%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 pick>pcross

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))-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

if pick>pcross

continue;

end

% random position of crossover

pick=rand;

while pick==0

pick=rand;

end

mask=2^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;

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 pick>pcross

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;

end

ret=chrom;

end

function 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 varibles

switch opts

case 'binary' % binary coding

for i=length(lenchrom):-1:1

data(i)=bitand(code,2^lenchrom(i)-1);

code=(code-data(i))/(2^lenchrom(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)=bitand(code,2^lenchrom(i)-1);

code=(code-data(i))/(2^lenchrom(i));

end

ret=bound(:,1)'+data./(2.^lenchrom-1).*(bound(:,2)-bound(:,1))';

case 'float' % float coding

ret=code;

end

function y=f1(t)

y=0.[***********][1**********]632*erfc((30000-30*t)./sqrt(4.8*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 : 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' % 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 pick>pmutation

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=rand;

while pick==0

pick=rand;

end

index=ceil(pick*sizepop);

pick=rand;

if pick>pmutation

continue;

end

% mutation position

pick=rand;

while pick==0

pick=rand;

end

pos=ceil(pick*sum(lenchrom));

v=chrom(i,pos);

v1=v-bound(pos,1);

v2=bound(pos,2)-v;

pick=rand;

if pick>0.5

delta=v2*(1-pick^((1-pop(1)/pop(2))^2)); chrom(i,pos)=v+delta;

else

delta=v1*(1-pick^((1-pop(1)/pop(2))^2)); chrom(i,pos)=v-delta;

end

end

ret=chrom;

end

function ret=select(individuals,sizepop,opts)

% In this function,it fulfils a selection among chromosomes

% individuals input : set of individuals

% sizepop input : size of population

% opts input : tag for choosing method of selection

% ret output : new set of individuals

switch opts

case 'roulette' % roulette wheel model

sumfitness=sum(1./individuals.fitness);

sumf=(1./individuals.fitness)./sumfitness;

index=[];

for i=1:sizepop

pick=rand;

while pick==0

pick=rand;

end

for i=1:sizepop

pick=pick-sumf(i);

if pick

index=[index i];

break;

end

end

end

individuals.chrom=individuals.chrom(index,:);

individuals.fitness=individuals.fitness(index);

ret=individuals;

case 'tournament' % tourment model allindex=[];

for i=1:sizepop

pick=rand(1,2);

while prod(pick)==0

pick=rand(1,2);

end

index=ceil(pick.*sizepop);

if individuals.fitness(index(1))

else

allindex=[allindex index(2)];

end

end

individuals.chrom=individuals.chrom(allindex,:);

individuals.fitness=individuals.fitness(allindex);

ret=individuals;

end

z=zeros(1);

z=quadl('f1',1000,9800);

disp