本人尝试利用蛙跳算法求解水库单目标调度，目标函数是极大值函数，在运算中适应值最大的解在全局迭代中一直未更新，调试算法 发现子群内更新时按适用度降序排列第一个最优解也未更新，不知道问题出在哪里，求大神指点。

%% 多目标混合智能蛙跳算法，先进行两库单目标调度优化，以一年为调度总时段，月为时段，以发电量最大为目标函数。

%%%定义相关参数

tic

clc

clear all

m=20; %种群分组数

n=10; %t每组青蛙包含的个数

Ne=10; %组内迭代数

smax = 5;  %最大步长

MAXGEN=30;  %种群总进化代数

d=24; %优化问题维数

pmax =980; %d维最大值

pmin =888;%d维最小值

F=m*n;

%% 生成初始解向量

a=[977,977,977,980,980,977,977,957,977,970,977,977,898,898,898,898,898,898,893,892,892,898,898,898];

b=[970,970,970,970,952,952,952,952,952,952,970,970,888,888,888,888,888,893,888,888,888,888,888,888];

X1=rand(F,1)*(a(1)-b(1))+b(1);    X13=rand(F,1)*(a(13)-b(13))+b(13);

X2=rand(F,1)*(a(2)-b(2))+b(2);    X14=rand(F,1)*(a(14)-b(14))+b(14);

X3=rand(F,1)*(a(3)-b(3))+b(3);    X15=rand(F,1)*(a(15)-b(15))+b(15);

X4=rand(F,1)*(a(4)-b(4))+b(4);    X16=rand(F,1)*(a(16)-b(16))+b(16);

X5=rand(F,1)*(a(5)-b(5))+b(5);    X17=rand(F,1)*(a(17)-b(17))+b(17);

X6=rand(F,1)*(a(6)-b(6))+b(6);    X18=rand(F,1)*(a(18)-b(18))+b(18);

X7=rand(F,1)*(a(7)-b(7))+b(7);    X19=rand(F,1)*(a(19)-b(19))+b(19);

X8=rand(F,1)*(a(8)-b(8))+b(8);    X20=rand(F,1)*(a(20)-b(20))+b(20);

X9=rand(F,1)*(a(9)-b(9))+b(9);    X21=rand(F,1)*(a(21)-b(21))+b(21);

X10=rand(F,1)*(a(10)-b(10))+b(10);X22=rand(F,1)*(a(22)-b(22))+b(22);

X11=rand(F,1)*(a(11)-b(11))+b(11);X23=rand(F,1)*(a(23)-b(23))+b(23);

X12=rand(F,1)*(a(12)-b(12))+b(12);X24=rand(F,1)*(a(24)-b(24))+b(24);

XX=[X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11,X12,X13,X14,X15,X16,X17,X18,X19,X20,X21,X22,X23,X24];%在限制水位内随机生成初始向量

%% %%全局迭代寻优

yy=zeros(1,MAXGEN);              %产生一行maxgen列的零矩阵

for i1=1:MAXGEN

for i2=1:F

fitness(i2)=maxf1(XX(i2,:)); %计算每一只青蛙的适应度值

end

%排序，找最好的，并分组

[fitsort,index]=sort(fitness,'descend');

for i3=1:F

x(i3,:)=XX(index(i3),:);

end

gx=x(1,:);                   %种群内最好的青娃

yy(i1)=fitsort(1);           %种群中最好青蛙的适应度值赋予yy(ii)中

% yy(ii)=maxf1(x(1,:));

% local=zeros(n,d);

%%

for i4=1:m

local = x(1:m:end,:);   %每隔m个分配一个青蛙到一个组

for  j=1:Ne             %每组青蛙迭代次数

j=4;

pb=local(1,:);     %组内最优

pw=local(n,:);     %组内最差

s1=rand.*(pb-pw)*2;  %采用组内最优更新

s1(find(s1>smax))=smax;

temp= pw+s1;

temp(find(temp>pmax))=pmax;

temp(find(temp

%%

if  maxf1(temp)

s1=rand.*(gx-pw)*2;      %采用全局最优更新

s1(find(s1>smax))=smax;

temp=pw+s1;

temp(find(temp>pmax))=pmax;

temp(find(temp

end

%%

if maxf1(temp)

s1=(pmax-pmin)*rands(1,d)*2;       %随机更新

s1(find(s1>smax))=smax;

temp=pw+s1;

temp(find(temp>pmax))=pmax;

temp(find(temp

end

%%

local(n,:) = temp;       %替换群中最差解

for  loc=1:n

fitlocal(loc)=maxf1(local(loc,:));

end

%%

[localsort,indexlocal]=sort(fitlocal,'descend');

for loc=1:n

localnew(loc,:) = local(indexlocal(loc),:);

end

local=localnew;

end   %结束Ne

x(i4:m:end,:) =local;

end      %结束m

%最好的青娃适配值

end         %结束MAXGEN

toc         %计算整个迭代过程持续时间

%% 结果分析

plot(yy)    %记录每次全局迭代过程中最优解

title('混合跬跳算法优化');

xlabel('总进化代数');ylabel('函数最优解');

toc

disp('yy')

目标函数为：

%% 定义梯级水库优化调度目标函数，此处以发电量最大为目标。

function E= maxf1(x)

%-----------------------------------

%--------基本数据------------------

%%%%%%%%%%%%%%%%%万家寨坝址，每月来水量，丰水年情况

qwan=[598.21,699.12,564.58,667.52,195.60,195.60,195.60,946.89,3235.13,1505.35,951.89,730.58];

%qwan=[598.21,699.12,564.58,667.52,195.60,195.60,195.60,946.89,3235.13,1505.35,951.89,730.58]; %丰水年

%qwan=[167.16,176.81,378.46,317.67,286.8,286.8,1061.56,1312.18,918.52,768.25,500.99,250.13];%平水年

%qwan=[332.3,295.18,664.44,365.42,26.15,45.32,122.58,950.26,622.13,141.65,236.98,334.18];%枯水年

%%%%%%%%%%%%%%%万龙区间每月来水，丰水年情况

qwl=[0,0,0,0,0,0,0,0,0,0,0,0];  %上下游水库区间流量w代表万家寨，l代表龙口

%%%%%%%%%%%%每月小时数

t(1)=31*24; t(3)= t(1); t(5)= t(1);

t(5)= t(1); t(7)= t(1); t(8)= t(1);

t(10)= t(1); t(12)= t(1);

t(2)=28*24;

t(4)=30*24; t(6)= t(4);t(9)= t(4);t(11)= t(4);

%% 依据水库上下游水位库容曲线关系，插值求上下游水位，计算发电量

for j=1:12        %%%%%%%万家寨水位推求库容

if 952<=x(j)<955

vw(j)=(4.508+(x(j)-952)/(955-952)*(4.548-4.508))*10^8;  %vw(j)：上游水位对应的库容

elseif 955<=x(j)<957

vw(j)=(4.548+(x(j)-955)/(957-955)*(4.739-4.548))*10^8;

elseif 957<=x(j)<960

vw(j)=(4.739+(x(j)-957)/(960-957)*(4.936-4.739))*10^8;

elseif 960<=x(j)<965

vw(j)=(4.936+(x(j)-960)/(965-960)*(5.546-4.936))*10^8;

elseif 965<=x(j)<970

vw(j)=(5.546+(x(j)-965)/(970-965)*(7.756-6.563))*10^8;

elseif 975<=x(j)<977.5

vw(j)=(7.756+(x(j)-975)/(977.5-975)*(8.35-7.756))*10^8;

elseif 977.5<=x(j)<980

vw(j)=(8.35+(x(j)-977.5)/(980-977.5)*(8.962-8.35))*10^8;

end

end

for j=1:11    %%%%%%%%%%%%%%%万家寨下泄流量推水位

qxw(j)=qwan(j)-(vw(j+1)-vw(j))/(t(j)*3600); % qxw(j)为出库流量

end

qxw(12)=qwan(12)-(vw(1)-vw(12))/(t(12)*3600);  %vw(13)=vw(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%万家寨下泄流量水位关系

for j=1:12

if 0<=qxw(j)<85

hxw(j)=898+(qxw(j)-0)/(85-0)*(899-898);  %下游流量对应的水位

elseif 85<=qxw(j)<174

hxw(j)=899+(qxw(j)-85)/(174-85)*(900-899);

elseif 174<=qxw(i)<286

hxw(j)=900+(qxw(j)-174)/(286-174)*(901-900);

elseif 286<=qxw(j)<571

hxw(j)=901+(qxw(j)-286)/(571-2860)*(902-901);

elseif 571<=qxw(j)<878

hxw(j)=902+(qxw(j)-571)/(878-571)*(903-902);

elseif 878<=qxw(j)<1320

hxw(j)=903+(qxw(j)-878)/(1320-878)*(904-903);

elseif 1320<=qxw(j)<1860

hxw(j)=904+(qxw(j)-1320)/(1860-1320)*(905-904);

elseif 1860<=qxw(j)<2480

hxw(j)=905+(qxw(j)-1860)/(2480-1860)*(906-905);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%万家寨上下游水位差

for j=1:11

hw(j)=(x(j)+x(j+1))/2-hxw(j); %上游平均水位-下游水位  %%%%%% 此处修改

end

hw(12)=(x(12)+x(1))/2-hxw(12);

%%%%%%%%%%%%%%%%%%%%%%%万家寨发电量，暂取A=8.0

for j=1:12

nw(j)=8.0*qxw(j)*hw(j);  %出力系数取为8.0，可以改进，求每个月出力值

ew(j)=8.0*qxw(j)*hw(j)*t(j); %求每个月电量值

end

sw=ew(1);

for j=2:12

sw=sw+ew(j);  %万家寨电量累计 求和

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%龙口

%%%%%%%%%%%%%%%%龙口库容曲线

for j=13:24

if 888<=x(j)<889

vl(j)=(0.0246+(x(j)-888)/(889-888)*(0.0448-0.0246))*10^8;%vl龙口

elseif 889<=x(zqun,j)<890

vl(j)=(0.0448+(x(j)-889)/(890-889)*(0.0713-0.0448))*10^8;

elseif 890<=x(j)<891

vl(j)=(0.0713+(x(j)-890)/(891-890)*(0.1089-0.0713))*10^8;

elseif 891<=x(j)<892

vl(j)=(0.1089+(x(j)-891)/(892-891)*(0.1715-0.1089))*10^8;

elseif 892<=x(j)<893

vl(j)=(0.1715+(x(j)-892)/(893-892)*(0.2502-0.1715))*10^8;

elseif 893<=x(j)<894

vl(j)=(0.2502+(x(j)-893)/(894-893)*(0.3373-0.2502))*10^8;

elseif 894<=x(j)<895

vl(j)=(0.3373+(x(j)-894)/(895-894)*(0.4314-0.3373))*10^8;

elseif 895<=x(j)<896

vl(j)=(0.4314+(x(j)-895)/(896-895)*(0.53-0.4314))*10^8;

elseif 896<=x(j)<897

vl(j)=(0.53+(x(j)-896)/(897-896)*(0.6317-0.53))*10^8;

elseif 897<=x(j)<898

vl(j)=(0.6317+(x(j)-897)/(898-897)*(0.7375-0.6317))*10^8;

end

end

for j=13:23     %%%%%%%%%%%%%%%%%%%%%%龙口下泄流量

qxl(j)=qxw(j-12)+qwl(j-12)-(vl(j+1)-vl(j))/(t(j-12)*3600);

end

qxl(24)=qxw(12)+qwl(12)-(vl(13)-vl(24))/(t(12)*3600);  %vl(25)=vl(13)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%龙口下泄流量、下游水位关系

for j=13:24

if 150<=qxl(j)<440

hxl(j)=861+(qxl(j)-150)/(440-150)*(861.5-861);

elseif 440<=qxl(j)<900

hxl(j)=861.5+(qxl(j)-440)/(900-440)*(862-861.5);

elseif 900<=qxl(j)<2180

hxl(j)=862+(qxl(j)-900)/(2180-900)*(863-862);

elseif 2180<=qxl(j)<3780

hxl(j)=863+(qxl(j)-2180)/(3780-2180)*(864-863);

elseif 3780<=qxl(j)<5820

hxl(j)=864+(qxl(j)-3780)/(5820-3780)*(865-864);

elseif 5820<=qxl(j)<8220

hxl(j)=865+(qxl(j)-5820)/(8220-5820)*(866-865);

elseif  8220<=qxl(j)<10800

hxl(j)=866+(qxl(j)-8220)/(10800-8220)*(867-866);

end

end

for j=13:23  %%%%%%%%%%%%%%%%%%%%%龙口上下游水位差

hl(j)=(x(j)+x(j+1))/2-hxl(j);

end

hl(24)=(x(24)+x(13))/2-hxl(24);

for j=13:24  %%%%%%%%%%%%%%%%%%%%%%%%%龙口发电量

nl(j)=8.6*qxl(j)*hl(j);  %出力系数取为8.6

el(j)=8.6*qxl(j)*hl(j)*t(j-12);

end

sl=el(13);

for j=14:24

sl=sl+el(j);  %龙口发电量累计

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%梯级电站的总体年发电量

E=(sw+sl); %两个水库总发电量

untitled.jpg

(16.04 KB, 下载次数: 3)

2015-6-29 11:13 上传

函数最优解