基于matpower潮流计算的风电并网粒子群无功优化实例

%---------------------------------说明----------------------------------------------- %--------------------------------风电并网无功补偿优化--------------------------------- %问题描述：对象为风电接入的IEEE33节点配电系统，已知风电分别接入在10节点（pw1）和17节点（pw2），采用粒子 %群优化算法求解无功补偿装置的最优补偿无功功率，使得系统的网损最小，潮流计算通过调用matpower计算。 %目标函数：确定无功补偿装置接入系统的最优无功补偿注入功率使得系统的运行网损最小（程序中的注释有详细的说明） %约束：风电出力上下限，关于粒子群迭代过程中，粒子位置越限处理注释中有说明。 %时间步长：1h；优化时域为一天（24小时） %------------------------------------------------------------------------------------ clear; close all clc; % F[T,1]为IEEE33节点的各时间段的各节点的负荷变化系数。 F[T,2]各时间段PW1的有功出力， F[T,3]为各时间段PW2的有功出力。 F=[1.1479 0.4241 0.50892 1.1382 0.4048 0.48576 1.1143 0.4805 0.5766 1.1062 0.4919 0.59028 1.1591 0.4248 0.50976 1.1869 0.3805 0.4566 1.2064 0.3972 0.47664 1.2301 0.4553 0.54636 1.2775 0.4188 0.50256 1.2901 0.382 0.4584 1.2913 0.3326 0.39912 1.2802 0.3413 0.40956 1.2175 0.3822 0.45864 1.2398 0.3033 0.36396 1.244 0.3352 0.40224 1.2413 0.3791 0.45492 1.283 0.4289 0.51468 1.301 0.4268 0.51216 1.4042 0.4255 0.5106 1.3374 0.461 0.5532 1.2984 0.4465 0.5358 1.2648 0.4096 0.49152 1.2272 0.4289 0.51468 1.1952 0.4005 0.4806 ]; %潮流计算参数 Sb=10; %基准功率 Ub=12.66; %基准电压 Zb=Ub^2/Sb; %基准阻抗 tic; %开始计时 for T=1:24 f=F(T,:); x=[0 0 0]; %无无功补偿时 % [U,ws,S0,S1]=flow(x,f); mpc=case33bw; %11节点-pw1 mpc.bus(11,3)=(0.045*f(1)-f(2))/Sb; mpc.bus(11,4)=(0.030*f(1)-f(2)*tan(acos(0.9)))/Sb; %17节点-pw2 mpc.bus(18,3)=(0.09*f(1)-f(3))/Sb; mpc.bus(18,4)=(0.040*f(1)-f(3)*tan(acos(0.9)))/Sb; %无功补偿装置 mpc.bus(22,4)=(0.040*f(1)-x(1))/Sb; %22节点 mpc.bus(25,4)=(0.200*f(1)-x(2))/Sb; %25节点 mpc.bus(33,4)=(0.040*f(1)-x(3))/Sb; %33节点 result_pf(T)=runpf(mpc); %将潮流结果存入result_pf结构体中 if result_pf(T).success==1 %判断潮流是否收敛 csv=result_pf(T).bus(:,8); % 优化前的各节点的电压 csdy=sum(abs(result_pf(T).bus(:,8)-ones(33,1))); % 优化前的各节点的电压偏差之和 csws=sum(abs(result_pf(T).branch(:,14)+result_pf(T).branch(:,16))); %优化前的有功网损 csgp=result_pf(T).gen(1,2); %优化前上级电网输入的有功功率 csgq=result_pf(T).gen(1,3); %优化前上级电网输入的无功功率 end %PSO参数设置 popsize=40; %粒子数 MAXITER=100;%迭代次数 dimension=3;%维数 w_max = 0.9;%最大惯性因子 w_min = 0.4; c1 = 2; %粒子个体学习因子,加速度常数 c2 = 2; %粒子社会学习因子，加速度常数 q1min=0; %接入节点21的无功补偿装置无功出力下限 q1max=0.8; %接入节点21的无功补偿装置无功出力上限 q2min=0; %接入节点24的无功补偿装置无功出力下限 q2max=0.8; %接入节点24的无功补偿装置的无功出力上限 q3min=0; %接入节点32的无功补偿装置无功出力下限 q3max=0.8; %接入节点32的无功补偿装置的无功出力上限 v=2*(rand(popsize,dimension)); %初始化粒子飞翔速度 Q1=(q1max-q1min)*rand(popsize,1); %初始化接入节点21无功补偿装置的无功出力 Q2=(q2max-q2min)*rand(popsize,1); %初始化接入节点24无功补偿装置的无功出力 Q3=(q3max-q3min)*rand(popsize,1); %初始化接入节点32无功补偿装置的无功出力 x=[Q1 Q2 Q3]; pbest=x; %初始化粒子群粒子局部最优解 gbest=zeros(1,dimension); %定义全局最优粒子 Ulim=zeros(32,1); %初始化支路电压 %% 随机初始化后，计算每个粒子的适应值，当前粒子的适应值作为初始化局部最优解，并找出全局最优解 for i=1:popsize % [U,ws,S0,S1]=flow(x(i,:),f); mpc.bus(22,4)=(0.040*f(1)-x(i,1))/Sb; %22节点 mpc.bus(25,4)=(0.200*f(1)-x(i,2))/Sb; %25节点 mpc.bus(33,4)=(0.040*f(1)-x(i,3))/Sb; %33节点 result_pf(T)=runpf(mpc); a=10000; %罚函数值 if result_pf(T).success==1 %判断潮流是否收敛，设置罚函数 U=result_pf(T).bus(:,8); % 各节点的电压 U_sum=sum(abs(result_pf(T).bus(:,8)-ones(33,1))); % 各节点的电压偏差之和 ws=sum(abs(result_pf(T).branch(:,14)+result_pf(T).branch(:,16))); %有功网损 S0=result_pf(T).gen(1,2); %上级电网输入的有功功率 S1=result_pf(T).gen(1,3); %上级电网输入的无功功率 f_pbest(i)=ws/csws; %初始状态时各粒子的适应度函数 else f_pbest(i)=a;%初始状态时各粒子的适应度函数 end end g=find(f_pbest==min(f_pbest(1:popsize)));%找出f_pbest中等于min(f_pbest(1:popsize))元素的下标 if length(g)~=1 g=g(1); end gbest=pbest(g,:); %得到初始状态时的全局最优粒子 f_gbest=f_pbest(g); %初始状态时种群的最小目标函数值 MINIUM=f_gbest; %% 迭代，更新 for t=1:MAXITER disp('MINIUM='); disp(MINIUM); disp('迭代次数t='); disp(t); w = w_max-(w_max-w_min)*t/MAXITER;%惯性权重自适应调整 for i=1:popsize %--------------------------------更新粒子位置--------------------------------- v(i,:) = w.*v(i,:)+c2.*rand.*(gbest(1:dimension)-x(i,1:dimension))+c1.*rand.*(pbest(i,1:dimension)-x(i,1:dimension));% 更新粒子速度 if abs(v(i,1))>0.003 %最大速度设置，粒子的范围宽度 v(i,1) = sign(v(i,1))*0.003; end if abs(v(i,2))>0.003 %最大速度设置，粒子的范围宽度 v(i,2) = sign(v(i,2))*0.003; end if abs(v(i,3))>0.003 %最大速度设置，粒子的范围宽度 v(i,3) = sign(v(i,3))*0.003; end x(i,1:dimension) = x(i,1:dimension)+v(i,1:dimension);%更新后的粒子 %----------------------------粒子越限处理----------------------------------- if x(i,1)>q1max x(i,1)=q1max; end if x(i,1)<0 x(i,1)=0; end if x(i,2)>q2max x(i,2)=q2max; end if x(i,2)<0 x(i,2)=0; end if x(i,3)>q3max x(i,3)=q3max; end if x(i,3)<0 x(i,3)=0; end %-----------------计算更新后的粒子的适应值---------------------------- mpc.bus(22,4)=(0.040*f(1)-x(i,1))/Sb; %22节点 mpc.bus(25,4)=(0.200*f(1)-x(i,2))/Sb; %25节点 mpc.bus(33,4)=(0.040*f(1)-x(i,3))/Sb; %33节点 result_pf(T)=runpf(mpc); a=10000; %罚函数值 if result_pf(T).success==1 %判断潮流是否收敛，设置罚函数 U=result_pf(T).bus(:,8); % 各节点的电压 U_sum=sum(abs(result_pf(T).bus(:,8)-ones(33,1))); % 各节点的电压偏差之和 ws=sum(abs(result_pf(T).branch(:,14)+result_pf(T).branch(:,16))); %有功网损 S0=result_pf(T).gen(1,2); %上级电网输入的有功功率 S1=result_pf(T).gen(1,3); %上级电网输入的无功功率 f_x(i)=ws/csws; %初始状态时各粒子的适应度函数 else f_x(i)=a;%初始状态时各粒子的适应度函数 end %-----------------------------与历史全局最优解比较---------------------------- if f_x(i)<f_pbest(i) pbest(i,:)=x(i,:);%若粒子的目标函数小于更新前的局部最优值，则局部最优值就为这个粒子（局部最优解对应粒子的最优解，即最小值） f_pbest(i)=f_x(i); end if f_pbest(i)<f_gbest gbest=pbest(i,:);%若更新后的局部最优值小于全局最优值，则全局最优值就是局部最优值 f_gbest=f_pbest(i); end fitness=f_gbest; end MINIUM = f_gbest; yy(t)=MINIUM; end G(T,:)=gbest; %得到各