IEEE33节点 基于遗传算的无功优化 MATLAB代码

控制变量为：电容组投切次数，范围是0~10。优化目标是网络损耗最小。代码如下，欢迎大家理性探讨，写得不好的地方，也请多多包涵。本代码远远不够完善，存在很多瑕疵，欢迎大佬指出我的错误。

// %%%%%采用粒子群算法，控制变量为ieee33节点中9和32处的电容器组投切数量，优化目标为网损最小
clc
clear;
clear all;
close all;
%% GA参数
maxgen=30;%迭代次数
sizepop=40;%种群规模
popmax=[10];%电容器组投入上限
popmin=[0];%电容器组投入下限
dim=[2];%有两个节点进行了电容器补偿
shuju=loadcase('ieee33');
fitness=0;%初始的适应度为0
D=zeros(maxgen,4);%用来记录每代的最优解,平均值，最差解,最优解是否为可行解
fv=inf;%初始最优值的适应度函数为正无穷
%% 种群初始化
pc=[zeros(sizepop,3)];
for i=1:sizepop
  pop(i,1:dim(1))=round(popmax(1)*rand(1,dim(1)));%初始化电容器投入组数，取整  
  shuju.bus([9,32],6)=pop(i,1:dim(1))*0.1;%%%%%%%%%%%%此处有疑问%%%%%%%%%% 是因为第六列代表了并联电容器的容量吗
  [MVAbase,bus,gen,branch,success,et]=runpf(shuju,mpoption('OUT_ALL',0,'VERBOSE',0));%不显示迭代过程
  if success==1 %判断是否能解潮流
      pc(i,1:dim(1))=pop(i,1:dim(1));
      fitness(i)=sum(branch(:,14)+branch(:,16));%计算网损大小
      pc(i,3)=fitness(i);
      pc(i,4)=1;%第一二列表示数目，第三列表示适应度，第四列表示是否为可行解
  else
      pc(i,1:dim(1))=pop;
      pc(i,3)=inf;
      pc(i,4)=0;
  end
end
[best, index]=min(fitness);
gtsite=pop;%个体位置
gtfit=fitness;%个体适应度
popsite=pop(index,:);%全局最佳位置
popfit=fitness(index); %全局最佳适应度
trace(1)=popfit;%记录每一代的最佳适应度
trace_site(1,:)=popsite;%记录每一代的最佳位置
D=zeros(maxgen,3);%记录每代最优解、最差解、是否为可行解
%% 进化

for i=1:maxgen
    %个体选择——锦标赛选择法
    M=zeros(sizepop,dim(1)+2);%用来存储优胜者的中间矩阵
    for j=1:sizepop
        A=randperm(sizepop,2);%随机获取1~30范围的两个数
        if pc(A(1),3)<=pc(A(2),3)%适应度小的的保留下来
            M(j,:)=pc(A(1),:);
        else
            M(j,:)=pc(A(2),:);%
        end        
    end
    pc=M;%将优胜者保存下来
    
    %交叉
     proc=0.8;%交叉概率——从所有个体中取80%进行交叉操作
     NC=floor(proc*sizepop);%交叉操作的个体数
     NC=NC-mod(NC,2);%保证是偶数
     tempCross=randperm(sizepop,NC);%执行交叉操作的个体,1*22
     pointCross=randperm(dim(1),1);%随机生成交叉点,1*1
     alpha=rand(pointCross,1); %2*1
     alpha=alpha';%转置
     for k=1:0.5*NC
         Cross1=alpha.*pc(2*k-1,1:pointCross)+(1-alpha).*pc(2*k,1:pointCross);
         Cross2=alpha.*pc(2*k,1:pointCross)+(1-alpha).*pc(2*k-1,1:pointCross);
         pc(2*k-1,1:pointCross)=Cross1;
         pc(2*k,1:pointCross)=Cross2;
         %%%交叉是否需要重新编码 %%%%%%%%
     end  
     
     %变异——每次取一定量个体进行变异
     prom=0.1;
     NM=floor(prom*sizepop)+1;%变异个体数
     idM=randperm(sizepop,NM);%变异个体编号
     pointMu=randperm(dim(1),1);%随机生成变异位置
     tempM=rand(1);%变异方向随机
     for j=idM
         pc(j,pointMu)=randperm(popmax(1),1);%随机获得一个值作为新的一代
%          if tempM>0.5
%              pc(j,pointMu)=pc(j,pointMu)+0.1*popmax(1);
%          else
%              pc(j,pointMu)=pc(j,pointMu)-0.1*popmax(1);
%          end
     end
     pc(:,1:dim(1))=round(pc(:,1:dim(1)));
     pop=(pc(:,1:dim(1)));
        

for j=1:sizepop
    %计算适应度
    shuju.bus([9 32],6) =pop(j,1:dim(1))*0.1; %改电容器容量
    [basemva, bus, gen, branch, success, et] =runpf(shuju, mpoption('OUT_ALL',0,'VERBOSE',0));
    if success==1
        fitness=sum(branch(:,14)+branch(:,16));%计算出网损大小
        pc(j,3)=fitness;
        pc(j,4)=1;
    else
        fitness=inf;
        pc(j,3)=fitness;
        pc(j,4)=0;
    end
    % fitness=f1(shuju);
    %个体适应度及位置更新
    if fitness<gtfit(j)
    gtfit(j)=fitness;
    gtsite(j,:)=pop(j,:);
    end
    %全局适应度及位置更新
    if fitness<popfit
    popfit=fitness;
    popsite=pop(j,:);
    end
end
[Q,IX]=sort(pc,1);
 D(i,2)=Q(sizepop,3);%sort默认升序，所以这是适应函数最大值
 D(i,1)=Q(1,3);%适应度函数最小值
 D(i,3)=pc(IX(1,3),4);%是否为可行解 
trace(i+1)=popfit;%记录每次迭代值
trace_site(i+1,:)=popsite;
if  D(i,1)<fv&&D(i,3)==1
    fv=D(i,1);
end
end

figure(1)
plot(trace)
xlabel('迭代次数')
ylabel('适应度函数')%%%%图有问题，哪个才是真实的？


k=2;
N(1,:)=[trace(1),trace(1),1];
for i=2:maxgen
    if D(i,3)==1
         N(k,:)=D(i,:);
         k=k+1;
    end
end
figure(2)
plot(N(:,1));
legend('最优值');