%MOPSO:63、311;example:42 %x = [ % x1_1, x1_2, x1_3, ..., x1_nVar; % 粒子1 % x2_1, x2_2, x2_3, ..., x2_nVar; % 粒子2 % ... % xN_1, xN_2, xN_3, ..., xN_nVar; % 粒子N %] % ----------------------------------------------------------------------- % 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,n*m,n个节点,每个节点是m次;m=8,n=10 %按照每个节点的1-8次电流来划分 %总数 MultiObj.nVar = 10*8; %这里再加入判断有几个反向接入的,用概率再取整; g=0.2; nz=80*g;%g为每时刻的反向接入的概率 %前几个的最小/最大 MultiObj.var_min(1:nz) = -15; MultiObj.var_max(1:nz) = 0; MultiObj.var_min(nz+1:80) = 0; MultiObj.var_max(nz+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 = 7; % 最大迭代次数 %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 % MOPSO REP = MOPSO(params,MultiObj);%MOPSO的迭代次数决定了REP的行数,而M函数的值是一个迭代下x的所有值决定了一行中所有列的值 %所以,每次输出80个(一行) %也即谐波计算要返回80个值 disp('REP.pos_fit 的大小:'); size(REP.pos_fit) disp('REP.pos_fit 的内容:'); disp(REP.pos_fit) disp('REP.pos 的大小:'); size(REP.pos) disp('REP.pos 的内容:'); disp(REP.pos) 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); a=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); %figure; %plot(REP.pos_fit(:,1), REP.pos_fit(:,2), 'ro'); %xlabel('f1 (THDu)'); %ylabel('f2 (Voltage Deviation)'); %title('Pareto Front'); %grid on; %iter作用, 对抗函数 function f= m(x) %%这个函数就是给一组粒子,返回一个sumTHDU! %% n=size(x,1);%粒子数。 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 f1 = zeros(n, 1); f2 = zeros(n, 1); % 获取充电桩节点信息 charge_bus = charge1(:,1); num_chargers = length(charge_bus); harmonic_num = size(charge1, 2) - 1; % 谐波次数数量 % 初始化目标函数值 total_charge_current = 0; % 用于累积谐波电流绝对值之和 % 提取所有充电桩的谐波电流值 for charger_idx = 1:num_chargers for harmonic_idx = 1:harmonic_num % 计算在x中的索引位置 % 假设x是行向量:每个充电桩的谐波按顺序排列 pos = (charger_idx - 1) * harmonic_num + harmonic_idx; % 获取谐波电流值(实数,表示幅值) current_val = x(pos); % 累加到总谐波电流(取绝对值) total_charge_current = total_charge_current + abs(current_val); end end %% 基于牛拉法的基波潮流计算 %基波下,充电桩的参数设定,与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(Y*Uacs));%公式不能用 % dQ = Qacs - imag(Uacs.*conj(Y*Uacs)); 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)+1i*Bus_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)))); %Ui*sum(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)+1i*ac_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))*(1i*ac_branch(k,5)*ZB)/2; I_branch(k, 3) = I_branch(k, 3) +Ui(to_node(k))*(1i*ac_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)+1i*Qii(:, 2); Iii(:, 2) = conj(Sii)./conj(Ui); %注意共轭的取值 Ii = Ybus*Ui; %Ii = Yi1*U1+Yi2*U2+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 = Sii*ac_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; %次数 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+1i*h*source_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-1i*imag(load_power(j))/(h*Bus_V(j, 2)^2); %yii = Pi/Vi^2-j*Qi/(h*Vi^2)??? % else % yloads(j, 1) = 0; % end % end %线路谐波模型 branch_h = ac_branch;%已经标幺化 branch_h(: ,3) = ac_branch(: ,3); branch_h(: ,4) = h*ac_branch(: ,4); branch_h(: ,5) = h*ac_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)/100; % 所有节点的 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 % 计算充电桩谐波电流绝对值之和 (f2) f22 = total_charge_current; f2(l,1)=f22; % 添加对抗性关系:当谐波电流很小时,强制增加THD % if f22 <500 % 阈值可根据实际情况调整 % f11 = f11 + 5 * (400 - f22); % 惩罚项 % f1(l,1)=f11; %end f=[f1,f2]; end %for iter=1:7 % 打印调试信息 iter=1; disp(['Iteration ', num2str(iter), ':']); disp(['f1: min = ', num2str(min(f1)), ', max = ', num2str(max(f1))]); disp(['f2: min = ', num2str(min(f2)), ', max = ', num2str(max(f2))]); R = corrcoef([f1, f2]); disp(['Correlation between f1 and f2: ', num2str(R(1,2))]); end %d=[U_h;I_h;THDU];%这里不是目标函数,不需要回到example中的REP中,只要在另一个程序中画图即可,那个程序在寻优之后,是运用寻优得到的确切各个电流值来画图即可 请用matlab回答
时间: 2025-07-23 17:55:48 浏览: 4
<think>我们分析用户提供的MATLAB代码,该代码用于实现多目标粒子群优化(MOPSO)来优化33节点系统的总谐波畸变率(THD)。代码主要分为两部分:主脚本和函数`m`。主脚本设置参数并调用MOPSO算法,函数`m`计算目标函数(THD和充电桩谐波电流总和)。
### 代码分析
#### 主脚本部分
1. **加载系统数据**:使用`loadcase`加载33节点系统('case33bw_60')。
2. **定义负荷曲线**:`Load`数组,但未在后续代码中使用。
3. **多目标函数设置**:
- 选择目标函数为`mubiao`(总THD)。
- 决策变量数量为80(10个充电桩节点,每个节点8次谐波电流)。
- 设置变量的上下界:前20%(即16个)变量范围为[-15,0],其余为[0,15](表示反向接入和正向接入的谐波电流)。
4. **MOPSO参数设置**:
- 种群规模`Np=200`,最大迭代次数`maxgen=7`。
- 惯性权重`W=0.7`,个体和社会认知因子`C1=1.0`,`C2=1.0`。
- 其他参数如存储库大小`Nr=200`,网格数`ngrid=10`等。
5. **运行MOPSO**:调用`MOPSO`函数,得到结果`REP`(存储库,包含Pareto最优解)。
6. **结果分析**:
- 输出`REP.pos_fit`(目标函数值)和`REP.pos`(决策变量)的大小和内容。
- 将80个决策变量按充电桩分组(每组8个谐波电流)。
- 找到目标函数值(假设第二个目标是最小网损?)最小的解并显示。
#### 函数`m`
该函数计算两个目标:总谐波畸变率(THD)和充电桩谐波电流绝对值之和。
1. **输入**:`x`为决策变量矩阵,每行代表一个粒子(即一组谐波电流配置)。
2. **循环处理每个粒子**:对每个粒子(即每个解)计算两个目标值。
3. **基波潮流计算**:
- 构建导纳矩阵`Ybus`。
- 进行牛顿-拉夫逊潮流计算,得到节点电压和相角。
4. **谐波潮流计算**:
- 设置谐波源(包括充电桩)的谐波电流。
- 对每个谐波次数(5,7,...,17次)计算谐波潮流:
- 建立负荷和发电机的谐波阻抗模型。
- 构建谐波导纳矩阵`Ybus_h`。
- 求解谐波电压和电流。
- 计算每个节点的总谐波畸变率(THDU)并求和得到总THD(`sum_THDU`)。
5. **第二个目标**:计算充电桩谐波电流绝对值之和(`total_charge_current`)。
6. **输出**:每个粒子的两个目标值(`f1`为总THD,`f2`为谐波电流绝对值之和)。
### 潜在问题
1. **代码效率**:谐波潮流计算使用迭代求解,且嵌套在粒子群优化中,计算量大。最大迭代次数设为7(`maxgen=7`)可能因计算时间而设小,但可能影响收敛。
2. **目标函数`f2`的意义**:`f2`是充电桩谐波电流绝对值之和,用于最小化谐波注入。但谐波电流有正负(相位),绝对值之和可能不是最佳指标(例如,若相位相反可抵消)。
3. **对抗关系**:注释中尝试添加对抗关系(当`f2`小时增加`f1`),但被注释掉。可能需要重新考虑两个目标的关系。
4. **数据未使用**:定义的`Load`数组未使用,可能遗漏负荷变化的影响。
5. **谐波模型**:
- 充电桩谐波电流直接使用决策变量,但未考虑相位(只取幅值?)。
- 谐波源模型(`source1`)中,谐波电流赋值方式(`I_h0`)可能需验证。
6. **潮流收敛**:基波潮流计算中,功率计算使用节点电压和导纳矩阵的显式计算,但收敛判据使用`dP`和`dQ`,需确保正确。
7. **谐波潮流迭代**:谐波潮流计算中,迭代收敛判据为谐波电压变化,但未设置最大迭代次数,可能不收敛。
8. **变量范围**:决策变量范围(-15至15)可能过大,导致谐波电压过大(THD超过限值)。
### 优化建议
1. **减少计算量**:
- 使用向量化操作替代循环。
- 预计算基波潮流的导纳矩阵及其分解(若系统不变)。
2. **目标函数调整**:
- 考虑用THD的最大值(而非总和)作为目标,更符合标准。
- 第二个目标可改为谐波电流的平方和(反映能量)。
3. **参数调整**:
- 增加MOPSO迭代次数(如50-100次)以获得更好的Pareto前沿。
- 调整MOPSO参数(如惯性权重线性递减)。
4. **模型修正**:
- 充电桩谐波电流应考虑相位(复数形式),决策变量可改为幅值和相位。
- 验证谐波源模型是否正确(特别是`I_h0`的计算)。
5. **代码结构**:
- 将基波潮流计算移出循环(若系统拓扑和负荷不变),只需计算一次。
- 谐波潮流计算中,不同谐波次数独立,可并行计算。
### 结论
当前代码框架正确,但需验证谐波潮流计算正确性,并优化计算效率。建议先在小规模系统测试谐波潮流(如固定谐波电流,验证THD结果)。
§§相关问题§§
1. 您希望如何平衡THD最小化和谐波电流绝对值最小化这两个目标?
2. 充电桩谐波电流的相位信息是否重要?是否考虑扩展为复数优化?
3. 系统负荷曲线已定义但未使用,是否需考虑时变负荷下的动态优化?
4. 当前MOPSO迭代次数较少(7代),是否满足收敛需求?是否需要调整参数?
5. 在谐波潮流计算中,如何处理不收敛的情况?是否有备用方案?</think>### 代码分析总结
您的代码实现了基于MOPSO的33节点系统谐波优化,核心是通过调整10个充电桩节点的8次谐波电流(共80个变量),最小化总谐波畸变率(THD)。以下是关键发现和改进建议:
---
### 主要问题
1. **目标函数冲突**
- 当前同时优化两个目标:`f1=总THD` 和 `f2=谐波电流绝对值之和`
- 但二者存在对抗性:降低谐波电流会减少THD,而您的对抗逻辑(注释部分)试图强制二者负相关,这与物理规律矛盾
- **建议**:
```matlab
% 删除对抗性代码(第202-206行),保留自然关系
f = [f1, f2]; % f1和f2实际是正相关
```
2. **谐波潮流计算缺陷**
- 充电桩谐波电流仅用幅值(`abs(current_val)`),未考虑相位角
- 导致谐波抵消效应未被利用(相反相位可降低THD)
- **修改建议**:
```matlab
% 将谐波电流改为复数形式(决策变量包含相位)
MultiObj.nVar = 10*8*2; % 幅值+相位
```
3. **变量范围不合理**
- 反向电流范围 `[-15, 0]` 和正向 `[0, 15]` 未与系统容量关联
- 可能导致电压越限(如 `UB=12.66kV` 系统电流基准值约 `IB=437A`)
- **修正方案**:
```matlab
% 根据节点容量约束范围(示例)
max_current = mpc.bus(:, 5) / (sqrt(3)*UB); % 节点最大电流
```
4. **计算效率瓶颈**
- 每个粒子需计算33节点基波潮流+5次谐波潮流
- 200粒子*7代=1400次全潮流计算,耗时过长
- **优化建议**:
```matlab
% 预计算基波潮流(系统不变时只需1次)
[base_V, base_angle] = calculate_base_power_flow();
% 谐波潮流中复用基波电压
```
---
### 关键BUG修复
**问题**:谐波迭代不收敛(第134-164行)
**原因**:`Zh_eq` 计算错误(未用标幺值)导致阻抗失配
**解决方案**:
```matlab
% 修正阻抗计算(第148行)
Zh_eq = (source_R1 + 1i*h*source_X1) / ZB; % 标幺化
```
**问题**:节点电压越限(`Bus_V` 未校验)
**风险**:优化结果可能违反安全约束
**添加代码**:
```matlab
% 在目标函数中增加惩罚项(第188行后)
violation = max(0, abs(Bus_V(:,2)) - 1.1) + max(0, 0.9 - abs(Bus_V(:,2)));
f1(l) = f11 + 1000 * sum(violation); % 惩罚系数
```
---
### 优化建议
1. **并行计算加速**
```matlab
parfor l = 1:n % 替换for循环
```
2. **灵敏度分析降维**
通过谐波阻抗矩阵识别关键节点,减少优化变量
3. **Pareto前沿筛选**
增加拥挤度计算避免解聚集:
```matlab
REP = updateRepository(REP, new_particles, params.ngrid);
```
---
### 执行验证
运行后检查 `REP.pos_fit`:
- 若 `f1`(THD)范围在1%~5%说明物理合理
- 若出现 `NaN` 检查谐波导纳矩阵奇异性
阅读全文
相关推荐


















