库群联合调度的轮库寻优程序.docx

库群联合调度的轮库寻优程序 摘要：按一定目标对库群进行联合调度，具有较大的实用价值州河地处津冀交界处，流域水资源对河北省地方经济发展及天津市用水十分重要运用轮库寻优算法编写程序，对州河流域三库进行联合调度，可有效地解决流域水资源优化利用问题 关键词：联合调度 库寻优算法 程序 1 前言随着地区经济的迅速发展，区域水资源短缺程度日趋严重，各地区在水资源分配与利用上的矛盾越发突出如何挖掘区域内水资源潜力，充分利用现有水利工程，实现库群水资源的联合优化调度，具有较高的科研与实用价值这里，以州河流域内三库水资源联合优化调度为例，编写轮库寻优计算程序，使现实问题得以解决，在库群联合调度方面进行了有益的探索2 流域状况与资料年限2.1 流域状况州河流域北部和东部与滦河为邻，南部和西部分别与泃河水系接壤流域呈扇形，东西长约56 km，南北宽近50 km，总面积2 114 km2流域内岩土以第四纪冲击层为主，北部山区多片麻岩，土壤为淋溶褐土，森林茂密，植被率达50%～60%；中部为冲积平原，土壤多为砂质粘土，农业发达流域属大陆性季风气候，多年平均降水量750 mm降水的年内、年际分配极不均匀年内降水主要集中在6～9月份，占全年的80%以上。

年际变化很大，丰水年多达1 256 mm ，枯水年仅有387 mm州河为蓟运河左支，上游干游为魏进河自东向西有黎河、沙河、淋河等3条主要支流，集水面积分别为560 km2、887 km2、252 km2三支流及干流经国河汇入于桥水库，出库后始称州河，向南汇流右支泃河后称蓟运河三支流均为山溪性河流，地处燕山迎风坡暴雨中心区，每遇汛期暴雨，洪水较大流域内主要水利工程为3座大、中型水库，即上关、般若院及于桥水库，详见图1图1 州河流域水利工程示意图2.2 资料年限所采用的水文资料及成果，为于桥水库1960—1999年、上关水库1974—2000年、般若院水库1973—2000年逐年各月入库径流量资料及其计算成果，流域内各水文站1960—1999年的观测资料资料年限最长达40年，有利于对州河流域三水库水资源进行联合调度分析3 联合调度的数学模型3.1 最优准则与目标函数研究库群联合优化调度必先确定最优准则这里根据州河流域三库水文特性及水资源利用现状，确定最优准则为：三水库控制并利用水资源量最大；三水库内水面蒸发量最小；于桥水库向下游排放水量最小上关、般若院及于桥水库自起始调度月起，水库各月平均水面面积分别记作：S11 ，S12，…，S112 ；S21 ，S22，…，S212 ；S31 ，S32，…，S312 。

各月平均水面蒸发系数与蒸发量分别记作：K1 ，K2，…，K12 E1 ，E2，…，E12 ；于桥水库各月出库水量分别记作：QO31 ，QO32，…，QO312 则由最优准则可写出求解三库水量损失的数学表达式，亦即目标函数如下：Min∑[(S1i+S2i+S3i)·Ki·Ei/10+ QO3i ]　　　　i=1，2，…，123.2 约束方程（1）蓄水量约束Wji≤wji≤Wji　　j=1，2，3； i=1，2，…，13式中，wji为j水库第i个月份的库蓄水量；Wji、Wji分别为j水库第i个月份的允许最小、最大库蓄水量2）需供水量约束Wji+Uji≤wji　　j=1，2，3； i=1，2，…，13式中，Uji为j水库第i个月份需供水量，其它同上3）起调蓄量约束wj1=wj13　　j=1，2，3式中，wj1、wj13分别为j水库第1个月份、第13个月份库蓄水量，均为起调蓄量3.3 数学模型先确定三库联合调度的决策变量为：起调蓄量wj1、入库水量QIwji、需供水量Uwji、蒸发系数Ki、蒸发量Ei，则综合目标函数及约束方程，可获得州河流域三库联合调度的优化数学模型如下：MINEQ=Min∑[(S1i+S2i+S3i)·Ki·Ei /10+QO3i]s.t. wj1=wj13　　 Wji+Uji≤wji≤Wji 式中，j=1，2，3；i=1，2，…，13。

4 轮库寻优算法与程序流程示意图4.1 轮库寻优算法由州河流域三库数学模型求取目标函数，可采用轮库寻优算法，其思路如下：（1）根据一般经验、分析判断或用其它简便方法，先给三水库定出起调蓄量Wj1由起调蓄量，可对三水库分别确定一条满足约束条件且各月不超过最高蓄量的初始调度线：Wj1 ，Wj2，…，Wj12，计算目标函数值2）固定般若院水库的初始调度线，再将一定步长的蓄量变化△W1作为上关水库向于桥水库的放水量QO1，重新确定上关及于桥水库相应的蓄量调度线，并计算其目标函数值之后，逐次对上关水库进行减量优化调度，比较各次目标函数值，记录最小水量损失及相应三库调度线，直至不满足上关水库的约束条件为止3）将一定步长的蓄量变化△W2作为般若院水库向于桥水库的放水量QO2，并固定般若院水库减量后的蓄量调度线再重复进行（2）中对上关水库的减量优化调度步骤这样反复轮换优化于桥水库上游的上关及般若院两水库，直至不满足般若院水库的约束条件为止4）按上述步骤，可计算出相应步长三库所有蓄量调度线的目标函数值，进而寻得最小目标函数值（即最少损失水量）及三库相应的蓄量调度线、放水量4.2 程序流程示意图对上述轮库寻优算法的求解步骤，设计其计算机程序流程示意图，见图2。

图2 州河流域三库优化调度程序流程示意图 4.3 程序清单分别用W[i][j]、S[i][j]、QI[i][j] 、QO[i][j]、UW[i][j]表示i水库第j个月份的蓄水量、库水面面积、入库水量、出库水量、需供水量，用DW[i] 、steps表示i水库的死库容、寻优步长，用K[j]、E[j]分别表示第j个月份的蒸发系数、蒸发量，用Min_EQ （MINEQ）、Min_E（MINE）、Min_Q（MINQ）分别表示三库最小水量损失、最小蒸发水量、于桥水库最小出库水量，用WL[i][j]、QOL[i][j]表示i水库第j个月份的最优蓄水量、出库水量对轮库寻优算法的求解步骤，则可用C语言编程如下：int i，j; //初始化寻优起始状态　　for(i=1;i<4;i++){　　　　for(j=1;j<14;j++)　　WL[i][j]=0;}　　MINEQ= 99999999; MINE= 99999999;　　 MINQ=99999999;　　W[1][1]=?，W[2][1]=?，W[3][1]=?; //轮库寻优算法　　int hh=0，tt=0;　　for(int fp=0;;fp++){//按一定步长要求生成般若院水库蓄水量调度线　　　　for(int g=2;g<14;g++) W[2][g]=0;　　　　for(int gg=1;gg<13;gg++)　　 QO[2][gg]=0;　　　　for(int h=2;h<14;h++){//过量蓄量调整　　　　　　QO[2][h-1]=steps*fp;　　　　　　W[2][h]=W[2][h-1]+QI[2][h-1]-steps*fp;　　　　　　if(W[2][h]>2800){　　　　　　　　QO[2][h-1]=QO[2][h-1]+W[2][h]-2800;　　　　　　　　W[2][h]=2800;}　　　　　　if(W[2][13]>W[2][1]){　　　　　　　　QO[2][12]=W[2][13]-W[2][1]+QO[2][12];　　　　　　　　W[2][13]=W[2][1];}}　　　　for(i=1;i<13;i++){//判断并结束般若院水库过程线　　　　　　if(W[2][i]< DW[2]+UW[i]) break;}　　　if(W[2][13]<> 　　　　for(int fs=0;;fs++){//按一定步长要求生成上关水库蓄水量调度线　　　　　　for(g=2;g<14;g++)　　　　W[1][g]=0;　　　　　　for(gg=1;gg<13;gg++)　　　　QO[1][gg]=0;　　　　　　for(int t=2;t<14;t++){　　　　　　　　QO[1][t-1]=steps*fs;　　　　　　　　W[1][t]=W[1][t-1]+ QI[1][t-1]-steps*fs;　　　　　　　　if(W[1][t]>2960){　　　　　　　　　　QO[1][t-1]=QO[1][t-1]+W[1][t]-2960;　　　　　　　　　　W[1][t]=2960;} 　　　　　　　　if(W[1][13]>W[1][1]){　　　　　　　　　　QO[1][12]=W[1][13]-W[1][1]+QO[1][12];　　　　　　　　　　W[1][13]=W[1][1];}}　　　　　　for(i=1;i<13;i++){//判断并结束上关水库过程线　　　　　　　　if(W[1][i]< DW[1]+ UW[i])　　break;}　　　　　　if(W[1][13]<> 　　　　　　for(g=2;g<14;g++)　　　　W[3][g]=0;　　　　　　for(gg=1;gg<13;gg++) QO[3][gg]=0;　　　　　　for(int r=2;r<14;r++){//生成于桥水库蓄水量过程线　　　　　　　W[3][r]=W[3][r-1]+app->QI[3][r-1]+QO[1][r-1]+QO[2][r-1];　　　　　　　　if(W[3][r]>59791){　　　　　　　　　　QO[3][r-1]=W[3][r]-59791;　　W[3][r]=59791;}　　　　　　　　if(W[3][13]>W[3][1]){　　　　　　　　　　QO[3][12]=QO[3][12]+(W[3][13]-W[3][1]);　　　　　　　　W[3][13]=W[3][1];}}　　　　　　 //计算并判断最小水量损失，记录相应库状态　　　　　　for(j=1;j<13;j++){//库水面面积与库容关系模拟　　　　　　　　S[1][j]= 0.1*pow(10，-9)*pow(W[1][j]，3) - 0. 6*pow(10，-6)*pow(W[1][j]，2) +0.0014*W[1][j]+ 0.2583; 　　　　　　　　S[2][j]=0.6*pow(10，-10)*pow(W[2][j]，3)-0.5*pow(10，-6)*pow(W[2][j]，2)+ 0.0021*W[2][j]+ 0.2345;　　　　　　　　S[3][j]=0.5*pow(10，-12)*pow(W[3][j]，3)-0.7*pow(10，-7)*pow(W[3][j]，2)+ 0.0043*W[3][j]+ 6.8841;}　　　　　　　　 Min_EQ=0，Min_E=0，Min_Q=0;　　　　　　for(int jj=1;jj<13;jj++){//计算目标函数　　　　　　　　Min_EQ=Min_EQ+(SO[1][jj]+SO[2][jj]+SO[3][jj])* K[jj]*E[。