CN111210141A - 基于约束机制粒子群算法的水库库容曲线修正方法 - Google Patents
基于约束机制粒子群算法的水库库容曲线修正方法 Download PDFInfo
- Publication number
- CN111210141A CN111210141A CN202010006199.4A CN202010006199A CN111210141A CN 111210141 A CN111210141 A CN 111210141A CN 202010006199 A CN202010006199 A CN 202010006199A CN 111210141 A CN111210141 A CN 111210141A
- Authority
- CN
- China
- Prior art keywords
- reservoir
- area
- formula
- water level
- function
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
- G06Q10/063—Operations research, analysis or management
- G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
- G06Q50/06—Electricity, gas or water supply
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Abstract
本发明公开了一种基于约束机制粒子群算法的水库库容曲线修正方法,包括下述步骤:1,构建决策变量内部约束机制;2,水库运行资料预处理;3,目标函数设计;4,基于约束机制的改进粒子群算法求解目标函数;步骤5,输出条件判别。本发明提出以不同水位对应的水面面积作为决策变量,并作用于符合实际库区地形变化规律的约束机制,解决了常规水量平衡法求解结果中可能出现的高水位水面面积小于低水位水面面积这种违背常理的问题,通过运行资料预处理,分离了部分水位和部分库容,阻断了误差累计和传播途径,通过筛选水库蓄水和消落对称过程的资料作为输入,有效降低了流量资料中系统误差对库容曲线计算结果的影响,提高了库容曲线求解精度。
Description
技术领域
本发明涉及水利工程中的已建水库库容曲线修正方法,尤其是涉及基于约束机制粒子群算法的水库库容曲线修正方法。
背景技术
水库库容曲线是水库工程规划设计和运行管理的必备参数和基础成果之一,其精度对水库调控方式及效益产生重大影响。目前,我国大约有9.8万余座已建成水库工程,仅有少数大型的、特别重要的、泥沙淤积问题严重的水库会经常性的对库区地形进行测量,但由于测量代价较大,大多数水库库容曲线极少进行复核测量,若在运行管理期间库区地形形态发生显著性变化,对水库运行调度决策将产生重大影响。因此,如何采用成本较低的方式复测水库库容曲线具有广阔的现实需求。通常获取水库库容曲线有两种途径:一是直接测量法:对库区陆地和水下地形进行全面、高精度的测量;二是水量平衡法:基于水库运行资料,采用水量平衡方程对库容要素进行推算。
直接测量法又分为地形图法和横断面法。地形图法通过全库区周边及水下范围的测量,“自下而上”进行累计计算,该方法所得成果可靠、精度高,但外业工作量巨大,对测量设备要求高,需投入大量人力、物力和财力,作业周期较长,因此该方法适用于库区范围较小、库区冲淤变化不显著的水库或库容影响重大的水库。
水量平衡法根据水库不同水位高程下的入库、出库流量、降雨、蒸发、渗漏等观测资料,依据水量平衡方程计算对应库容差,综合形成库容曲线。该方法几乎不需开展外业工作,对库区范围、形态也没有限制,仅对运行资料的精度、时段、完整性等有一定要求,因而适用范围相对较广,成本相对较低。但是,水量平衡法也存在一定的局限性和问题:首先,由于水文测验精度有限,库周汇流难以精确测算,入库流量计算可能存在一定的误差;其次,水库蓄泄运用时,库岸坍塌、库区泥沙冲淤将引起不同高程下的水面面积发生变化,影响了库面降雨直接入库水量及库面蒸发水量的测算,进而增大入库、出库水量的计算误差;再次,由于库容值是由库容差值由低水位向高水位累加求和,导致库容误差被不断累计和传播;最后,水量平衡法仅得到库容值,尽管通过试算法可反推对应的水面面积,但会出现高水位水面面积小于低水位水面面积这种违背常理的不协调现象。因此,如何克服上述水量平衡法存在的问题,提高水库库容曲线精度是一个亟需解决的问题。
发明内容
本发明目的在于提供一种基于约束机制粒子群算法的水库库容曲线修正方法,以不同水位对应的水面面积作为决策变量,通过构建符合变量内部相邻成员间变化规律的约束机制,运用改进的粒子群算法,求解与运行资料偏差最小(即离差平方和最小)的库容曲线。
为实现上述目的,本发明采取下述技术方案:
本发明所述基于约束机制粒子群算法的水库库容曲线修正方法,包括下述步骤:
步骤1,构建决策变量内部约束机制:通过分析水库库区地形的一般变化规律,根据水库水位~面积函数的微分特性及实际应用时水库水位~面积~库容关系为离散形式,以不同水位对应的水面面积作为决策变量,构建符合变量内部相邻成员间变化规律的约束机制;
步骤2,水库运行资料预处理:为便于分析和比较求解成果与运行资料的匹配程度,根据目标水库的运行资料、气象资料、库区地质资料以及近期库容曲线,由水量平衡原理,逐时段计算初、末水位及其对应的时段库容差k=1,2,……,P,并以此系列作为基准库容差;
步骤3,目标函数设计:针对某一备选库容曲线,在曲线上逐时段查得运行资料中的初、末水位及其对应的时段库容差值dVk,以此系列作为备选方案库容差;以时段备选方案库容差与基准库容差的离差平方和最小为所述目标函数,即:
步骤4,基于约束机制的改进粒子群算法求解目标函数,步骤如下:
子步骤4.1,设定算法基本参数,初始粒子数量M、粒子维度N、更新迭代周期上限S,学习因子c1、c2均取2,速度限制常数vmax取值为1/(2N),惯性权重w采用线性递减自适应调整策略;
子步骤4.2,建立粒子与面积函数决策变量的映射关系;
子步骤4.3,基于步骤1中相关的约束机制,随机生成满足要求的初始粒子群体;
子步骤4.5,根据改进粒子群算法,粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置;
子步骤4.6,重复子步骤4.4、4.5,直至粒子群更新迭代周期上限S时,输出群体最好位置,经换算后即为求解结果;
步骤5,输出条件判别:设定允许距离阈值为δ*,计算步骤4中求解的库容曲线结果与步骤2中水库运行资料中库容曲线的欧氏距离δ,若δ≤δ*,则输出步骤4计算的库容曲线结果;反之,用步骤4的库容曲线替代步骤2中近期库容曲线,并重复步骤2、3、4,直至满足输出条件为止。
所述步骤1中,构建符合变量内部相邻成员间变化规律的约束机制步骤为:
子步骤1.1,通过分析水库库区地形的变化规律,提出水库水位~面积函数的微分特性为:
设水库水位变量Z,对应水面面积函数A=f(Z),根据库区地形的一般变化规律可知,随着水库水位的抬高,相应水面面积越大,即面积函数的一阶导数大于零,有:
A′=f′(Z)=dA/dZ>0 式(1)
另外,根据水位越高水面面积增长率也越来越大,即面积函数的二阶导数也大于零,因此有:
A″=f″(Z)=d2A/dZ2>0 式(2)
式2可以看出,面积函数为严格凹函数,由其性质可知,对于x1<x2,存在0<α<1,有:
f(αx1+(1-α)x2)<αf(x1)+(1-α)f(x2) 式(3)
子步骤1.2,基于子步骤1.1的面积函数微分特性,等差水位相邻离散点距的相互制约关系为:
设库容曲线中离散水位值为Zj,(j=0,1,2,……,N),从小到大排列的序列,对应水面面积系列为Aj=f(Zj),由式(1)有:
(Aj-Aj-1)/(Zj-Zj-1)>0,
因Zj>Zj-1,有Aj-Aj-1>0,即:
0≤Aj-1<Aj(j=1,2,……,N) 式(4)
由式(2)有:
对等差水位序列,简化为:
子步骤1.3,针对某些库区的水位~面积曲线存在面积函数局部区间为非凹函数时,约束条件变更判别式及非凹函数情况下的制约关系为:
对于遵循凹函数规律的水位~面积曲线而言,式(5)是严格成立的,有(2×Aj-1-Aj-2)<(Aj+1+Aj-1)/2,整理得:
Aj+1>(3×Aj-1-2×Aj-2)(j=2,3,......,N-1) 式(6)
因此,式(6)是判别Aj有无满足凹函数性质解的条件;对于部分特殊位置的地形,存在不满足凹函数的情况,即式(5)无解或式(6)不成立,所以式(6)也是约束条件变更的判别式,若式(6)不成立,说明在该库区高程范围对应的面积增长率非正值,此时的面积函数应变更为非凹函数约束,即:
(Aj+1+Aj-1)/2≤Aj<Aj+1(j=1,2,......,N-1) 式(7)
子步骤1.4,初始化水位~面积曲线时,由于Aj+1未知,无法利用式(5)~式(7)来约束Aj,根据式(3),其离散表达形式为:
Aj<A0+(AN-A0)×(Zj-Z0)/(ZN-Z0) 式(8)
对等差水位序列,进一步简化为:
Aj<A0+(AN-A0)×j/N(j=1,2,……,N-1) 式(9)
由此,在初始化水位~面积曲线时,利用式(8)或式(9)来代替式5中的约束条件上边界。
所述步骤2中,所述目标水库的运行资料包括库水位、入库流量、出库流量;所述气象资料包括库区降雨量、蒸发量。
本发明优点体现在以下方面:
1、本发明提出以不同水位对应的水面面积作为决策变量,并作用于符合实际库区地形变化规律的约束机制,计算结果更为合理,解决了常规水量平衡法求解结果中可能出现的高水位水面面积小于低水位水面面积这种违背常理的问题。
2、本发明提出以不同水位对应的水面面积作为决策变量,解决了常规水量平衡法计算过程中无法精确计算时段蒸发量和渗漏量的问题,在一定程度上有效保证了库容曲线求解精度。
3、本发明通过运行资料的预处理,分离了部分水位和部分库容,阻断了误差累计和传播途径,将误差限定到一定影响范围内;通过筛选水库蓄水和消落这种对称过程的资料作为输入,有效降低了流量资料中系统误差对库容曲线计算结果的影响,进一步提高了库容曲线求解精度。
附图说明
图1是本发明方法的流程框图。
具体实施方式
下面结合附图对本发明的实施例作详细说明,本实施例在以本发明技术方案为前提下进行实施,给出了详细的实施方式和具体的操作过程,但本发明的保护范围不限于下述实施例。
如图1所示,本发明所述基于约束机制粒子群算法的水库库容曲线修正方法,包括以下步骤:
步骤1,构建决策变量约束机制
本发明以水库不同水位对应的水面面积作为决策变量,因此满足决策变量内部相邻成员间的约束条件是确保求解结果符合库容曲线自然变化特性的基础和前提。本发明通过分析水库库区地形的一般变化规律,总结提出水库水位~面积函数的微分特性,考虑实际应用时水库水位~面积~库容关系通常为离散值,提出等差水位条件下相邻离散点距的相互制约关系;具体包括以下子步骤:
子步骤1.1,通过分析水库库区地形的一般变化规律,总结提出水库水位~面积函数的微分特性;
设水库水位变量Z,对应水面面积函数A=f(Z),根据库区地形的一般变化规律可知,随着水库水位的抬高,相应水面面积越大,即面积函数的一阶导数大于零,有:
A′=f′(Z)=dA/dZ>0 式(1)
另外,地形越向高处通常也越开阔,也就是说水位越高,水面面积增长率也越来越大,即面积函数的二阶导数也大于零,因此有:
A″=f″(Z)=d2A/dZ2>0 式(2)
式2可以看出,面积函数为严格凹函数,由其性质可知,对于x1<x2,存在0<α<1,有:
f(αx1+(1-α)x2)<αf(x1)+(1-α)f(x2) 式(3)
子步骤1.2,基于子步骤1.1的面积函数微分特性,提出一般情况下等差水位相邻离散点距的相互制约关系;
设库容曲线中离散水位值为Zj(j=0,1,2,……,N),为从小到大排列的序列,对应水面面积系列为Aj=f(Zj),由式(1)有:
(Aj-Aj-1)/(Zj-Zj-1)>0
因Zj>Zj-1,有Aj-Aj-1>0,即:
0≤Aj-1<Aj(j=1,2,……,N) 式(4)
由式(2)有:
对等差水位序列,可简化为:
子步骤1.3,针对某些库区的水位~面积曲线可能存在局部“反常”情况,即面积函数局部区间为非凹函数时,提出约束条件变更判别式及非凹函数情况下的制约关系;
对于遵循凹函数规律的水位~面积曲线而言,式(5)是严格成立的,有(2×Aj-1-Aj-2)<(Aj+1+Aj-1)/2,整理得:
Aj+1>(3×Aj-1-2×Aj-2)(j=2,3,......,N-1) 式(6)
因此,式(6)是判别Aj有无满足凹函数性质解的条件;对于部分特殊位置的地形,存在不满足凹函数的情况,即式(5)无解或式(6)不成立,所以式(6)也是约束条件变更的判别式,若式(6)不成立,说明在该库区高程范围对应的面积增长率非正值,此时的面积函数应变更为非凹函数约束,即:
(Aj+1+Aj-1)/2≤Aj<Aj+1(j=1,2,......,N-1) 式(7)
子步骤1.4,初始化水位~面积曲线时,由于Aj+1未知,无法利用式(5)~式(7)来约束Aj;根据式(3),其离散表达形式为:
Aj<A0+(AN-A0)×(Zj-Z0)/(ZN-Z0) 式(8)
对等差水位序列,可进一步简化为:
Aj<A0+(AN-A0)×j/N(j=1,2,......,N-1) 式(9)
由此,在初始化水位~面积曲线时,可利用式(8)或式(9)来代替式5中的约束条件上边界;
步骤2,水库运行资料预处理
为便于分析和比较求解成果与运行资料的匹配程度,根据目标水库的运行资料(库水位、入库流量、出库流量)、气象资料(库区降雨量、蒸发量)、库区地质资料及近期的库容曲线,由水量平衡原理,逐时段计算初、末水位及其对应的时段库容差,如式(10)所示,并以此系列作为基准库容差;
其中:k为当前计算时段标识,k=1,2,......,P,P为总时段数;
Qs,k、Qe,k分别为第k时段初、末的入库流量;
qs,k、qe,k分别为第k时段初、末的出库流量;
Δtk为第k时段的时段长度,时段划分时应保证时段内入库、出库流量呈线性变化;
Wpre,k、Wevap,k、Wseep,k分别为第k时段库区的降水量、蒸发量及渗漏量;
步骤3,目标函数设计
首先,针对某一备选水位~面积曲线,根据椎体体积计算公式,计算与之对应的库容曲线,计算公式如下:
其中:c、j均为库容曲线离散点序列标识,c,j=1,2,......,D,D为离散点总数;
Zc、Ac、Zj分别为第c个序列的水位值、面积值和第j个序列的库容值;
dVk的计算公式如下:
dVk=Ve,k-Vs,k 式(12)
其中:dVk为从备选库容曲线上查得的第k时段库蓄水量变化,即时段备选方案库容差;
Vs,k、Ve,k分别为从式(11)计算的备选库容曲线上查得的第k时段初、末的水库蓄水量;
f=Min(Φ(dV)) 式(13)
其中:f为目标函数;
P为运行资料的时段数量;
Min()为最小值函数;
步骤4,基于约束机制改进粒子群算法的目标函数求解
具体包括以下子步骤:
子步骤4.1,设定算法基本参数,初始粒子数量M、粒子维度N、更新迭代周期上限S的具体取值应视求解问题而定,学习因子c1、c2均取2,速度限制常数vmax取值不宜过大,建议取值vmax=1/(2N);
惯性权重w采用典型的线性递减自适应调整策略,计算公式如下:
w=wmax-(wmax-wmin)×s/S 式(15)
其中:wmax、wmin分别为设定的最大、最小惯性权重;
s为当前迭代周期;
子步骤4.2,建立粒子与决策变量间的映射关系;
设第i个粒子xi=(xi,0,xi,1,…,xi,j,…,xi,N)与决策变量Ai=(Ai,0,Ai,1,…,Ai,j,…,Ai,N)的映射关系如下:
Ai,j=A0+(AN-A0)×xi,j 式(16)
其中:i=1,2,......,M,xi,j∈[0,1];
式(16)反函数为:
xi,j=(Ai,j-A0)/(AN-A0) 式(17)
子步骤4.3,基于步骤1的相关约束机制(式(5)及式(8)或式(9)),随机生成满足约束要求的初始粒子群体;
初始化的粒子内各维度变量间应满足如下关系:
因此,建议通过如下方式生成第i个粒子第j个维度的初始位置和速度:
其中:Sgn为符号函数,取值为{-1,0,1};
Rnd为随机函数,且Rnd∈[0,1];
对每个粒子,首先根据式(16)和式(11)转换成对应的备选库容曲线,再由式(12)、(14)及其他数据计算对应的以备选方案库容差,以备选方案库容差与基准库容差的离差平方和作为该粒子的适应度值。因此,迭代周期s时各粒子自身最好位置对应的适应度为:
其中:Min()为最小值函数;。
子步骤4.5,根据改进粒子群算法粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置;
改进粒子群算法粒子的速度和位置更新公式如下:
其中:r1、r2为介于[0,1]区间的随机数;
根据式4、式5和式16,更新后的粒子各相邻维度位置之间应满足如下两个公式条件:
根据式6和式16,更新后的粒子当前维度位置存在满足凹函数解区间的判别公式如下:
式(23)用于确保水位~面积曲线的合理性,须强制性满足;
式(24)用于确保水位~面积曲线符合凹函数要求,存在解空间时应强制性满足;
式(25)是粒子当前维度是否存在凹函数解空间的判定条件;
因此,粒子当前维度位置约束过程如下:
式(23)成立时,判断式(24)是否成立:
若式(24)成立,说明粒子当前维度位置满足要求,可转入下一步骤;
若式(24)不成立,需进一步判断式(25)是否成立:
若式(25)成立,说明粒子当前维度存在凹函数解空间,当j=1或j=N时,采用下式重新生成粒子当前维度的位置和速度:
当j=2,3,......,N-1时,采用下式重新生成粒子当前维度的位置和速度:
若式(25)不成立,说明粒子当前维度不存在凹函数解空间,此时应变更粒子维度位置的约束机制,当j=1或j=N时,采用式(26)重新生成粒子当前维度的位置和速度;当j=2,3,......,N-1时,采用下式重新生成粒子当前维度的位置和速度:
式(23)不成立时,判断式(25)是否成立:
若式(25)成立,说明粒子当前维度存在凹函数解空间,采用式(26)或式(27)重新生成粒子当前维度的位置和速度;
若式(25)不成立,说明粒子当前维度不存在凹函数解空间,此时应变更粒子维度位置的约束机制,采用式(26)或式(28)重新生成粒子当前维度的位置和速度;
子步骤4.6,重复子步骤4.4、4.5,直至粒子群更新迭代周期至上限S时,群体最好位置g(S)经式(16)和式(11)转换后即为求解结果;
步骤5,输出条件判别
设定允许距离阈值δ*,计算步骤4中求解的库容曲线结果与步骤2水库运行资料中的库容曲线的欧氏距离δ,计算公式如下:
若δ≤δ*,则输出步骤4库容曲线结果;反之,用步骤4的库容曲线替代步骤2中近期库容曲线,并重复步骤2、3、4,直至满足输出条件。
本发明方法各步骤说明如下:
步骤1说明:
1)水库库容曲线的离散格式
水库库容曲线在理论上应是连续函数,即A=f(Z)和V=g(Z)
在实际应用当中,由于这种函数关系难以精确获取,因此通常情况下,水库库容曲线均采用离散格式来标识,如表1所示。
表1某水库库容曲线
变量名称 | Z | A | V |
序号 | 水位(m) | 面积(km<sup>2</sup>) | 库容(百万m<sup>3</sup>) |
0 | 110 | 0.0 | 0 |
1 | 120 | 3.3 | 11 |
2 | 130 | 7.3 | 62 |
3 | 140 | 14.3 | 169 |
4 | 150 | 25.4 | 365 |
表中,Z0=110,Z1=120,Z2=130,……;
A0=0.0,A1=3.3,A2=7.3,……
V0=0,V1=11,V2=62,……
这种离散形式也满足A1=f(Z1),V1=g(Z1),A2=f(Z2),V2=g(Z2),……
2)水库库容曲线具备性质的依据
水库库容曲线所具备的函数性质,主要依据天然河谷的地形特性:
A、水都是从高处向低处流,河道都是正比降;
B、由于风化作用和山体稳定的要求,随着高度的增加,多数山体的自然坡度会先增加再变小。
在上述两方面地形特性的作用下,反映到库容曲线形式上是水位~面积函数f(Z)的1阶导数大于零,2阶导数决大多数情况下大于零。
由此根据微分函数的导数特性,推导出离散形式下的面积变量Aj(j=0,1,2,……,N)的相互制约关系,包括式(4)、(5)和式(9)。
例如上述库容曲线中,A2>A1,即7.3>3.3。
以及(2×A1-A0)<A2<(A3+A1)/2,即(2×3.3-0)<7.3<(14.3+3.3)/2。
以上均用于说明对面积离散点Aj,与其邻近点之间存在本专利提出的这种约束关系,为建模寻优提供了理论和物理基础。
步骤2说明:
本步骤主要为资料预处理,核心为水量平衡原理,利用该原理求解分段库容,得到(时段初水位~时段末水位~库容差)这样一系列数组;不直接累加,可以避免误差的累积。
需要说明的是这里的库容差是利用水量平衡原理计算出来的,并非从库容曲线查出来的。
步骤3说明:
本步骤主要是构建目标函数,若已知水库库容曲线,根据步骤2中的时段初、末水位,可以直接在曲线中查得库容差值,如果所有测量数据没有误差,那么步骤2与步骤3这2个库容差值应是等值的。
但由于测量误差的存在,2个库容差值不相等,在此以离差平方和最小为目标函数,认为哪个库容曲线代入目标函数的离差平方和最小,即为目标库容曲线。
步骤4说明:
本步骤主要为目标函数求解;
子步骤4.1:设定算法基本参数
举例说明设置如下:粒子数量M=500,更新迭代次数K=40,粒子维度N=7,wmax=1.3,wmin=0.35,vmax=1/14,最低水位Z0=600m,对应面积A0=0km2,最高水位Z7=628m,假定对应面积上限为A7=50km2。
Z1~Z7分别为604、608、612、616、620、624、628,等4m间距。
求解目标为最优的一组A1~A7。
子步骤4.2:建立粒子与决策变量间的映射关系,如说明公式所示。
子步骤4.3:初始化粒子
例如随机生成的初始化粒子位置:
……
初始化粒子速度为0。
子步骤4.4:逐个粒子i计算相应的适应度
根据式(16)和式(11),转换为对应的备选水位~面积曲线:
……
由式(11)转换成对应库容曲线后,再根据(12)、(14)计算对应的目标函数值,作为粒子的适应度值。
子步骤4.5:根据标准PSO算法粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置。
根据式(22),更新粒子的位置和速度,
对于速度超限制,调整至边界上。
……
其余粒子更新和约束调整方式类似,不再赘述。
子步骤4.6:重复步骤(4)、(5),得到:
第1代群体最好位置:
g(1)=(0,0.017,0.036,0.060,0.103,0.224,0.428,0.931),
第10代群体最好位置:
g(10)=(0,0.013,0.047,0.090,0.140,0.203,0.353,0.614),
第16代群体最好位置:
g(16)=(0,0.029,0.059,0.097,0.135,0.176,0.356,0.601),
继续更新迭代群体最好位置不变,说明算法的收敛速度很快。
得到的库容曲线如下:
水位(m) | 面积(km<sup>2</sup>) | 库容(百万m<sup>3</sup>) |
600 | 0 | 0 |
604 | 1.45 | 1.93 |
608 | 2.95 | 10.56 |
612 | 4.85 | 26.00 |
616 | 6.75 | 49.10 |
620 | 8.8 | 80.11 |
624 | 17.8 | 132.26 |
628 | 30.05 | 226.90 |
步骤5说明:
分析步骤4得到的库容曲线与目标库容曲线(上一次计算结果)的欧式距离:
若欧氏距离在设代的阈值范围内,说明此次计算的库容曲线满足要求,可直接输出作为结果;
若在阈值范围之外,可用本次计算得到的库容曲线代替目标库容曲线,重复步骤2、3、4,直至满足输出条件为止。
Claims (3)
1.一种基于约束机制粒子群算法的水库库容曲线修正方法,其特征在于:包括下述步骤:
步骤1,构建决策变量内部约束机制:通过分析水库库区地形的一般变化规律,根据水库水位~面积函数的微分特性及实际应用时水库水位~面积~库容关系为离散形式,以不同水位对应的水面面积作为决策变量,构建符合变量内部相邻成员间变化规律的约束机制;
步骤2,水库运行资料预处理:为便于分析和比较求解成果与运行资料的匹配程度,根据目标水库的运行资料、气象资料、库区地质资料以及近期库容曲线,由水量平衡原理,逐时段计算初、末水位及其对应的时段库容差 并以此系列作为基准库容差;
步骤3,目标函数设计:针对某一备选库容曲线,在曲线上逐时段查得运行资料中的初、末水位及其对应的时段库容差值dVk,以此系列作为备选方案库容差;以时段备选方案库容差与基准库容差的离差平方和最小为所述目标函数,即:
步骤4,基于约束机制的改进粒子群算法求解目标函数,步骤如下:
子步骤4.1,设定算法基本参数,初始粒子数量M、粒子维度N、更新迭代周期上限S,学习因子c1、c2均取2,速度限制常数vmax取值为1/(2N),惯性权重w采用线性递减自适应调整策略;
子步骤4.2,建立粒子与面积函数决策变量的映射关系;
子步骤4.3,基于步骤1中相关的约束机制,随机生成满足要求的初始粒子群体;
子步骤4.5,根据改进粒子群算法,粒子更新方式及步骤1中提出的变量约束机制更新粒子速度和位置;
子步骤4.6,重复子步骤4.4、4.5,直至粒子群更新迭代周期上限S时,输出群体最好位置,经换算后即为求解结果;
步骤5,输出条件判别:设定允许距离阈值为δ*,计算步骤4中求解的库容曲线结果与步骤2中水库运行资料中库容曲线的欧氏距离δ,若δ≤δ*,则输出步骤4计算的库容曲线结果;反之,用步骤4的库容曲线替代步骤2中近期库容曲线,并重复步骤2、3、4,直至满足输出条件为止。
2.根据权利要求1所述基于约束机制粒子群算法的水库库容曲线修正方法,其特征在于:所述步骤1中,构建符合变量内部相邻成员间变化规律的约束机制步骤为:
子步骤1.1,通过分析水库库区地形的变化规律,提出水库水位~面积函数的微分特性为:
设水库水位变量Z,对应水面面积函数A=f(Z),根据库区地形的一般变化规律可知,随着水库水位的抬高,相应水面面积越大,即面积函数的一阶导数大于零,有:
A′=f′(Z)=dA/dZ>0 式(1)
另外,根据水位越高水面面积增长率也越来越大,即面积函数的二阶导数也大于零,因此有:
A″=f″(Z)=d2A/dZ2>0 式(2)
式2可以看出,面积函数为严格凹函数,由其性质可知,对于x1<x2,存在0<α<1,有:
f(αx1+(1-α)x2)<αf(x1)+(1-α)f(x2) 式(3)
子步骤1.2,基于子步骤1.1的面积函数微分特性,等差水位相邻离散点距的相互制约关系为:
设库容曲线中离散水位值为Zj,(j=0,1,2,……,N),从小到大排列的序列,对应水面面积系列为Aj=f(Zj),由式(1)有:
(Aj-Aj-1)/(Zj-Zj-1)>0,
因Zj>Zj-1,有Aj-Aj-1>0,即:
0≤Aj-1<Aj(j=1,2,……,N) 式(4)
由式(2)有:
对等差水位序列,简化为:
子步骤1.3,针对某些库区的水位~面积曲线存在面积函数局部区间为非凹函数时,约束条件变更判别式及非凹函数情况下的制约关系为:
对于遵循凹函数规律的水位~面积曲线而言,式(5)是严格成立的,有(2×Aj-1-Aj-2)<(Aj+1+Aj-1)/2,整理得:
Aj+1>(3×Aj-1-2×Aj-2)(j=2,3,......,N-1) 式(6)
因此,式(6)是判别Aj有无满足凹函数性质解的条件;对于部分特殊位置的地形,存在不满足凹函数的情况,即式(5)无解或式(6)不成立;
所以式(6)也是约束条件变更的判别式,若式(6)不成立,说明在该库区高程范围对应的面积增长率非正值,此时的面积函数应变更为非凹函数约束,即:
(Aj+1+Aj-1)/2≤Aj<Aj+1 (j=1,2,......,N-1) 式(7)
子步骤1.4,初始化水位~面积曲线时,由于Aj+1未知,无法利用式(5)~式(7)来约束Aj,根据式(3),其离散表达形式为:
Aj<A0+(AN-A0)×(Zj-Z0)/(ZN-Z0) 式(8)
对等差水位序列,进一步简化为:
Aj<A0+(AN-A0)×j/N(j=1,2,......,N-1) 式(9)
由此,在初始化水位~面积曲线时,利用式(8)或式(9)来代替式5中的约束条件上边界。
3.根据权利要求1所述基于约束机制粒子群算法的水库库容曲线修正方法,其特征在于:所述步骤2中,所述目标水库的运行资料包括库水位、入库流量、出库流量;所述气象资料包括库区降雨量、蒸发量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010006199.4A CN111210141B (zh) | 2020-01-03 | 2020-01-03 | 基于约束机制粒子群算法的水库库容曲线修正方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010006199.4A CN111210141B (zh) | 2020-01-03 | 2020-01-03 | 基于约束机制粒子群算法的水库库容曲线修正方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111210141A true CN111210141A (zh) | 2020-05-29 |
CN111210141B CN111210141B (zh) | 2023-07-14 |
Family
ID=70788089
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010006199.4A Active CN111210141B (zh) | 2020-01-03 | 2020-01-03 | 基于约束机制粒子群算法的水库库容曲线修正方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111210141B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112418491A (zh) * | 2020-11-06 | 2021-02-26 | 黄河勘测规划设计研究院有限公司 | 一种水库剩余拦沙库容动态配置方法 |
CN115907536A (zh) * | 2022-11-22 | 2023-04-04 | 国能大渡河流域水电开发有限公司 | 一种基于区间谬误率的水库特征曲线校正方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104613943A (zh) * | 2013-11-04 | 2015-05-13 | 中国水利水电科学研究院 | 水库蓄水量遥感与地面协同监测方法 |
CN105178240A (zh) * | 2015-06-08 | 2015-12-23 | 武汉大学 | 一种绘制p-iii型分布频率曲线的优化方法 |
US20180357584A1 (en) * | 2017-06-12 | 2018-12-13 | Hefei University Of Technology | Method and system for collaborative scheduling of production and transportation in supply chains based on improved particle swarm optimization |
CN110598983A (zh) * | 2019-08-08 | 2019-12-20 | 华中科技大学 | 一种自适应改进粒子群算法的梯级水库优化调度方法 |
-
2020
- 2020-01-03 CN CN202010006199.4A patent/CN111210141B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104613943A (zh) * | 2013-11-04 | 2015-05-13 | 中国水利水电科学研究院 | 水库蓄水量遥感与地面协同监测方法 |
CN105178240A (zh) * | 2015-06-08 | 2015-12-23 | 武汉大学 | 一种绘制p-iii型分布频率曲线的优化方法 |
US20180357584A1 (en) * | 2017-06-12 | 2018-12-13 | Hefei University Of Technology | Method and system for collaborative scheduling of production and transportation in supply chains based on improved particle swarm optimization |
CN110598983A (zh) * | 2019-08-08 | 2019-12-20 | 华中科技大学 | 一种自适应改进粒子群算法的梯级水库优化调度方法 |
Non-Patent Citations (4)
Title |
---|
KAZEM SABERCHENARI & HIRAD ABGHARI & HOSSEIN TABARI: "Application of PSO algorithm in short-term optimization of reservoir operation" * |
王芳;杨虎山;: "基于混合DE-PSO算法的水库优化调度模型及其应用" * |
苏业助,洪为善: "根据水量平衡原理修正水库库容曲线方法" * |
许辉熙;薛万蓉;何政伟;张东辉;: "基于DEM和"土地平整法"的水库面积-库容曲线计算系统设计与开发" * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112418491A (zh) * | 2020-11-06 | 2021-02-26 | 黄河勘测规划设计研究院有限公司 | 一种水库剩余拦沙库容动态配置方法 |
CN112418491B (zh) * | 2020-11-06 | 2021-08-03 | 黄河勘测规划设计研究院有限公司 | 一种水库剩余拦沙库容动态配置方法 |
CN115907536A (zh) * | 2022-11-22 | 2023-04-04 | 国能大渡河流域水电开发有限公司 | 一种基于区间谬误率的水库特征曲线校正方法 |
CN115907536B (zh) * | 2022-11-22 | 2023-10-31 | 国能大渡河流域水电开发有限公司 | 一种基于区间谬误率的水库特征曲线校正方法 |
Also Published As
Publication number | Publication date |
---|---|
CN111210141B (zh) | 2023-07-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111582755B (zh) | 一种基于多维度集合信息山洪灾害综合风险动态评估方法 | |
CN107730151B (zh) | 一种基于概念性水文模型的流域设计洪水推求方法 | |
CN112149314B (zh) | 一种基于虚拟库容修正的多沙水库库容冲淤模拟方法 | |
CN107045568B (zh) | 基于动态规划逐次逼近法的河道糙率反演方法 | |
CN111241758B (zh) | 基于可溶性污染物在水环境中输移扩散模型的评估方法 | |
CN110598290A (zh) | 考虑气候变化的流域未来水电发电能力预测方法和系统 | |
CN111210141A (zh) | 基于约束机制粒子群算法的水库库容曲线修正方法 | |
CN111625922A (zh) | 一种基于机器学习代理模型的大规模油藏注采优化方法 | |
CN112287436B (zh) | 一种多沙河流水库淤积断面和有效库容设计方法及系统 | |
CN113095694B (zh) | 一种适用于多地貌类型区的降雨输沙模型构建方法 | |
CN112418491A (zh) | 一种水库剩余拦沙库容动态配置方法 | |
CN112949895A (zh) | 一种基于动态可扩展神经网络模型的风速预测方法 | |
CN113486556B (zh) | 一种改进的油气藏高效自动历史拟合方法 | |
CN108427654B (zh) | 一种中型以上淤地坝已淤积库容快速演算方法 | |
CN102236746B (zh) | 无测风记录区风资源模拟推算方法 | |
CN110059443B (zh) | 一种分层水库取水下泄水温的快速预测方法 | |
CN110569622B (zh) | 一种基于多目标优化的挡土墙优化设计方法 | |
Khalaf et al. | Multi-objective groundwater management using genetic algorithms in Kerbala desert area, Iraq | |
CN115169243A (zh) | 基于ga-pso-glssvm算法的土岩复合地层深基坑变形时序预测方法 | |
CN113283183A (zh) | 基于dpso-anfis的弃渣坝变形预测方法 | |
CN113139692A (zh) | 基于dpso-anfis的弃渣坝变形预测方法 | |
CN108664711A (zh) | 锚杆轴力变化趋势预测方法 | |
CN112347700A (zh) | 一种泥石流发生概率及规模预报方法 | |
Wang et al. | A new treatment of depression for drainage network extraction based on DEM | |
CN117332544B (zh) | 矢量与栅格水文计算单元协同的城市雨洪模型建模方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |