作者简介
很高兴认识您！
我叫卢家波，河海大学水文学及水资源博士研究生，研究兴趣为高效洪水淹没预测、洪水灾害预警、机器学习、替代模型和降阶模型。
变化环境下，极端洪水事件多发，我希望能通过研究为水灾害防御做出贡献，为人民服务。
我很乐意与您交流和研究合作，vx Jiabo_Lu。
主页 https://lujiabo98.github.io
简历 https://lujiabo98.github.io/file/CV_JiaboLu_zh.pdf
博客 https://blog.csdn.net/weixin_43012724?type=blog
来信请说明博客标题及链接，谢谢。
一、基本条件
- 以溪洛渡-葛洲坝梯级电站发电量最大为优化目标；
- 约束仅为各库下游生态基流，取消原保证出力约束；
- 来水条件为1998年实测来水；
- 算法采用缩放因子 F 0 F_0 F0和交叉率 C r Cr Cr自适应的差分进化算法。
二、稳定性测试
2.1 测试方法
运行优化程序20次，统计发电量平均值、方差，观察各次调度方案差异。
2.2 测试结果
| 平均值（亿kW·h） | 标准差(亿kW·h) |
|---|---|
| 2040.522 | 12.807 |
其中一次的结果
三、结果分析
仅考虑生态基流约束，发电量较为稳定，且各库都能满足约束。下一步考虑对不确定条件下水库群多目标优化问题建模。
四、程序代码
GroupSingleStableTest.m
%-----------------------备注-------------------------
%2020.1.16
%溪洛渡、向家坝、三峡、葛洲坝梯级水库群1998年来水资料
%参数F0、CR的自适应调整
%稳定性测试
%目标函数：发电量最大
%约束：生态基流
%max f = K Sigma Sigma H(i,j) Q(i,j) delta t
%调度期为1年，每旬为一个时段，共36个时段
%-----------------------参数定义及初始化-------------------------
close all;
clc,clear;
NP=70; %种群个体数量
D=37*4; %染色体长度，这里有36旬，共37个节点蓄水量，共四库
G=2000; %迭代代数
F0=0.5; %缩放因子[0,2],一般取0.5
Fl=0.1;Fu=0.7;
CR=0.1; %交叉率[0,1],一般取0.3，减小会更快收敛，造成早熟
CRi=zeros(NP,4);
CRl=0.3;CRu=0.7;
punish=10^13; %罚因子
trace=zeros(1,G+1);
K=[8.5 8.5 8.5 8.5]; %水库出力系数
Nmin=[379.5*10^4 200.9*10^4 499*10^4 76.8*10^4]; %四库保证出力 单位：kW
% Nmin=[1 1 1 1]; %测试
Nmax=[1260*10^4 600*10^4 2250*10^4 271.5*10^4]; %四库装机容量 单位：kW
Nmax=kron(Nmax,ones(1,36)); %扩充Nmax向量
qmax=[7740 7132 25900 17900]; %水轮机组最大过水流量 单位：m3/s
qmax=kron(qmax,ones(1,36)); %扩充qmax向量
%出库流量过程
q1=zeros(NP,D/4-1); %溪洛渡出库流量
q2=zeros(NP,D/4-1); %向家坝出库流量
q3=zeros(NP,D/4-1); %三峡出库流量
q4=zeros(NP,D/4-1); %葛洲坝出库流量
q=[q1 q2 q3 q4]; %四库出库流量过程
%出力过程
N=zeros(NP,D-4);
%尾水曲线的三个参数 Zd=a2*q^2+b2*q+d2 q单位：立米秒
a2=[-0.0000000064650 -0.0000000111000 0.0000000011101 -0.0000000026237 ];
b2=[0.001181 0.00112 0.0000836 0.00045252];
d2=[356 265 65.5 38.052 ];
pH=[0.61 1.19 1.07 5.72]; %四库平均单耗 单位：m3/s/MW
pH=kron(pH,ones(1,36)); %扩充qmax向量
%单耗曲线三参数 K=a3*H^2+b3*H+d3 K：单耗，H：水头
a3=[0.0000793 0.0007590 0.0008990 0.0883000];
b3=[-0.0432 -0.199 -0.219 -5.120];
d3=[7.48 16.3 16.9 89.2];
%水库的水位库容曲线 Zu=a1*V^2+b1*V+d1 V单位：亿立米
zu=[130 135 140 145 150 155 160 165 170 175 185];
V=[103.3 124.1 147 171.5 196.9 228 262 300.2 344 393 445.7];
a1=[-0.0022098 -0.0065343 -0.00017085 -0.34965];
b1=[1.2867 1.7 0.233 8.09];
d1=[481 312 110 26.4 ];
%四库各旬初的最低蓄水量与最高蓄水量
Vmin1=[50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18 50.18];
Vmax1=[115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 93.51 69.75 69.75 69.75 69.75 69.75 69.75 69.75 102.90 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33 115.33];
Vmin2=[40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39 40.39];
Vmax2=[49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 40.39 40.39 40.39 40.39 40.39 40.39 40.39 41.24 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37 49.37];
Vmin3=[228,228,228,228,228,228,228,228,228,228,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,171.5,228];
Vmax3=[393,393,393,393,393,393,393,393,393,393,393,393,393,393,344,228,179.12,179.12,179.12,179.12,179.12,179.12,179.12,179.12,179.12,300.2,300.2,393,393,393,393,393,393,393,393,393,393];
Vmin4=[6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17 6.17];
Vmax4=[7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03 7.03];
Vmin=[Vmin1 Vmin2 Vmin3 Vmin4];
Vmax=[Vmax1 Vmax2 Vmax3 Vmax4];
%龙头水库溪洛渡2013年来水，其余三库区间入流 单位：立米秒
% Q1=[1859.597656 1868.240723 1594.209229 1439.341187 1449.626099 1567.817505 1754.30127 1770.719971 1747.469116 2129.174072 2598.531738 1869.13916 1930.608521 1982 2462.727295 3956 4361 3721 4666 6517 10163.63672 7371 6394 5144.54541 6390 9605 5673 5251 4917 4402.727051 3107 3038 2774 2281 1965 1792];
% Q2=[16.88837433 16.96686935 14.47818661 13.07171631 13.16512108 14.23850346 15.9321003 16.08121109 15.87005234 19.33659554 23.59917831 16.97502708 17.5332756 18 22.365839 35.92734528 39.60544968 33.7931366 42.37537766 59.18567276 92.30345917 66.94147491 58.06861877 46.72140121 58.03229141 87.23007202 51.5206871 47.68819427 44.65489197 39.98440552 28.21695328 27.59031296 25.19273376 20.71543884 17.84560966 16.27447128 ];
% Q3=[3447.013916 3466.192383 2858.130859 2514.487061 2537.308838 2799.569092 3213.366699 3249.798828 3198.206299 4045.189453 5086.668945 3468.185791 5044.358398 9091 10432.54297 10785.77246 11165.59473 13304.80664 23107.42383 26221.71484 23595.24219 16008.3584 13596.63184 9604.551758 11330.56836 11944.87012 11541.67969 6455.312012 5481.544922 4576.561035 4589.683105 4178.209961 3361.707275 3121.884521 3488.054199 3441.925537];
% Q4=[16.88837433 16.96686935 14.47818661 13.07171631 13.16512108 14.23850346 15.9321003 16.08121109 15.87005234 19.33659554 23.59917831 16.97502708 17.5332756 18 22.365839 35.92734528 39.60544968 33.7931366 42.37537766 59.18567276 92.30345917 66.94147491 58.06861877 46.72140121 58.03229141 87.23007202 51.5206871 47.68819427 44.65489197 39.98440552 28.21695328 27.59031296 25.19273376 20.71543884 17.84560966 16.27447128];
%龙头水库溪洛渡1998年来水，其余三库区间入流 单位：立米秒
data= xlsread('三峡梯级联合调度数据分析.xlsm',1, 'BH3:BH119')';
Q1 = data(:,1:36);
Q2 = data(:,41:76);
Q3 = data(:,82:117);
Q4= zeros(1,36);
Q= [Q1 Q2 Q3 Q4];
%四库各旬的最小下泄流量
qmin1=[1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600];
qmin2=[1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600];
qmin3=[6000 6000 6000 6000 6000 6000 6000 6000 6000 6000 6000 6000 5700 5700 5700 5700 5700 5700 ...
5700 5700 5700 5700 5700 5700 5700 5700 10000 10000 8000 8000 5700 5700 5700 5700 5700 5700];
qmin4=[5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700 5700];
qmin=[qmin1 qmin2 qmin3 qmin4];
%各旬的时段长，h
dt=[240 240 264 240 240 192 240 240 264 240 240 240 240 240 264 240 240 240 ...
240 240 264 240 240 264 240 240 240 240 240 264 240 240 240 240 240 240];
Z0=[600 380 175 66]; %四库调度期初水位约束
ZT=[600 380 174.5 66]; %四库调度期末水位约束
%由库容曲线求根公式得到
V0=zeros(1,4);VT=zeros(1,4);
for i = 1 : 4
V0(i)=(-b1(i)+sqrt(b1(i)^2-4*a1(i)*(d1(i)-Z0(i))))/(2*a1(i)); %四库调度期初蓄量约束 单位：亿立米
VT(i)=(-b1(i)+sqrt(b1(i)^2-4*a1(i)*(d1(i)-ZT(i))))/(2*a1(i)); %四库调度期末蓄量约束 单位：亿立米
end
x=zeros(NP,D); %初始种群
v=zeros(NP,D); %变异种群
u=zeros(NP,D); %选择种群
%种群初始化
for i=1:1:NP
for j=2:1:D-1
%随机产生本时段初蓄水量
x(i,j)=Vmin(j)+rand*(Vmax(j)-Vmin(j));
end
end
%调度期初、末蓄量约束
for i = 1 : 1 : 4
x(:,1+37*(i-1))=V0(i); x(:,37*i)=VT(i);
end
[ob,ob_ind,q,N,E] = Cfit(x,dt,Q1,Q2,Q3,Q4,K,Nmax,Nmin,qmin,a1,b1,d1,a2,b2,d2,punish); %适应度矩阵每一行的和
trace(1)=max(ob); %记录下适应度最大值
%差分进化循环
for gen=1:G
% 变异操作
for m=1:NP
% 产生不同的r1,r2,r3
r1=randi([1,NP],1,1); %从种群NP中随机选择一个个体
while(r1==m)
r1=randi([1,NP],1,1);
end
r2=randi([1,NP],1,1);
while(r2==r1)||(r2==m)
r2=randi([1,NP],1,1);
end
r3=randi([1,NP],1,1);
while(r3==m)||(r3==r2)||(r3==r1)
r3=randi([1,NP],1,1);
end
temp = [x(r1,:);x(r2,:);x(r3,:)]; %临时矩阵存储三个随机选择的个体
[tempfit,tempfit_ind,~,~,~] = Cfit(temp,dt,Q1,Q2,Q3,Q4,K,Nmax,Nmin,qmin,a1,b1,d1,a2,b2,d2,punish); %适应度矩阵每一行的和
[tempfit_ind,index]=sort(tempfit_ind); %适应度排序
Fi=Fl+(Fu-Fl)*((tempfit_ind(2,:)-tempfit_ind(3,:))./(tempfit_ind(1,:)-tempfit_ind(3,:))); %参数的自适应调整
for i = 1 : 4
k=1+37*(i-1);l=37*i; %计算状态量循环索引
% 这个变异是针对单个个体的变异
v(m,k:l)=temp(index(3,i),k:l)+Fi(i)*(temp(index(2,i),k:l)-temp(index(1,i),k:l)); %变异
end
end
%在交叉操作前先初始化新种群u
u = x;
max_ob=max(ob_ind);min_ob=min(ob_ind);mean_ob=mean(ob_ind);
for m =1:NP
for j = 1:4
if ob_ind(m,j)<mean_ob(j)
CRi(m,j)=CRl;
else
CRi(m,j)=CRl+(CRu-CRl)*((ob_ind(m,j)-min_ob(j))/(max_ob(j)-min_ob(j)));
end
end
end
%四库交叉操作
for m = 1:1:NP
for i = 1 : 4
k=1+37*(i-1);l=37*i; %计算状态量循环索引
% 这个变异是针对单个个体的变异
for n=k+1:l-1
cr=rand;
if (cr>=CRi(m,i))
u(m,n)=v(m,n);
end
end
end
end
%边界条件处理
%中间各旬蓄量约束
for m=1:NP
for i = 1 : 4
k=1+37*(i-1);l=37*i; %计算状态量循环索引
for n=k+1:l-1
if (u(m,n)<Vmin(n))||(u(m,n)>Vmax(n))
u(m,n)=Vmin(n)+rand*(Vmax(n)-Vmin(n));
end
% if u(m,n)<Vmin(n)
% u(m,n)=Vmin(n);
% end
% if u(m,n)>Vmax(n)
% u(m,n)=Vmax(n);
% end
end
end
end
% 自然选择
% 计算新的适应度
[ob_1,~,~,~,~] = Cfit(u,dt,Q1,Q2,Q3,Q4,K,Nmax,Nmin,qmin,a1,b1,d1,a2,b2,d2,punish); %适应度矩阵每一行的和
for m=1:NP
if ob_1(m)>ob(m)
x(m,:)=u(m,:);
end
end
% 现在x为经过选择后的种群
[ob,ob_ind,q,N,E] = Cfit(x,dt,Q1,Q2,Q3,Q4,K,Nmax,Nmin,qmin,a1,b1,d1,a2,b2,d2,punish); %适应度矩阵每一行的和
trace(gen+1)=max(ob); %记录下适应度最大值
end
[obmax,p]=max(ob); %求出最大适应度和位置
xfit=x(p,:); %记录最适合的库蓄量调度序列
qfit=q(p,:); %记录最适合的出库流量序列
% qfit(qfit<qmin)=0;
for i = 1 : 4
k=1+37*(i-1);l=37*i; %计算循环索引
Zufit(:,k:l)=a1(i)*xfit(:,k:l).^2+b1(i)*xfit(:,k:l)+d1(i); %记录最适合的库水位调度序列
end
Nfit=N(p,:)/10^4; %记录最适合的出力序列 单位：万kW
Nmin=kron(Nmin/10^4,ones(1,36)); %扩充Nmin向量 单位：万kW
Efit=E(p,:)/10^8; %记录最适合的发电量序列 亿kW·h
Esum=sum(Efit); %一年的发电量 亿kW·h
%绘制寻优过程
figure(1);
plot(trace);
title(['发电量最大值: ', num2str(Esum),'亿kW·h']);
xlabel('迭代次数');
ylabel('适应度值/亿kW·h');
print(gcf,'-djpeg','梯级水库群迭代图');
for j = 1 : 4
m=1+36*(j-1);n=36*j; %计算过程量循环索引
k=1+37*(j-1);l=37*j; %计算状态量循环索引
%判断绘图图名
if j == 1
name = '溪洛渡';
else
if j == 2
name = '向家坝';
else
if j == 3
name = '三峡';
else
name = '葛洲坝';
end
end
end
% %绘制解集
% figure(2+5*(j-1));
% plot(xfit(:,k:l));
% title([name,'水库蓄量变化过程 ']);
% xlabel('旬');
% ylabel('蓄量/亿m^3');
% print(gcf,'-djpeg',[name,'水库蓄量变化过程 ']);
%绘制出库流量过程
figure(3+5*(j-1));
plot(qfit(:,m:n));
title([name,'水库出库流量过程 ']);
xlabel('旬');
ylabel('流量(m^3/s)');
fig=qfit(:,m:n);
locat=find(fig<qmin(:,m:n));
for i = 1 : length(locat)
tag=text(locat(i),fig(locat(i)),['(',num2str(locat(i)),',',num2str(round(fig(locat(i)))),')'],'color','b');
set(tag,'FontSize',8) %设置字号8
end
print(gcf,'-djpeg',[name,'水库出库流量过程 ']);
% %绘制水位变化过程
% figure(4+5*(j-1));
% plot(Zufit(:,k:l));
% title([name,'水库水位变化过程 ']);
% xlabel('旬');
% ylabel('水位/m');
% print(gcf,'-djpeg',[name,'水库水位变化过程 ']); %保存图片
% %绘制出力过程
% figure(5+5*(j-1));
% plot(Nfit(:,m:n));
% title([name,'水库出力变化过程，保证出力： ',num2str(round(Nmin(m),1)),'万kW']);
% xlabel('旬');
% ylabel('出力/万kW');
% fig=Nfit(:,m:n);
% locat=find(fig<Nmin(:,m:n));
% for i = 1 : length(locat)
% tag=text(locat(i),fig(locat(i)),['(',num2str(locat(i)),',',num2str(round(fig(locat(i)))),')'],'color','b');
% set(tag,'FontSize',8)% 设置字号8
% end
% print(gcf,'-djpeg',[name,'水库出力变化过程 ']);
% %绘制发电量变化过程
% figure(6+5*(j-1));
% plot(Efit(:,m:n));
% title([name,'水库发电量变化过程 ']);
% xlabel('旬');
% ylabel('发电量/亿kW·h');
% print(gcf,'-djpeg',[name,'水库发电量变化过程 ']);
end
%适应度计算函数，Nmin替换为~
function [fit,fit_ind,q,N,E] = Cfit(x,dt,Q1,Q2,Q3,Q4,K,Nmax,~,qmin,a1,b1,d1,a2,b2,d2,punish)
%矩阵x行数
NP=size(x,1);
%矩阵x列数
D=size(x,2);
%四库出库流量过程
q=zeros(NP,D-4);
%出力过程
N=zeros(NP,D-4);
%发电量过程
E=zeros(NP,D-4);
%各库单独的适应度
fit_ind=zeros(NP,4);
%q为四库出库流量
for i = 1:1:NP
q(i,1:36)=-(diff(x(i,1:37)))*10^8./(3600*dt) + Q1; %q1(i,:)
j = 2;
q(i,1+36*(j-1):36*j)=-(diff(x(i,1+37*(j-1):37*j)))*10^8./(3600*dt) + q(i,1+36*(j-2):36*(j-1)) + Q2;
j = 3;
q(i,1+36*(j-1):36*j)=-(diff(x(i,1+37*(j-1):37*j)))*10^8./(3600*dt) + q(i,1+36*(j-2):36*(j-1)) + Q3;
j = 4;
q(i,1+36*(j-1):36*j)=-(diff(x(i,1+37*(j-1):37*j)))*10^8./(3600*dt) + q(i,1+36*(j-2):36*(j-1)) + Q4;
end
for j =1 : 4
m=1+36*(j-1);n=36*j; %计算循环索引
Vmean(:,m:n)=(conv2(x(:,1+37*(j-1):37*j)',[1;1],'valid'))'/2; %用conv2求相邻旬蓄量平均值
Zu(:,m:n)=a1(j)*Vmean(:,m:n).^2+b1(j)*Vmean(:,m:n)+d1(j); %根据库容曲线内插得到上游平均水位
Zd(:,m:n)=a2(j)*q(:,m:n).^2+b2(j)*q(:,m:n)+d2(j); %根据尾水曲线公式生成下游水位
N(:,m:n)=K(j)*q(:,m:n).*(Zu(:,m:n)-Zd(:,m:n)); %各旬的出力 单位：kW
N(N>Nmax)=Nmax(m); %出力超过装机容量的旬时段改为装机容量
E(:,m:n)=N(:,m:n).*dt; %各旬的发电量 单位：kW·h
Fit(:,m:n)=E(:,m:n)+punish*min((q(:,m:n)-qmin(:,m:n))./qmin(:,m:n),0); %适应度矩阵
fit_ind(:,j)=sum(Fit(:,m:n),2); %各库的适应度
end
fit=sum(Fit,2); %适应度矩阵每一行的和
end
105
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
