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

摘要：按一定目标对库群进行联合调度，具有较大的实用价值。州河地处津冀交界处，流域水资源对河北省地方经济发展及天津市用水十分重要。运用轮库寻优算法编写程序，对州河流域三库进行联合调度，可有效地解决流域水资源优化利用问题。

关键词：联合调度、库寻优算法、程序

1前言

随着地区经济的迅速发展，区域水资源短缺程度日趋严重，各地区在水资源分配与利用上的矛盾越发突出。如何挖掘区域内水资源潜力，充分利用现有水利工程，实现库群水资源的联合优化调度，具有较高的科研与实用价值。这里，以州河流域内三库水资源联合优化调度为例，编写轮库寻优计算程序，使现实问题得以解决，在库群联合调度方面进行了有益的探索。

2流域状况与资料年限

2.1流域状况

州河流域北部和东部与滦河为邻，南部和西部分别与泃河水系接壤。流域呈扇形，东西长约56km，南北宽近50km，总面积2114km2。流域内岩土以第四纪冲击层为主，北部山区多片麻岩，土壤为淋溶褐土，森林茂密，植被率达50%～60%;中部为冲积平原，土壤多为砂质粘土，农业发达。流域属大陆性季风气候，多年平均降水量750mm。降水的年内、年际分配极不均匀。年内降水主要集中在6～9月份，占全年的80%以上。年际变化很大，丰水年多达1256mm，枯水年仅有387mm。

州河为蓟运河左支，上游干游为魏进河。自东向西有黎河、沙河、淋河等3条主要支流，集水面积分别为560km2、887km2、252km2。三支流及干流经国河汇入于桥水库，出库后始称州河，向南汇流右支泃河后称蓟运河。三支流均为山溪性河流，地处燕山迎风坡暴雨中心区，每遇汛期暴雨，洪水较大。流域内主要水利工程为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，…，12

3.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。

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语言编程如下：

inti,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]=?;//轮库寻优算法

inthh=0,tt=0;

for(intfp=0;;fp++){//按一定步长要求生成般若院水库蓄水量调度线

for(intg=2;g<14;g++)w[2][g]=0;

for(intgg=1;gg<13;gg++)　　qo[2][gg]=0;

for(inth=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]

if(w[2][13]<>

for(intfs=0;;fs++){//按一定步长要求生成上关水库蓄水量调度线

for(g=2;g<14;g++)　　　　w[1][g]=0;

for(gg=1;gg<13;gg++)　　　　qo[1][gg]=0;

for(intt=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]

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(intr=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(intjj=1;jj<13;jj++){//计算目标函数

min_eq=min_eq+(so[1][jj]+so[2][jj]+so[3][jj])*k[jj]*e[jj]/10+qo[3][jj];

min_e=min_e+(so[1][jj]+so[2][jj]+so[3][jj])*k[jj]*e[jj]/10;

min_q=min_q+qo[3][jj];}

if(min_eqb[1]<{>

mineq=min_eq;　　mine=min_e;minq=min_q;

for(i=1;i<4;i++){

for(j=1;j<13;j++){

wl[i][j]=w[i][j];

qol[i][j]=qo[i][j];}}}}//上关水库优化轮库结束

}//般若院水库优化轮库结束

//输出寻优成果：mineq,mine,minq,steps,wl[i][j],qol[i][j];

5三库联合优化调度成果

由州河流域三库多年水文资料，可确定三水库在枯水年、平水年、丰水年的起调蓄量、入库水量、需供水量、蒸发系数及蒸发量等参数。这里仅以丰水年的调度为例进行运算，并最终确定三水库在最优调度时的相应状态。根据已知资料，可确定上关、般若院及于桥水库在丰水年的起调蓄量分别为：1880、2455、29777(104m3)，三库入库水量(表1)，需供水量(表.2)，蒸发系数及蒸发量(表.3)。

运用轮库寻优算法，可算得实现联合调度时三库丰水年最小水量损失为81915.23(104m3)，最小蒸发量10960.25(104m3)，于桥水库最小放水量70954.98(104m3)。三库优化调度方案中蓄量、水位及水库放水量(表.4)。

6结语

州河流域尽管有较强的产流能力，却因流域面积以及上关、般若院两水库库容相对较小，在防洪中同于桥水库进行联合调度的难度较大。这里，从径流调节及流域水资源联合调度角度出发，运用数学模型与轮库寻优算法，编写计算机程序对流域内三库进行寻优计算，使流域水资源损失减至最小，也为库群联合调度提供了实用的计算方法。