一种水库排沙发电多目标优化调度方法

摘要

一种水库排沙发电多目标优化调度方法，包括以下步骤：一、以系统发电量最大为目标建立目标函数；以水库输沙量最大为目标建立目标函数；二、将下泄流量作为决策变量，产生多个序列，利用遗传算法求解决策变量；通过上述步骤得到最优的下泄流量。在求解决策变量时，引入水库系统运行的约束条件，包括库水位约束、水库下泄流量约束和流量日变幅约束。本发明提供的一种水库排沙发电多目标优化调度方法，通过综合考虑排沙效率最优和发电效益最优这两个目标函数，实现对水库群长期运行综合效益最大化的多目标表达，并采用多目标优化遗传算法，给出了水库调度的Pareto最优解集。

权利要求

1.一种水库排沙发电多目标优化调度方法，其特征是包括以下步骤：一、以系统发电量最大为目标建立目标函数；

系统发电量最大目标函数为：

其中，E是年发电总量； 是编号第i台发电机组在第t个时段内的单机平均出力；Δt为计算时长；i为水电站的各机组编号；ni为各机组的台数；t代表时段长度；

以水库输沙量最大为目标建立目标函数；

水库输沙量最大目标函数为：

其中，S是水库的输沙量，单位为kg；QFt为t时段长度的泄洪流量，单位为m3/s， 为泄洪3

口平均含沙率，单位为kg/m；

泄洪口平均含沙率由下式计算得到：

式中，n为泄洪口个数，y为当前高程，h为泄洪口底部高程，L为泄洪口宽度，H为泄洪口高，Sva为悬移质泥沙在距床面为a处的体积比含沙浓度；Z为悬浮指标；

泄洪流量由下式计算得到：

式中，为泄洪口平均流速；

二、将下泄流量作为决策变量，产生多个基于目标函数的序列，利用基于非支配解的多目标遗传算法求解决策变量；

在求解决策变量时，引入水库系统运行的约束条件，包括库水位约束、水库下泄流量约束和流量日变幅约束；

随机产生p个染色体，即产生p个水库下泄流量Q的序列，组成初始种群P0，种群规模＝p；

以一个时间段为划分标准，相应每个染色体具有多个的发电引用流量元素和多个泄洪流量元素的基因，这些基因组成数组；

通过上述步骤得到最优的下泄流量。

2.根据权利要求1所述的一种水库排沙发电多目标优化调度方法，其特征是：以月为单位，每个染色体上有12个发电引用流量元素和4个泄洪流量元素共16个基因，组成16×p的数组。

3.根据权利要求1或2所述的一种水库排沙发电多目标优化调度方法，其特征是：在p个染色体中选取较优的染色体的m个值进行交叉变异计算，生成新种群P1；

计算新种群中的染色体对应的目标函数值，进入下一次迭代。

4.根据权利要求3所述的一种水库排沙发电多目标优化调度方法，其特征是：遍历种群中任意两对染色体，比较两个染色体的适应度函数值，选择较优的染色体；

将当前染色体的m个值进入交叉变异；

若m＝0，重新生成初始种群；若m＞0，进入交叉变异；

根据设定的交叉概率，若交叉发生，则从m个值中任意选择两个染色体，随机生成一个基因个数区间内的非零整数，大于此数的基因位置发生交换；

根据设定的变异概率，若变异发生，则对于上述染色体，随机生成一个基因个数区间内的非零整数，在该位置上的基因突变成m个染色体中该位置上最大值和最小值之间的任意数值；

交叉变异后新生成的p个染色体，与上一次迭代得到的较优的染色体一起构成新种群P1，种群规模为m+p；

计算新种群P1中对应的目标函数值，进入下一次迭代。

5.根据权利要求1所述的一种水库排沙发电多目标优化调度方法，其特征是：保持水库运行多目标中的矛盾性和不可公度性，而非划归为单目标求解；所得最优运行策略为一组最优解集，其中各最优解之间不可相互比较，即该最优解集中任意两个解之间比较时，若某个解的目标之一优于另一个解，则其目标之二必劣于另一个解。

说明书

一种水库排沙发电多目标优化调度方法

技术领域

[0001] 本发明属于水利水电领域的水力调度领域，特别是一种水库排沙发电多目标优化调度方法。

背景技术

[0002] 河流建设水库以后，由于水位升高、过水面积加大、流速减缓，从而使水流挟沙能力降低，导致水库内发生淤积。而挟沙能力与流速的高次方成比例，过水断面的些许改变，常会引起挟沙能力大幅度降低，对水库综合效益的发挥造成非常不利的影响。主要表现为：

[0003] 1、水库有效库容不断减少，减少了水库的使用寿命，影响水库的调能力。兴利效益逐年下降，使得水库防洪、发电、灌溉和供水等兴利效益大打折扣。

[0004] 2、受水库蓄水的影响，水流进入库区后，水流中的泥沙会在水库的回水末端，并且逐渐发展到上游，使库区尾水位产生的“翅尾巴”现象。抬高上游的河床高程和水库周围的地下水位，这不仅会降低上游河道的防洪效益，使通航条件恶化，而且会对两岸居民生活造成不利影响，威胁当地的农业生产安全。

[0005] 3、泥沙会造成水力机械和泄流设施磨损，影响机械效率，缩短水力机械寿命，增加材料损耗和维修时间，加速工程的老化损坏。

[0006] 4、下游河道原有的水沙条件发生变化。在非汛期，清水下泄，冲刷下游河道；水库进入汛期后，泥沙集中下泄，又会使下游回流淤积，使坝下游河道发生很大改变，对河道整治及两岸引水工作造成不利影响。

[0007] 因此，水库的调度运行不仅需要调节径流，蓄盈补缺，满足兴利、防洪等各方面的要求，还要调节泥沙，减少有效库容游损，使水库长期，有效，安全的发挥效益。其本质上是一个具有复杂约束条件的多目标优化问题。目前存在的多目标优化方法主要分为两大类：传统多目标优化方法和多目标智能优化方法。

[0008] 传统多目标优化方法以约束法、层次分析法、线性加权法、理想点法为主要的处理手段，根据决策者的主观因素，将多目标优化问题转化为单目标优化问题进行求解。其本质上还是单目标优化问题。

[0009] 多目标优化问题的难点实质在于目标间的不可公度性和矛盾性。所谓不可公度性是指各个目标没有统一的度量标准，因而难以进行比较。例如：在水库水沙联调的多目标决策中，发电量的单位为度，挟沙力的单位为kg/m3，无法进行比较。目标间的矛盾性是指如果采用某一种方案去改进某一目标值，可能会使另一目标的值变好或变坏。还以水沙联调为例，水库尽量保持最高水位，减少泄流，可以增大单位水体的发电量，但泄流降低将会使得水库泥沙淤积加重。增发电量和减少淤积是一对矛盾的运行目标。

[0010] 由于多个目标之间的矛盾性和不可公度性，不能简单地把多个目标归为单个目标，并使用单目标决策问题的方法去求解多目标决策问题。

发明内容

[0011] 本发明所要解决的技术问题是提供一种水库排沙发电多目标优化调度方法，可以实现对水库长期运行综合效益最大化的调度。

[0012] 为解决上述技术问题，本发明所采用的技术方案是：一种水库排沙发电多目标优化调度方法，包括以下步骤：

[0013] 一、以系统发电量最大为目标建立目标函数；

[0014] 以水库输沙量最大为目标建立目标函数；

[0015] 二、将下泄流量作为决策变量，产生多个基于目标函数的序列，利用基于非支配解的多目标遗传算法求解决策变量；

[0016] 通过上述步骤得到最优的下泄流量。

[0017] 在求解决策变量时，引入水库系统运行的约束条件，包括库水位约束、水库下泄流量约束和流量日变幅约束。

[0018] 优选的，步骤2中，随机产生p个染色体，即产生p个水库下泄流量Q的序列，组成初始种群P0，种群规模＝p；

[0019] 优选的，以一个时间段为划分标准，相应每个染色体具有多个的发电引用流量元素和多个泄洪流量元素的基因，这些基因组成数组。

[0020] 优选的，以月为单位，每个染色体上有12个发电引用流量元素和4个泄洪流量元素共16个基因，组成16×p的数组。

[0021] 进一步优选的，在p个染色体中选取较优的染色体的m个值进行交叉变异计算，生成新种群P1；

[0022] 计算新种群中的染色体对应的目标函数值，进入下一次迭代。

[0023] 进一步优选的，遍历种群中任意两对染色体，比较两个染色体的适应度函数值，选择较优的染色体；

[0024] 将当前染色体的m个值进入交叉变异；

[0025] 若m＝0，重新生成初始种群；若m＞0，进入交叉变异；

[0026] 根据设定的交叉概率，若交叉发生，则从m个值中任意选择两个染色体，随机生成一个基因个数区间内的非零整数，大于此数的基因位置发生交换；

[0027] 根据设定的变异概率，若变异发生，则对于上述染色体，随机生成一个基因个数区间内的非零整数，在该位置上的基因突变成m个染色体中该位置上最大值和最小值之间的任意数值；

[0028] 交叉变异后新生成的p个染色体，与上一次迭代得到的较优的染色体一起构成新种群P1，种群规模为m+p；

[0029] 计算新种群P1中对应的目标函数值，进入下一次迭代。

[0030] 保持水库运行多目标中的矛盾性和不可公度性，而非划归为单目标求解；所得最优运行策略为一组最优解集，其中各最优解之间不可相互比较，即该最优解集中任意两个解之间比较时，若某个解的目标之一优于另一个解，则其目标之二必劣于另一个解。

[0031] 系统发电量最大目标函数为：

[0032]

[0033] 其中，E是年发电总量； 是编号第i台发电机组在第t个时段内的单机平均出力；Δt为计算时长；i为水电站的各机组编号；ni为各机组的台数；t代表时段长度。

[0034] 水库输沙量最大目标函数为：

[0035]

[0036] 其中，S是水库的输沙量，单位为kg；QFt为t时段长度的泄洪流量，单位为m3/s， 为泄洪口平均含沙率，单位为kg/m3。

[0037] 泄洪口平均含沙率由下式计算得到：

[0038]

[0039] 式中，n为泄洪口个数，y为当前高程，h为泄洪口底部高程，L为泄洪口宽度，H为泄洪口高。Z为悬浮指标；

[0040] 泄洪流量由下式计算得到：

[0041]

[0042] 式中，n为泄洪口个数，L为泄洪口宽度，H为泄洪口高，为泄洪口平均流速。

[0043] 发明人发现多目标问题对应的最优解不是某个单一的最优解，而是一个集合的概念，是一组Pareto最优解。

[0044] 发明人对多目标优化算法的研究发现，进化算法非常适合求解这类多目标优化问题。进化算法是将每一个可能的解作为群体中的一个个体元素，通过简单的交叉、变异等操作对群体进行进化，把生存和繁殖的机会留给适应性较强的个体。它包括遗传算法、蚁群算法、粒子群算法等。进化算法，与其它智能计算方法相比，更适合用于处理那些不能用传统方法求解的复杂优化问题。随着智能计算的发展，具有并行特性的进化算法体现出对解决多目标优化问题具有很大优势。

[0045] 本发明提供的一种水库排沙发电多目标优化调度方法，通过综合考虑排沙效率最优和发电效益最优这两个目标函数，实现对水库群长期运行综合效益最大化的多目标表达，并采用多目标优化遗传算法，给出了水库调度的Pareto最优解集。

附图说明

[0046] 下面结合附图和实施例对本发明作进一步说明：

[0047] 图1为本发明中坝前汛期悬移质含沙量垂线分布曲线。

[0048] 图2为本发明中坝前汛期流速垂线分布曲线。

[0049] 图3为本发明中迭代后非支配解对应的目标函数曲线，图中日均流量变化幅值为3
50m/s。

[0050] 图4为本发明中迭代后非支配解对应的目标函数曲线，图中日均流量变化幅值为150m3/s。

[0051] 图5为本发明中迭代后非支配解对应的目标函数曲线，图中日均流量变化幅值为1000m3/s。

[0052] 图6为本发明中迭代后非支配解对应的目标函数曲线，图中日均流量变化幅值为5000m3/s。

具体实施方式

[0053] 实施例1：

[0054] 一种水库排沙发电多目标优化调度方法，包括以下步骤：

[0055] 一、以系统发电量最大为目标建立目标函数；

[0056] 以水库输沙量最大为目标建立目标函数；

[0057] 二、将下泄流量作为决策变量，产生多个基于上述目标函数的序列，利用基于非支配解的多目标遗传算法求解决策变量；

[0058] 通过上述步骤得到最优的下泄流量。

[0059] 在求解决策变量时，引入水库系统运行的约束条件，包括库水位约束、水库下泄流量约束和流量日变幅约束。

[0060] 优选的，步骤2中，随机产生n个染色体，即产生n个水库下泄流量(Q)的序列，组成初始种群P0，种群规模＝n；

[0061] 优选的，以一个时间段为划分标准，例如月或周，相应每个染色体具有多个的发电引用流量元素和多个泄洪流量元素的基因，这些基因组成数组。

[0062] 优选的，以月为单位，每个染色体上有12个发电引用流量元素和4个泄洪流量元素共16个基因，组成16×n的数组。

[0063] 进一步优选的，在n个染色体中选取较优的染色体的m个值进行交叉变异计算，生成新种群P1；

[0064] 计算新种群中的染色体对应的目标函数值，进入下一次迭代。

[0065] 进一步优选的，遍历种群中任意两对染色体，比较两个染色体的适应度函数值，选择较优的染色体；

[0066] 将当前染色体的m个值进入交叉变异；

[0067] 若m＝0，重新生成初始种群；若m＞0，进入交叉变异；

[0068] 根据设定的交叉概率，若交叉发生，则从m个值中任意选择两个染色体，随机生成一个基因个数区间内的非零整数，大于此数的基因位置发生交换；

[0069] 根据设定的变异概率，若变异发生，则对于上述染色体，随机生成一个基因个数区间内的非零整数，在该位置上的基因突变成m个染色体中该位置上最大值和最小值之间的任意数值；

[0070] 交叉变异后新生成的n个染色体，与上一次迭代得到的较优的染色体一起构成新种群P1，种群规模为m+n；

[0071] 计算新种群P1中对应的目标函数值，进入下一次迭代。

[0072] 优选的，所述的交叉概率为0.6～0.8；

[0073] 所述的变异概率为0.01～0.03。

[0074] 保持水库运行多目标中的矛盾性和不可公度性，而非划归为单目标求解；所得最优运行策略为一组最优解集，其中各最优解之间不可相互比较，即该最优解集中任意两个解之间比较时，若某个解的目标之一优于另一个解，则其目标之二必劣于另一个解。

[0075] 系统发电量最大目标函数为：

[0076]

[0077] 其中，E是年发电总量；是单机平均出力；Δt为计算时长；i为水电站的各机组编号；ni为各机组的台数；t代表时段长度。

[0078] 水库输沙量最大目标函数为：

[0079]

[0080] 其中，S是水库的输沙量，单位为kg；QF为泄洪流量，单位为m3/s， 为泄洪口平均含沙率，单位为kg/m3。

[0081] 泄洪口平均含沙率由下式计算得到：

[0082]

[0083] 式中，n为泄洪口个数，y为当前高程，h为泄洪口底部高程，L为泄洪口宽度，H为泄洪口高。

[0084] 泄洪流量由下式计算得到：

[0085]

[0086] 式中，n为泄洪口个数，L为泄洪口宽度，H为泄洪口高，为泄洪口平均流速。

[0087] 实施例2：

[0088] 1.1悬移质含沙量计算

[0089] 多沙河流中的泥沙输运大部分是以悬疑质运动的形式进行的。悬移运动的泥沙颗粒具有较细的粒径，可以跟随水流的紊动在水体中随机运动。在垂向上，泥沙颗粒的运动可以看做是两种运动的叠加，即重力驱动下的沉降运动和水流紊动驱动下的随机运动。当颗粒的数量很大时，将形成泥沙垂向运动的宏观动态平衡，此时泥沙浓度在垂向上有一个稳定的分布。

[0090] 1)悬移质含沙量的垂线分布

[0091] 泥沙的重力沉降使得含沙水流中沿垂线形成上清下浑的浓度分布。紊流中沿水深不同高度处各层水体之间存在水团的紊动交换，同时引起各水层间泥沙的交换，但上清下浑的浓度分布使得向上运动的水团所挟带的沙量大于乡下运动水团所挟带的沙量，所以紊动交换的结果是形成一个向上运动的泥沙通量qs1。另一方面，由于泥沙比水重，势必往下沉降并形成一个向下运动的泥沙净通量qs2。如果悬移质含沙量沿垂线出现稳定的时均泥沙浓度分布，则说明qs1与qs2达到了动平衡状态。此时，泥沙的紊动扩散过程是均匀、恒定的，二维扩散方程成为：

[0092]

[0093] 式(1)对y积分一次，令常数为零，得到

[0094]

[0095] 求解式(2)时的主要问题是必须首先知道εy沿垂线的分布。如假定εy为常数(意味着在垂线上紊动是均匀的)，则上式的解为

[0096]

[0097] (3)；

[0098] 其中，Sva为悬移质泥沙在距床面为a处的体积比含沙浓度。

[0099] 研究表明，泥沙的扩散系数εy不是常数而是空间位置的函数，但现有的理论还不能给出εy沿垂线的分布规律。最常用的方法是假定泥沙扩散系数εy与动量交换系数εm相等。剪切紊流中相邻流层间的剪切应力，主要是由于流体脉动导致的相邻流层间动量的交换所引起的，可以仿照分子粘性应力的表达方式给出其表达式：

[0100] 或

[0101] 由明渠恒定流中剪切力沿垂线的线性分布 和纵向流速沿垂线的对数型分布 (其中τ0＝γhJ，h为水深)，可以得到：

[0102]

[0103] 其中，k＝0.4为Karman常数。

[0104] 将式(4)代入式(2)得：

[0105]

[0106] 式(5)又可以写成：

[0107]

[0108] 对上式积分得：

[0109]

[0110] 令y＝a为参考点，该点泥沙浓度记为Sva，最后可得：

[0111]

[0112] 其中，指数Z的表达式为：

[0113]

[0114] Z又称为悬浮指标。式(8)中悬浮指标Z的大小决定了泥沙在垂线上分布的均匀程度。Z值越小，悬移质分布越均匀。

[0115] 2)水库悬移质含沙量垂线分布曲线

[0116] 本发明基于悬移质含沙量的垂线分布公式(8)，利用Matlab软件，开发出一套悬移质输沙率的交互式计算(SSDC,Suspended Sediment Discharge Calculation)软件，并利用该软件，结合实测含沙量的垂线分布数据，计算出水库悬移质含沙量的垂线分布曲线。本发明软件设计的主要思路为：输入实测含沙量垂线分布n个数据点，设置水深数据h，利用含沙量垂线分布函数(8)拟合上述n个实测数据点，获得分布函数中的参数值Sva,a和Z。采用基于最小二乘法的非线性曲线拟合算法实现对数据点的拟合。

[0117] 以某水库2011年主汛期坝前河段主流垂线含沙量的实测数据为例，如表1中所示。由于是主汛期，上游来流量较大，坝前水位低，水库表现为河道特征，水流含沙量明显增加，主流垂线表、中、底层含沙量有分层现象，近底层及底层含沙量较大。

[0118] 表1某水库坝前垂线含沙量实测值(2011年8月6日、7日测量)

[0119]

[0120] 确定含沙量垂线分布实测点个数n，水深数据h，并选择含沙量垂线分布拟合函数，在SSDC软件中的输入，然后依次输入n个实测数据点，实测含沙量数据沿垂线分布的拟合参数显示在软件的交互界面上，并作为共享数据存储在程序后台，以供其他程序模块调动。含沙量垂线分布拟合曲线如图1所示。此程序模块实现了从若干离散的测点数据得到连续的含沙量垂线分布曲线的功能。

[0121] 1.2坝前流速垂线分布

[0122] 本发明采用紊流断面流速分布公式计算水库坝前流速的垂线分布，如式(10)所示：

[0123]

[0124] 该式是紊流断面流速分布的对数公式。虽然它是根据壁面附近的条件推导出来的，但实验研究表明，该式适用于除粘性底层以外的整个过流断面。

[0125] 上面由于混合长度理论得到了紊流时均切应力表达式和流速对数分布规律。但这一理论的基本假设不够严谨，如认为流体质点经过混合长度后，才一次性与周围质点进行动量交换，然而由于流体是连续介质，质点在横向运动过程中不断地与周围质点进行动量交换。又如公式中的常数k,C都需要由实验确定。尽管如此，由于这一理论是从紊流的基本特征出发，同时推导简单，理论结果与实验比较一致，故至今仍是工程上得到广泛应用的紊流阻力理论。

[0126] 本研究采用公式(10)中的对数型流速分布公式，在SSDC软件中建立流速垂线分布计算模块，利用实测流速数据，拟合得出式(10)中的参数，并进而得到水库坝前流速垂线分布曲线。

[0127]

[0128] 某水库2011年主汛期坝前河段主流流速实测数据如表2所示。由于是主汛期，上游来流量较大，坝前水位低，水库表现为河道特征，水流流速均明显增加，坝前河段主流表层水流流速大于其他层。

[0129] 表2水库坝前垂线流速实测值(2011年8月6日、7日测量)

[0130]

[0131] 输入实测流速垂线分布n个数据点，利用流速垂线分布函数(11)拟合上述n个实测数据点，获得分布函数中的参数值U’和ks。采用基于最小二乘法的非线性曲线拟合算法实现对数据点的拟合。

[0132] 确定含流速垂线分布实测点个数n，并在SSDC中依次输入n个实测数据点，实测流速数据沿垂线分布的拟合参数显示在软件的交互界面上，并作为共享数据存储在程序后台，以供其他程序模块调动。流速垂线分布拟合曲线如图2所示。此程序模块实现了从若干离散的测点数据得到连续的流速垂线分布曲线的功能。

[0133] 1.3水库输沙率计算

[0134] 1)单宽输沙率

[0135] 根据前面介绍的理论方法，可以求出悬移质含沙量Sv及流速U沿垂线的分布。据此可计算得出单位时间内，河道过水断面中高程y处单位时间、单位断面面积上通过的悬移质沙量为USv，将其沿垂线积分即可得出全水深上的悬移质单宽输沙率。

[0136] 在之前的分析中，只得出了悬移质沿垂线的相对分布，当推求悬移质输沙率时，必须借助于其他途径得出参考高度a处的含沙量Sva，才能求出垂线上各点的含沙量。此外，求积分时首先还必须确定积分的上下限值。最简单的方法是从床面积分到水面处。但在已有的理论公式中，当y＝0时，流速和含沙量将分别趋于正、负无穷大。另一方面，从输沙的物理过程看，直接积分到床面也是不合理的，因为在床面附近运动的泥沙，其重量是由床面支持的，而不是由水流的紊动能量支持的，所以这一层运动的泥沙属于推移质的范畴。

[0137] 水深范围[h,h+H]内河道的单宽输沙率为该范围内的悬移质含沙量的垂线分布与紊流对数型流速分布之乘积在该深度范围内的积分。基于悬移质含沙量的垂线分布公式(8)和紊流对数型流速分布公式(11)计算单宽输沙率qsw，如式(12)所示：

[0138]

[0139] 2)水库输沙率计算

[0140] 水库主汛期上游来流量较大，坝前水位低，水库表现为河道特征，水流含沙量、流速均明显增加，但水流流速与含沙量的分布特征并不一致，坝前河段主流表层水流流速大于其他层，主流垂线表、中、底层含沙量有分层现象，近底层或底层含沙量较大，主汛期为排沙的主要时期，并主要依靠泄洪口排沙。

[0141] 以某水库为例，分别计算水库泄洪口全开状态下的单宽输沙率qsw、输沙率Qsw、泄洪口平均流速 泄洪口平均含沙率 泄洪流量QF，如式(13)～式(16)所示。

[0142] 输沙率计算：

[0143] Qsw＝n×L×qsw   (13)；

[0144] 泄洪口平均流速：

[0145]

[0146] 泄洪口平均含沙率：

[0147]

[0148] 泄洪流量：

[0149]

[0150] 以上各式中，n为泄洪口个数，h为泄洪口底部高程，L为泄洪口宽度，H为泄洪口高。

[0151] 该水库共有n＝23个泄洪口，泄洪口底部高程h＝90m，泄洪口宽L＝7m，泄洪口高H＝9m。

[0152] SSDC程序将计算该水库汛期泄洪口平均流速、泄洪口平均含沙率、泄洪流量、单宽输沙率及泄洪口总输沙率。程序中变量量纲统一，本算例使用国际标准单位制。

[0153] 1.4水库排沙发电多目标优化模型

[0154] 泥沙淤积与蓄水发电是水库运行中的一对非常突出的矛盾。若要减少淤积和有效排沙，需降低库水位并加大泄流，这会影响发电效益；若要提高水库的发电效益，需提高水库的蓄水位以及减少弃水，这又会导致泥沙淤积的发展。为协调水库运用过程中泥沙淤积与发电效益的矛盾关系，本研究以排沙效率最大和发电效益最大为基本目标，建立水库水沙电多目标优化调度模型，以合理安排水库的运行方案，从而提高枢纽的综合效益。

[0155] 1)目标函数

[0156] 目标函数一：系统发电量最大；

[0157] 以调度期内发电量最大作为一个目标，目的是尽量利用水电站的调节能力，加大发电平均水头同时减少弃水，尤其是充分利用梯级电站的水力联系和流域间的电力联系，实现在考虑系统用电负荷的情况下，最大程度的利用水能资源。目标函数可写为：

[0158]

[0159] 其中，E是年发电总量；是单机平均出力；Δt为计算时长；i为水电站的各机组编号；ni为各机组的台数；t代表时段长度。

[0160] 目标函数二：水库输沙量最大；

[0161] 以水库输沙量最大作为另一个目标，目的在于加大泄洪时的输沙。水库的入库泥沙主要来自于汛期，因此根据水库的输沙率计算公式，增大水库汛期的排沙量。目标函数可写为：

[0162]

[0163] 其中，S是水库的输沙量，单位为kg；QF为泄洪流量，由式(16)计算可得，单位为m3/s， 为泄洪口平均含沙率，由式(15)计算而得，单位为kg/m3。

[0164] 2)决策变量

[0165] 由前分析可知，泄流量既为水电站出力的影响因素，又为水库输沙量的决定因素。因此，选取水库的泄流量为决策变量，同时，将发电引用流量和泄洪流量区分开来。在非汛期的决策变量为水库的泄流量，即发电引用流量Q；在汛期的决策变量为水库的发电引用流量Q和泄洪流量QF。

[0166] 以某水库为例，由于该水库的汛期为6月至9月共4个月，取月为时段步长，决策变量u可表示如下：

[0167] u＝{Q1,Q2,…,Q12,QF1,QF2,QF3,QF4}   (19)；

[0168] 即应用遗传算法时，每个染色体有16个基因。

[0169] 3)约束条件

[0170] 水库系统运行的约束条件分为以下几种：

[0171] 库水位约束

[0172] 为了保证大坝安全和保护下游地区，水库必须在汛期预留出一定的防洪库容，将水位降到汛限水位以下。在非汛期也必须遵循调度规程。水库的安全运行通过对库水位的控制实现：

[0173]

[0174] 其中， 和 分别为第t时段允许的最低和最高水位。

[0175] 水库下泄流量约束

[0176] 出库流量的下限是为了保证下游通航和满足生态运作的需要，其上限取水库最大允许泄流量。因此，优化需要在大范围内寻优。这会产生大量不可行解，导致水库放空或大坝漫顶的极端结果，遗传算法的种群内充斥着无效个体。

[0177] 为了改进遗传算法的效率，根据入库流量和实际调度过程来设定出库流量的上下限，其变化范围为：

[0178]

[0179] 其中， 和 是当月入库流量的最小和最大值； 和 是当月的实际出库流量的最小和最大值。

[0180] 流量日变幅约束

[0181] 作为中国最大的河流，长江承担着重要的航运任务，这不仅要求河道的水位达到一定的高度，也要求水位变化不能太剧烈。

[0182] 所有的约束条件通过两种方式实现：要么在产生初始种群和新种群时自动满足，要么在计算后检验是否被满足，不可行解被剔除，不参与进化。

[0183] 4)状态方程

[0184] 状态方程是指计算过程中水力电力条件和联系，主要有以下几类：

[0185] 水库水量平衡方程：

[0186] Vt＝Vt-1+It×Δt-(Qt+QFt)×Δt   (22)；

[0187] 其中，V为库容，I为上游入库流量；Q为水电站发电引用流量；QF为泄洪流量；(Q+QF)为水库总泄流量，包括发电引用流量和泄洪流量。

[0188] 水位库容曲线

[0189] 采用分段线性插值函数f(x)来表达水库的水位库容关系：

[0190] Vt＝f(Lt)

[0191] Lt＝f-1(Vt)   (23)；

[0192] 其中，L为水库库水位。

[0193] 水库尾水位-出库流量关系曲线

[0194] 对水库尾水位-出库流量关系曲线同样采用分段线性函数插值函数g(x)表达：

[0195] Tt＝g(Qt+QFt)   (24)；

[0196] 其中，T为下游尾水水位；(Q+QF)为水库总泄流量。

[0197] 机组出力特性曲线

[0198] 电站机组出力曲线采用二元函数全微分的方式进行线性化，当水头一定时，机组出力可表达为单机引用流量的线性函数：

[0199] N＝C1(H)×q+C2(H)   (25)；

[0200] 其中，C1(H)和C2(H)是某一水头下的线性系数。

[0201] 流速垂线分布曲线

[0202] 本研究采用对数型流速分布公式，利用实测流速数据，拟合得到水库坝前流速垂线分布曲线，流速U可表示为相对深度dr的函数：

[0203] U＝h(dr)   (26)；

[0204] 悬移质含沙量垂线分布曲线

[0205] 本研究采用公式(8)中的基于扩散理论的悬移质含沙量垂线分布公式，利用实测流速数据，拟合得到水库悬移质含沙量垂线分布曲线，泥沙浓度Sv可表示为相对深度dr的函数：

[0206] Sv＝k(dr)   (27)；

[0207] 1.5水库排沙发电多目标优化流程

[0208] 1)染色体编码

[0209] 随机产生n个染色体，即产生n个水库下泄流量Q的序列，组成初始种群P0(种群规模＝n)；以月为单位，每个染色体上有12个发电引用流量元素和4个泄洪流量元素共16个基因。其中，第i个染色体的编码如式(28)所示：

[0210] Qi＝{QNi[1]，QNi[2]，...，QNi[t]，...，QNi[12]，

[0211] QFi[1+5]，QFi[2+5]，QFi[3+5]，QFi[4+5]}   (28)；

[0212] 其中，QN为发电引用流量，用于计算水电站出力：

[0213] N(N＝K·QN·H)；QF为泄洪流量；

[0214] 用于计算水库排沙量：

[0215] S(S＝QF·QS)，QTi(QTi[t]＝QNi[t]+QFi[t])为总泄流量，用于水量平衡方程，计算下一时刻的库容水位等参数。

[0216] 需要说明的是，由于一般来说，只有在汛期6-9月份的四个月会发生水库泄洪，因此，泄洪流量QF只编码了四个变量，且与总泄流量QT和发电引用流量QN对应时，需要考虑时间5个月的时间差。

[0217] 为了减少搜索空间，提高优化效率，采用实际的泄流序列为基准，在一定的范围内上下波动一定流量，构成搜索空间。即首先获取了当年黄陵庙站检测得到的月均泄流量和由水电站得到的月均总发电引用流量(汛期二者之差为泄洪流量)，然后设定日均流量波动范围值Qdr，在[-Qdr，Qdr]范围内随机产生波动值，此波动值与检测获得的日均流量值相加即为程序生成的日均流量值。按以上算法产生种群规模为n的初始种群，存储在16×n的数组中。

[0218] 2)目标函数值计算

[0219] 计算目标函数值，包括水电站年总出力和水库年总排沙量。

[0220] 年总出力计算

[0221] 编写求解第一个目标函数(水电站年总出力)的函数Calculate_Obj_E。算法如下：

[0222] ●计算库容(利用水量平衡方程)：

[0223] Vi[t]＝Vi[t-1]+Li[t]×0.0264-QNi[t]×0.0264-QFi|[t]×0.0264   (29)；

[0224] 其中0.0264为单位换算系数0.0264＝30.5×24×3600/108，保证库容的单位为亿立方米。初始库容值设置为Vi[0]＝393亿立方米。

[0225] ●计算上游水位

[0226] ●编写内插值函数Li[t]＝f-1(Vi[t])，水位-库容关系数据通过内部插值方法确定任意给定库容值下的水位数据。

[0227] ●计算下游水位Ti[t](尾水位Tailwater Level)；

[0228] 下游水位与下泄流量有关，Ti[t]＝g(QNi[t]+QF|i[t])；

[0229] ●计算月均出力,先令K[t]＝const；

[0230] Ni[t]＝K[t]×QNi[t]×(Li[t]-Ti[t])   (30)；

[0231] ●计算水电站年总出力，为第i个染色体对应的第一个目标函数值，即适应度fitness1.

[0232] Obj1i＝Ni＝Ni[1]+Ni[2]+…+Ni[t]+…+Ni[12]   (31)；

[0233] 水库年总排沙量Si

[0234] 编写求解第二个目标函数(水库年总排沙量)的函数Calculate_Obj_S。算法如下：

[0235] ●确定含沙量

[0236] 汛期泄洪口含沙量由垂线含沙量分布曲线在泄洪口高度范围积分得到。汛期过机含沙量和非汛期过机含沙量由实测数据得到。汛期过机含沙量用汛期的月平均值来表示。

[0237] ●计算月排沙量

[0238] 计算第i个染色体上第t个月的排沙量：

[0239] Si[t]＝QFi[t]×QSi[t]×0.0264   (32)；

[0240] 其中0.0264为单位换算系数。

[0241] ●计算年排沙量

[0242] 计算第i个染色体对应的水库年总排沙量Si，为第i个染色体对应的第二个目标函数值，即适应度fitness2。

[0243] Obj2i＝Si＝Si[1]+Si[2]+…+Si[t]+…+Si[12]   (33)；

[0244] 3)多目标遗传优化算法优化

[0245] 在n个染色体中选取较优的m个值进行交叉变异计算，生成新种群P1：

[0246] ●确定non-dominated solutions，即非支配解：

[0247] 遍历种群中任意两对染色体，比较两个染色体i和j的适应度函数值Obj1和Obj2,如果(Obj1i>Obj1j&&Obj2i>Obj2j)，删除染色体j。

[0248] ●将当前m个non-dominated solutions进入交叉变异

[0249] 若m＝0，重新生成初始种群；若m>0，进入交叉变异。

[0250] 交叉：设置交叉概率为0.75。若交叉发生，则从m个non-dominated solutions中任意选取两个染色体，随机产生一个[1,16]之间的整数，大于此数的基因位置发生交换。

[0251] 变异：设置变异概率为0.02。若变异发生，则对于上述染色体，随机产生一个[1,16]之间的整数，在该位置上的基因突变为m个染色体中该位置上最大值和最小值之间的任意数。

[0252] ●交叉变异后新生成n个染色体，与上一次迭代得到的非支配解一起构成新种群P1，种群规模为(m+n)。

[0253] ●计算新种群中的染色体对应的目标函数值，进入下一次迭代。

[0254] 实施例3；

[0255] 以某水库为例，采用本发明方法建立的水库排沙发电多目标优化方法进行实例计算，参数设置如下：初始种群规模为：100，迭代步数：100步，汛期泄洪口含沙量：0.066738kg/m3，汛期过机含沙量：0.026722kg/m3，非汛期过机含沙量：0.009015kg/m3。交叉概率：0.75，变异概率：0.02，日均流量变化幅值Qdr从50m3/s变化到5000m3/s。

[0256] 以下图3至图6中，散点表示初始种群对应的目标函数；从左下到右上的数据点集合链分别表示第N次迭代后非支配解对应的目标函数，即Prato front曲线；图中分别显示了1次迭代、5次迭代、10次迭代、25次迭代、100次迭代后的非支配解对应的目标函数，由此目标函数，决策者可以方便地选择系统发电量最大和水库输沙量最大的多目标优化结果，例如从图中最靠近右上角的数据点中进行选择。

[0257] 图3至图6中，纵轴表示年排沙量，单位为亿吨/每年，横轴代表年发电量，单位为万千瓦。图3至图6中日均流量变化幅值Qdr为依次为50m3/s、150m3/s、1000m3/s、5000m3/s。