上一篇:https://blog.csdn.net/m0_43401436/article/details/106564397
我自己的代码都差点认不出来了。。。
完整代码如下,安装matpower后可直接运行。注释写得比较清楚,结合上一篇应该能看明白。
clc
clear
%%
mpc=loadcase('case6ww');
bus=mpc.bus;%节点参数
branch=mpc.branch;%支路参数
gen=mpc.gen;%发电机参数
NG=size(gen,1);%发电机数目
%%
pgMax=1*gen(:,9);
pgMin=1*gen(:,10);%机组最大最小出力
gencost=mpc.gencost;%发电机花费参数
costa=gencost(:,5);
costb=gencost(:,6);
costc=gencost(:,7);
%% GA参数
size=30;%种群规模
NP=50;%进化多少代
cost=zeros(1,size);%每个个体的花费
fv=0;%初始最优值的适应函数值取0
D=zeros(NP,4);%用来记录每代的最优解,平均值,最差解,最优解是否为可行解
%% GA Process
%% 种群初始化
pg=(pgMax-pgMin).*rand(NG,size)+pgMin;
pg=[pg;zeros(3,size)];%再加3行,分别是4适应函数值、5是否为可行解、6燃料花费
%% 判断是否为可行解
for i=1:size
for j=1:NG
mpc.gen(j,2)=pg(j,i);
end
results=runpf(mpc);
if results.success==1
if (sum(results.gen(:,2)>=pgMin)==NG)&&(sum(results.gen(:,2)<=pgMax)==NG)%机组出力不越限
pg(5,i)=1;
pg(1:NG,i)=results.gen(1:NG,2);%实际计算结果赋值过来
cost(i)=sum((pg(1:NG,i).^2)'*costa+pg(1:NG,i)'*costb+ones(NG,1)'*costc);%燃料花费
pg(4,i)=1/cost(i);%适应函数值
pg(6,i)=cost(i);%燃料花费
else
pg(5,i)=0;
pg(4,i)=0;
end
else
pg(5,i)=0;
pg(4,i)=0;%非可行解适应函数值为0
end
end
%% 进化(包括:个体选择、交叉、变异)
for i=1:NP
%个体选择——锦标赛选择法
M=zeros(NG+3,size);%用来存储优胜者的中间矩阵
for j=1:size
A=randperm(size,2);%随机获取1~N范围的两个数
if pg(4,A(1))<=pg(4,A(2))%适应度大的保留下来
M(:,j)=pg(:,A(2));
else
M(:,j)=pg(:,A(1));
end
end
pg=M;%直接把优胜者种群作为新种群;或者也可以在竞争过程中把竞争失败的替换成胜利的
%交叉
proc=0.75;%交叉概率——从所有个体中取75%进行交叉操作
NC=floor(proc*size);%交叉操作的个体数
NC=NC-mod(NC,2);%保证是偶数
tempCross=randperm(size,NC);%执行交叉操作的个体
pointCross=randperm(NG,1);%随机生成交叉点
alpha=rand(pointCross,1);
for k=1:0.5*NC
pgCross1=alpha.*pg(1:pointCross,2*k-1)+(1-alpha).*pg(1:pointCross,2*k);
pgCross2=alpha.*pg(1:pointCross,2*k)+(1-alpha).*pg(1:pointCross,2*k-1);
pg(1:pointCross,2*k-1)=pgCross1;
pg(1:pointCross,2*k)=pgCross2;
end
%变异——每次取一定量个体进行变异
prom=0.05;
NM=floor(prom*size)+1;%变异个体数
idM=randperm(size,NM);%变异个体编号
pointMu=randperm(NG,1);%随机生成变异位置
tempM=rand(1);%变异方向随机
for j=idM
if tempM>0.5
pg(pointMu,j)=pg(pointMu,j)+0.01*pgMax(pointMu);
else
pg(pointMu,j)=pg(pointMu,j)-0.01*pgMax(pointMu);
end
end
%判断子代是否可行
for k=1:size
mpc.gen(1:NG,2)=pg(1:NG,k);
results=runpf(mpc);
if results.success==1
if (sum(results.gen(:,2)>=pgMin)==3)&&(sum(results.gen(:,2)<=pgMax)==3)
pg(5,k)=1;
pg(1:NG,k)=results.gen(1:NG,2);%实际计算结果赋值过来
cost(k)=sum((pg(1:NG,k).^2)'*costa+pg(1:NG,k)'*costb+ones(NG,1)'*costc);%燃料花费
pg(4,k)=1/cost(k);%适应函数值
pg(6,k)=cost(k);%燃料花费
else
pg(5,k)=0;
pg(4,k)=0;
end
else
pg(5,k)=0;
pg(4,k)=0;%非可行解适应函数值为0
end
end
%排序获取当前代的结果
[Q,IX]=sort(pg,2);%sort默认竖着排,1——竖着排;2——横着排
D(i,1)=Q(4,size);%sort默认升序,所以这是适应函数最大值
D(i,2)=mean(Q(4,:));%适应函数平均值
D(i,3)=Q(4,1);%适应函数最小值
D(i,4)=pg(5,IX(4,size));%IX(4,size)得到的是最优解在原矩阵E中的索引,5对应的是表示是否为可行解的那一列
D(i,5)=pg(6,IX(4,size));%最优解燃料花费
if D(i,1)>fv && D(i,4)==1
fv=D(i,1);%更新最优解目标函数值
costopt=1/fv;
xm=[pg(1:3,IX(4,size));pg(5,IX(4,size))];%最优解对应的自变量值,是否为可行解
end
end
%所有代进化完毕
%% 画图
k=1;
for i=1:NP %各代都画出来
if D(i,4)==1
N(k,:)=D(i,:);%如果第i代最优解是可行解,存储改最优解的全部信息,这里的N矩阵是未预先定义的,因为D里面存储的未必都是可行解
k=k+1;
end
end
plot(1./N(:,1),'r*');%这里取了个导数,变成了花费,看得更清楚一些,其实实际上直接用燃料花费作为适应度函数还方便一些,花费小的保留下来
legend('最优值');
输出的效果图如下,不是很好看,但大致意思有了,要用的话需要再自己调参数。
T h e e n d . The \ end. The end.