这是我的代码，帮我分析下，我的粒子是什么
%MOPSO:63、311;example:42
% ----------------------------------------------------------------------- %
clear
clc
mpc=loadcase(‘case33bw_60’);
n_Bus=mpc.bus(end,1);
Load=[ 0.8, 0.805, 0.81, 0.818, 0.83, 0.91, 0.95, 1, 1.05, 1.08, 1.03, 1.01, 0.95, 0.9, 0.905, 0.91, 0.93, 1.02, 1.1, 1.18, 1.2, 1.08, 0.9, 0.84]/ 1.4;
% 多目标函数
MultiObjFnc = ‘mubiao’;
switch MultiObjFnc
case ‘mubiao’ % 总的THDu
MultiObj.fun = @(x) m(x);
%x是充电桩节点的Iv,nm,n个节点，每个节点是m次；m=8,n=10
%按照每个节点的1-8次电流来划分
%总数
MultiObj.nVar = 108;
%这里再加入判断有几个反向接入的，用概率再取整；
g=0.2;
n=80*g;%g为每时刻的反向接入的概率
%前几个的最小/最大
MultiObj.var_min(1:n) = -15;
MultiObj.var_max(1:n) = 0;
MultiObj.var_min(n+1:80) = 0;
MultiObj.var_max(n+1:80) = 15;
end
%% MOPSO参数
params.Np = 200; % 种群规模
% 增加种群多样性的参数设置
params.W = 0.7; % 提高惯性权重
params.C1 = 1.0; % 降低个体认知因子
params.C2 = 1.0; % 降低社会认知因子
params.maxvel = 5; % 增大最大速度限制
params.u_mut = 0.5; % 增加变异概率
params.Nr = 200; % 存储库大小
params.maxgen = 25; % 最大迭代次数
%params.W = 0.4; % 惯性权重
%params.C1 = 1.4; % Individual confidence factor
%params.C2 = 1.4; % Swarm confidence factor
params.ngrid = 10; % Number of grids in each dimension
%params.maxvel = 10; % Maxmium vel in percentage
%params.u_mut = 0.8; % Uniform mutation percentage
function f= m(x)
%%这个函数就是给一组粒子，返回一个sumTHDU!
%%
n=size(x,1);%粒子数
f=zeros(n);
for l=1:n %l的意思是每一种粒子情况下的一个sumTHDU,最终会合在一个f中，输出给example的REP中
%80个谐波电流值用x(l,j)的形式带入，表示第l个粒子情况下的第j个谐波电流的大小，10*8的一行排布
%以上是构建粒子群的循环代码，以下是在一个粒子值下的运算结果（得到一个sum）
%%
ac_data = case33bw_60; %这里未使用loadcase %case33bw case_ieee34 ac_case8(阻抗较小)
ac_baseMVA = ac_data.baseMVA;
ac_bus = ac_data.bus;
ac_branch = ac_data.branch;
ac_gen = ac_data.gen;
ac_dg = ac_data.dg; %%电网数据
source1 = ac_data.source; %谐波源
charge1=ac_data.charge;
%充电桩需要再单独在33中去设置
ac_bn = find(ac_bus(:, 2)==3); %平衡节点
ac_pv = find(ac_bus(:, 2)==2); %PV节点
ac_pq = find(ac_bus(:, 2)==1); %PQ节点
UB = ac_bus(1,10);
ZB =UB^2/ac_baseMVA;
IB = ac_baseMVA/sqrt(3)/UB;
m = length (ac_pq);
% 初始化变量：保存各节点各次谐波注入电流幅值
harmonic_num = length(source1(1, 😃)-1;
bus_num = length(ac_bus(:, 1));
I_bus_h_total_abs = zeros(harmonic_num, bus_num); % harmonic_num x bus_num
%% 基于牛拉法的基波潮流计算
%基波下，充电桩的参数设定，与x包含的80值无关，但与正反向充电时的功率有关，这块的思路
Ybus = createYbus(ac_baseMVA, ac_bus, ac_branch);
Pacs = -ac_bus(:,3); %有功负荷消耗，看做负值
Pacs(ac_gen(:,1)) = ac_gen(:,2) - ac_bus(ac_gen(:,1),3); %发电机节点的有功功率
% 删不删去dg,充电桩的有功功率要不要算，有功注入与无功注入
Pacs(ac_dg(:,1)) = Pacs(ac_dg(:,1))+ ac_dg(:,2) ; %接入dg电源用功注入
Qacs = -ac_bus(:,4);
Qacs(ac_gen(:,1)) = ac_gen(:,3) - ac_bus(ac_gen(:,1),4); %发电机节点的无功功率
Qacs(ac_dg(:,1)) = Qacs(ac_dg(:,1))+ ac_dg(:,3) ; %接入dg电源无功注入
Pacs = Pacs/ac_baseMVA; Qacs = Qacs /ac_baseMVA; %标幺化,bus不用标幺化
kac = 0;%迭代次数
%计算节点功率Pi，Qi
busNum = length(ac_bus(:, 1));
Bus_V = ac_bus(:, [1,8:9]); %仅用于形成雅克比矩阵
Bus_V(ac_gen(:,1),2) = ac_gen(:,6); %平衡节点电压
Pii=zeros(busNum, 2); %所有节点(包括平衡节点)的Pi
Qii=zeros(busNum, 2);
Pii(:,1)=Bus_V(:, 1); %没有特殊的含义，复制编号
Qii(:,1)=Bus_V(:, 1);
Y = Ybus([ac_pq;ac_pv],[ac_pq;ac_pv]); %除去平衡节点后的Y阵，用于雅克比矩阵，直接去掉该行列
Pacs = Pacs([ac_pq;ac_pv]😅;
Qacs = Qacs([ac_pq;ac_pv]😅;
while true
for i=1:busNum for j=1:busNum detal_ij=(Bus_V(i,3)-Bus_V(j,3))*pi/180; %角度转弧度 Pii(i, 2)=Pii(i, 2)+Bus_V(j, 2)*(real(Ybus(i, j))*cos(detal_ij)+imag(Ybus(i, j))*sin(detal_ij)); %Pi从2开始 Qii(i, 2)=Qii(i, 2)+Bus_V(j, 2)*(real(Ybus(i, j))*sin(detal_ij)-imag(Ybus(i, j))*cos(detal_ij)); end Pii(i, 2) = Bus_V(i, 2)*Pii(i, 2); Qii(i, 2) = Bus_V(i, 2)*Qii(i, 2); end %重新编号1~m为PQ节点，m+1~n-1+m为PV节点 Pi = Pii([ac_pq;ac_pv], :); %除去平衡节点且重置节点位置 Qi = Qii([ac_pq;ac_pv], :); %除去平衡节点且重置节点位置 V = Bus_V([ac_pq;ac_pv], :); Uac_m = V(:, 2); Uac_a = V(:, 3); Uacs = Uac_m .* exp(1j * (Uac_a*pi/180)); dP = Pacs - Pi(:,2); dQ = Qacs - Qi(:,2);
% dP = Pacs - real(Uacs.conj(YUacs));%公式不能用
% dQ = Qacs - imag(Uacs.conj(YUacs));
if max(abs(dP)) <= 1e-9 && max(abs(dQ(1:m))) <= 1e-9 %收敛判据 break end [J,H,N,K,L] = Jacobi(V, Y, ac_pq ,ac_pv,Pi,Qi); %返回雅克比矩阵各子矩阵 Ud2 = diag(V(:,2)); dPQ = [dP ; dQ(1:m)]; ang_U = -J^-1*dPQ; dang = ang_U(1:busNum-1); %弧度结果 dU = diag(Uac_m(1:m))*ang_U(busNum:end); %数据修正 Bus_V(V(1:m,1),2) = Bus_V(V(1:m,1),2)+dU; %修正原编号下PQ节点电压幅值和相角 Bus_V(V(:,1),3) = Bus_V(V(:,1),3)+dang*180/pi; %弧度转角度,空载情况角度增加 Pii(:,2) =0; Qii(:,2) = 0; kac = kac+1;
end
%% 平衡节点功率计算Sn
%默认只有一个平衡节点
Ui = Bus_V(:, 2).cos(Bus_V(:, 3)pi/180)+1iBus_V(:, 2).sin(Bus_V(:, 3)pi/180); %电压复数形式
Sn=0;
for i=1:busNum
Sn = Sn+ Ui(ac_bn)(conj(Ybus(i, ac_bn))(conj(Ui(i)))); %Uisum(conj(Yij)*conj(Uj))
end
%方法一：计算支路电流
branch_num = length(ac_branch(:, 11));
I_branch = zeros(branch_num, 3);
S_branch = zeros(branch_num, 5);
I_branch(:, 1:2) = ac_branch(:, 1:2);
S_branch(:, 1:2) =ac_branch(:, 1:2);
from_node = ac_branch(:, 1); to_node = ac_branch(:, 2);
z_branch = (ac_branch(:, 3)+1iac_branch(:, 4))/ZB;
I_branch(:, 3) = (Ui(from_node)-Ui(to_node))./(z_branch); %支路流过的电流
for k = 1:length(I_branch(:,1)) %考虑电容注入电流
I_branch(k, 3) = I_branch(k, 3) +Ui(from_node(k))(1iac_branch(k,5)ZB)/2;
I_branch(k, 3) = I_branch(k, 3) +Ui(to_node(k))(1iac_branch(k,5)*ZB)/2;
end
S_branch(:, 3) = Ui(from_node).*conj(I_branch(:, 3))*ac_baseMVA; %支路ij流过的容量
S_branch(:, 4) = Ui(to_node).*conj(-I_branch(:, 3))*ac_baseMVA; %支路ji流过的容量
S_branch(:, 5) =real(S_branch(:, 3)+S_branch(:, 4)); %基波潮流有功损耗
ac_loss = sum(S_branch(:, 5));
%计算节点注入电流
bus_num = length(ac_bus(:, 1));
Iii = zeros(bus_num, 2);
Iii(:, 1) = ac_bus(:, 1);%标点
Sii = Pii(:,2)+1iQii(:, 2);
Iii(:, 2) = conj(Sii)./conj(Ui); %注意共轭的取值
Ii = YbusUi; %Ii = Yi1U1+Yi2U2+Yi3*U3+……
%% 诺顿谐波潮流计算
harmonic_U_h = zeros(bus_num, 3);%3列，1节点、2谐波幅值、3相角
harmonic_U_h(:, 1)=ac_bus(:, 1);%第一列是节点数
%harmonic_U_h算的是全部33节点的
source_bus = source1(:, 1);%确定谐波源的节点数作为第一列
%单独谐波源的节点编号
%Charge_bus = [5,10,13,16,20,24,11,28,31,17]; % 10个特定节点的编号,作为充电桩
charge_bus=charge1(:,1);
%单独充电桩的节点编号
%但是，应该在33节点中设置
harmonic_U_h(source_bus, 2) = 0.05; %谐波源的谐波幅值设为0.05
%全体先设为了0.05吗？
harmonic_U_h(:, 3) = 0; %谐波相角设为0
harmonic_num = length(source1(1, 😃)-1; %此处考虑5,7,11,13,17次谐波，17也全为0
load_power = (ac_bus(:, 3)+1i*ac_bus(:, 4)); %负荷功率MVA
Zload = zeros(bus_num, harmonic_num);
V_bus_h_total_abs = zeros(harmonic_num , bus_num);
I_branch_h_total_abs = zeros(harmonic_num , branch_num);
loss = zeros(branch_num ,harmonic_num );
Pi_fam = -Pii(:,2)ac_baseMVA; Qi_fam = -Qii(:,2)ac_baseMVA;Sii_fam = Siiac_baseMVA;%采用有名值计算
source_R1 = Pi_fam(source_bus).(Bus_V(source_bus, 2)UB).^2./(abs(Sii_fam(source_bus))).^2;%阻抗串联模型
source_X1 = Qi_fam(source_bus).(Bus_V(source_bus, 2)*UB).^2./(abs(Sii_fam(source_bus))).^2;
for i = 1:harmonic_num %计算各次的 %谐波源负荷模型被代替
h = 2*i+1; %次数
k=i;
I_h = zeros(bus_num, 1); %所有节点电流向量,表示一个特定次数的所有电流值（列向量）
charge1(:, i+1)=0;
if charge1(:, i+1)==0 %存在一个充电桩发出该次谐波
%这是在假设充电桩也在33以后的形式
Charge_Ih = zeros(bus_num, 1); % 初始化已知谐波电流矩阵（复数，标幺值）
for j=1:10
index = (j-1)8 + i; % 第j个节点的第i次谐波在x中的位置
Charge_Ih(charge_bus(j)) = x(l, index); % charge_bus(j)是第j个充电桩的节点编号
end
I_h = Charge_Ih;
%上述是对已知量的先赋值给I_h,对于剩下的谐波源及其他节点再计算他们的谐波电流
if source1(:, i+1)~=0 %存在谐波源发出该次谐波
I_h0= -Iii(source1(:,1), 2).source1(:, i+1); %用节点注入电流约为0.0075时谐波电压畸变率2%
%什么意思 基准的意思
harmonic_Uhs = harmonic_U_h(:,2). exp(1j * (harmonic_U_h(:,3)pi/180));%谐波电压初始值
Zh_eq = source_R1+1ihsource_X1;
while true %解耦迭代过程
I_eq = harmonic_Uhs(source1(:,1))./Zh_eq; I_eq = I_eq/IB; %标幺化 I_h(source1(:,1)) = I_h0+I_eq; %发电机谐波模型 XG_h = 1i*h*ac_gen(:, 23)/ZB; %负序电抗0.2；XG=1i*h*0.2*sqrt(real(Sn)^2+imag(Sn)^2)； yG = 1./XG_h; Zload([ac_bn;ac_pv] , i) = XG_h; %负荷谐波阻抗模型 Num_L= find( load_power(ac_pq) ~=0); %pq节点中非零负荷的位置 RL = real(load_power(ac_pq(Num_L))).*(Bus_V(ac_pq(Num_L), 2)*UB).^2./(abs(load_power(ac_pq(Num_L)))).^2;%有名值计算 RL = RL/ZB; %33节点系统需要标幺化,下同 XL = imag(load_power(ac_pq(Num_L))).*(Bus_V(ac_pq(Num_L), 2)*UB).^2./(abs(load_power(ac_pq(Num_L)))).^2; XL = XL/ZB; Zload(ac_pq(Num_L) ,i) = RL+1i*h*XL; %Z=U^2*S/|S|^2 Zload(source1(:,1) ,i) = Zh_eq/ZB; %电源节点和负荷节点自导纳 yloads=zeros(size(Zload)); for j=1:bus_num if Zload(j, i) ~=0 yloads(j, i) = 1/Zload(j, i); end end
% yloads=zeros(bus_num,1); %各节点的负荷等值模型
% for j=1:bus_num
% if load_power(j) ~=0
% yloads(j, 1)=real(load_power(j))/Bus_V(j, 2)^2-1iimag(load_power(j))/(hBus_V(j, 2)^2); %yii = Pi/Vi^2-jQi/(hVi^2)？？？
% else
% yloads(j, 1) = 0;
% end
% end
%线路谐波模型
branch_h = ac_branch;%已经标幺化
branch_h(: ,3) = ac_branch(: ,3);
branch_h(: ,4) = hac_branch(: ,4);
branch_h(: ,5) = hac_branch(: ,5);
%形成网络谐波导纳矩阵
bus_h = ac_bus;%bus没有标幺化
bus_h(:, 5) = real(yloads(:, i))+bus_h(:, 5)*ZB; %导纳乘以阻抗基准值
bus_h(:, 6) = imag(yloads(:, i))+bus_h(:, 6)*ZB;
Ybus_h = createYbus(ac_baseMVA, bus_h, branch_h);
%由负荷等值模型和发电机等值模型组成
U_h = (Ybus_h^-1)*I_h;
%迭代收敛判据
dU_h = harmonic_Uhs -U_h;
if max(abs(dU_h)) <= 1e-7%收敛判据
break
end
harmonic_Uhs = U_h;
end
V_bus_h_total_abs (i,:) = abs(U_h); I_bus_h_total_abs(i, :) = abs(I_h); % 保存各节点注入电流幅值 z_branch_h = branch_h(:, 3)+1i*branch_h(:, 4); I_branch_h = (U_h(from_node)-U_h(to_node))./(z_branch_h); % for k = 1:length(I_branch_h) I_branch_h(k) = I_branch_h(k)+U_h(from_node(k))*(1i* branch_h(k,5))/2; I_branch_h(k) = I_branch_h(k)+U_h(to_node(k))*(1i* branch_h(k,5))/2; end I_branch_h_total_abs(i, :) = abs(I_branch_h); Sij_branch_h = U_h(from_node).*conj(I_branch_h); %支路ij流过的容量 Sji_branch_h = U_h(to_node).*conj(-I_branch_h); %支路ji流过的容量 loss(: ,i) =abs(real(Sij_branch_h+Sji_branch_h))*ac_baseMVA*10^6; %谐波潮流有功损耗
% for k= 1:branch_num
% yij = abs(1/(branch_h(k, 3)+1i*branch_h(k, 4)));
% I_branch_h_total_abs(i, k)=abs(U_h(branch_h(k ,1))-U_h(branch_h(k,2)))yij;%Iij = |Ui-Uj||yij|
% end
else %fprintf('******%dth harmonic is not considered******\n',h); end
end
end
p=size(abs(U_h));
THDU = sqrt(sum(V_bus_h_total_abs.^2, 1)) ./ Bus_V(:, 2)';%各个节点的谐波畸变率
sum_THDU = sum(THDU)/10000; % 所有节点的 THDU 总和
f11=sum_THDU;
f1(l,1)=f11;
% 计算电压偏差（例如：所有节点电压与额定电压的偏差的最大值）
nominal_voltage = 1.0; % 假设额定电压为1.0 p.u.
voltage_deviation = max(abs(Bus_V(:, 2) - nominal_voltage))/0.1;
f22 = voltage_deviation; % 第二个目标电压偏差最小化
f2(l,1)=f22;
end
f=[f1,f2];
end
% MOPSO
REP = MOPSO(params,MultiObj);%MOPSO的迭代次数决定了REP的行数，而M函数的值是一个迭代下x的所有值决定了一行中所有列的值
%所以，每次输出80个（一行）
%也即谐波计算要返回80个值
REP.pos(:,1:8);
REP.pos(:,9:16);%每个充电桩的8个谐波电流
REP.pos(:,17:24);
REP.pos(:,25:32);
REP.pos(:,33:40);
REP.pos(:,41:48);
REP.pos(:,49:56);
REP.pos(:,57:64);
REP.pos(:,65:72);
REP.pos(:,73:80);
h=REP.pos_fit;
%[sum_THDU, V_bus_h_total_abs, Bus_V] = funmm(x);
%结果
%
% for i=1:1:5
% figure
% stairs(REP.pos(i,9:32),‘ro-’)
% figure
% stairs(REP.pos(i,33:56),‘ro-’)
% end
%Figure(REP,12)
%最终的图是24小时的THD的变化图展示抑制谐波的效果；
%给目标，以各个的Iv初始值先运行一次，求出一次REP的矩阵（THDu在f1在M中，M在REP中），在进入MOPSO的循环，每次求一次REP
%REP是每次包含多个矩阵是因为M中原先两个函数，所以双维，现在只一个函数，则单维，但仍有多个矩阵
% 找到网损最小的解的索引
[~, minLossIndex] = min(REP.pos_fit(:, 2));
% 提取该解的决策变量和目标函数值
bestSolution = REP.pos(minLossIndex, 😃; % 决策变量
bestFitness = REP.pos_fit(minLossIndex, 😃; % 目标函数值
y=bestFitness;
h=REP.pos;
% 显示结果
%disp('对应的目标函数值（REP.pos_fit）😂;
%disp(bestFitness);
最新发布
扫一扫
4412
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
