在Octave以及Matlab上，仿真了使用粒子群PSO实现MPPT的过程。粒子数为4。太阳能电池为4个串联。
2019年4月24日更新matlab代码。
目录
本文主要是代码。
我的软件环境是winxp(32bit)，Octave4.4.0。Matlab2010b。
Octave和Matlab的函数有些差异。要在Octave上运行，请使用1.1和1.2的代码。要在Matlab上运行 ，请使用2.1和2.2的代码。
1.1 先绘制出PV曲线（Octave）
光伏电阻串联的方法如下图。每个光伏电池都有各自的反向并联二极管。此处二极管的作用是保护光伏电池避免经过大电流，造成过度发热。
代码中仿真得到了，由4个光伏电池串联成的光伏电池组特性曲线，如下图。
红色线代表理想情况下的曲线。而蓝色线是串联起来的光伏电池组的PV特性曲线。没达到理想曲线的原因是，每个光伏电池都受到不同程度的遮挡。以下是代码。
新建一个文件：mySolarPannel_4inseries.m。把以下代码复制进去：
clear
clc
%-----------------------------------------------
%-----------------------------------------------
%pannel in series
%first pannel
S_1=1000;
Tair_1=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_1 = Tair_1 + 0.028*S_1;
T_delta_1 = T_1 - Tref;
S_delta_1 = S_1/Sref - 1;
Isc_comp_1 = Isc*S_1/Sref*(1+a*T_delta_1);
Uoc_comp_1 = Uoc*(1-c*T_delta_1)*log(e+b*S_delta_1);
Im_comp_1 = Im*S_1/Sref*(1+a*T_delta_1);
Um_comp_1 = Um*(1-c*T_delta_1)*log(e+b*S_delta_1);
C2_1=(Um_comp_1/Uoc_comp_1-1)*(log(1-Im_comp_1/Isc_comp_1))^(-1);
C1_1=(1-Im_comp_1/Isc_comp_1)*exp(-Um_comp_1/(C2_1*Uoc_comp_1));
U_1=0:0.01:Uoc_comp_1;
Iph_1=Isc_comp_1*(1-C1_1*(exp(U_1/(C2_1*Uoc_comp_1))-1));
%-----------------------------------------------
%second pannel
S_2=800;
Tair_2=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_2 = Tair_2 + 0.028*S_2;
T_delta_2 = T_2 - Tref;
S_delta_2 = S_2/Sref - 1;
Isc_comp_2 = Isc*S_2/Sref*(1+a*T_delta_2);
Uoc_comp_2 = Uoc*(1-c*T_delta_2)*log(e+b*S_delta_2);
Im_comp_2 = Im*S_2/Sref*(1+a*T_delta_2);
Um_comp_2 = Um*(1-c*T_delta_2)*log(e+b*S_delta_2);
C2_2=(Um_comp_2/Uoc_comp_2-1)*(log(1-Im_comp_2/Isc_comp_2))^(-1);
C1_2=(1-Im_comp_2/Isc_comp_2)*exp(-Um_comp_2/(C2_2*Uoc_comp_2));
U_2=0:0.01:Uoc_comp_2;
Iph_2=Isc_comp_2*(1-C1_2*(exp(U_2/(C2_2*Uoc_comp_2))-1));
%-----------------------------------------------
%third pannel
S_3=600;
Tair_3=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_3 = Tair_3 + 0.028*S_3;
T_delta_3 = T_3 - Tref;
S_delta_3 = S_3/Sref - 1;
Isc_comp_3 = Isc*S_3/Sref*(1+a*T_delta_3);
Uoc_comp_3 = Uoc*(1-c*T_delta_3)*log(e+b*S_delta_3);
Im_comp_3 = Im*S_3/Sref*(1+a*T_delta_3);
Um_comp_3 = Um*(1-c*T_delta_3)*log(e+b*S_delta_3);
C2_3=(Um_comp_3/Uoc_comp_3-1)*(log(1-Im_comp_3/Isc_comp_3))^(-1);
C1_3=(1-Im_comp_3/Isc_comp_3)*exp(-Um_comp_3/(C2_3*Uoc_comp_3));
U_3=0:0.01:Uoc_comp_3;
Iph_3=Isc_comp_3*(1-C1_3*(exp(U_3/(C2_3*Uoc_comp_3))-1));
%-----------------------------------------------
%forth pannel
S_4=400;
Tair_4=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_4 = Tair_4 + 0.028*S_4;
T_delta_4 = T_4 - Tref;
S_delta_4 = S_4/Sref - 1;
Isc_comp_4 = Isc*S_4/Sref*(1+a*T_delta_4);
Uoc_comp_4 = Uoc*(1-c*T_delta_4)*log(e+b*S_delta_4);
Im_comp_4 = Im*S_4/Sref*(1+a*T_delta_4);
Um_comp_4 = Um*(1-c*T_delta_4)*log(e+b*S_delta_4);
C2_4=(Um_comp_4/Uoc_comp_4-1)*(log(1-Im_comp_4/Isc_comp_4))^(-1);
C1_4=(1-Im_comp_4/Isc_comp_4)*exp(-Um_comp_4/(C2_4*Uoc_comp_4));
U_4=0:0.01:Uoc_comp_4;
Iph_4=Isc_comp_4*(1-C1_4*(exp(U_4/(C2_4*Uoc_comp_4))-1));
plot(U_1,Iph_1)
plot(U_2,Iph_2)
plot(U_3,Iph_3)
plot(U_4,Iph_4)
%-----------------------------------------------
% 4 in series
% U=C2*Uoc*log((Isc-I)/(C1*Isc)+1)
%Iph_1 > Iph_2 > Iph_3
U_s = 0:0.01:Uoc_comp_1*4;
Iph_s=Isc_comp_1*(1-C1_1*(exp(U_s/(C2_1*Uoc_comp_1*4))-1));
plot(U_s,Iph_s,'k')
U_ss = zeros(1,size(U_1)(2)+size(U_2)(2)+size(U_3)(2)+size(U_4)(2));
Iph_ss = U_ss;
%for i = 1 : size(U_1)(2)
% U_ss(i) = U_1(i);
% Iph_ss(i) = Iph_1(i);
%end
for i = 1 : size(U_1)(2)
if Iph_1(i)>=Iph_2(1)
U_ss(i) = U_1(i);
Iph_ss(i) = Iph_1(i);
step1 = i;
else
break;
end
end
for i = 1 : size(U_2)(2)
if Iph_2(i)>Iph_3(1)
U_ss(step1+i) = U_2(i) + U_ss(step1);
Iph_ss(step1+i) = Iph_2(i);
step2 = step1+i;
else
break;
end
end
for i = 1 : size(U_3)(2)
if Iph_3(i)>Iph_4(1)
U_ss(step2+i) = U_3(i) + U_ss(step2);
Iph_ss(step2+i) = Iph_3(i);
step3 = step2+i;
else
break;
end
end
for i = 1 : size(U_4)(2)
U_ss(step3+i) = U_ss(step3) + U_4(i);
Iph_ss(step3+i) = Iph_4(i);
end
%plot(U_1,Iph_1) I-U
%plot(U_2,Iph_2)
%plot(U_ss,Iph_ss,'+')
%plot(U_ss,Iph_ss,'+')
figure(1)
plot(U_ss,Iph_ss,'+')
xlabel('U/V')
ylabel('I/A')
title('U-I')
figure(2)
P_ss = U_ss .* Iph_ss;
plot(U_ss,P_ss)
xlabel('U/V')
ylabel('P/W')
title('U-W')
hold on
P_1 = U_1*4 .* Iph_1;
plot(U_1*4,P_1)
执行完以后，会得到电压电流特性曲线见figure(1)和电压功率特性曲线见figure(2)，还有好多参数在工作区内。
1.2 PSO算法（Octave）
算法的作用是，在PV曲线中均匀布防4只个体。4只个体会共享数据，并能引导群体往最优值靠近。因此并不会陷入局部最优，可以找到整体最优解。
新建一个文件：PSO.m。把以下内容复制进去：（20190414更新了代码。我是电气工程背景的，编的代码没有那些开源代码那么清晰请多多见谅。）
gbest = 0;
# initialize
currentN = 1;
maxN = 100;
#2.
N=4;
Uoc_module = Uoc_comp_4;
Uoc_array = Uoc_comp_4 * N;
v_i = 0.25/(N-1)*Uoc_module;
x_max = 0.95 * Uoc_array;
v_max = Uoc_module;
i=1;
Vo_1 = 0.7* (i-1)*Uoc_module; #i<N
i=2;
Vo_2 = 0.7* (i-1)*Uoc_module; #i<N
i=3;
Vo_3 = 0.7* (i-1)*Uoc_module; #i<N
i=4;
Vo_4 = 0.7* Uoc_array; #i=N;
Uoc_array = Uoc * N;
v_i = 0.25/(N-1)*Uoc_module;
x_max = 0.95 * Uoc_array;
v_max = Uoc_module;
P_1_max = 0;
P_2_max = 0;
P_3_max = 0;
P_4_max = 0;
Power_1_max = 0;
Power_2_max = 0;
Power_3_max = 0;
Power_4_max = 0;
f_g = 0;
f_g_voltage =0;
c1_begin = 2.7;
c1_end = 1.2;
c2_begin = 0.5;
c2_end = 2.2;
v_1 = 0.25/(N-1)*Uoc_module;
v_2 = 0.25/(N-1)*Uoc_module;
v_3 = 0.25/(N-1)*Uoc_module;
v_4 = 0.25/(N-1)*Uoc_module;
#----------------------------------------------------------------------
#Power_1 = GetResult_U_to_P(Vo_1);
if round(Vo_1*100)==0
Power_1 = P_ss(1);
else
Power_1=P_ss(round(Vo_1*100));
end
#Power_2 = GetResult_U_to_P(Vo_2);
if round(Vo_2*100)==0
Power_2 = P_ss(1);
else
Power_2=P_ss(round(Vo_2*100));
end
#Power_3 = GetResult_U_to_P(Vo_3);
if round(Vo_3*100)==0
Power_3 = P_ss(1);
else
Power_3=P_ss(round(Vo_3*100));
end
#Power_4 = GetResult_U_to_P(Vo_4);
if round(Vo_4*100)==0
Power_4 = P_ss(1);
else
Power_4=P_ss(round(Vo_4*100));
end
Power_N = [Power_1,Power_2,Power_3,Power_4]
Vo_N = [Vo_1, Vo_2, Vo_3, Vo_4];
disp(currentN)
disp(Vo_N);
disp(Power_N);
for currentN = 1 : maxN-1
r1 = rand();
r2 = rand();
c1 = c1_begin + (c1_end - c1_begin)*(1-acos(-2*currentN/maxN+1)/pi);
c2 = c2_begin + (c2_end - c2_begin)*(1-acos(-2*currentN/maxN+1)/pi);
#f_avg = sum(f_i)/N; #paticle power average
#f_avg = (P_1+ P_2 + P_3 + P_4)/4;
f_avg = sum(Power_N)/N;
if f_g < max(Power_N)
f_g = max(Power_N);
for j=1:N
if Power_N(j) == max(Power_N)
f_g_voltage = Vo_N(j);
end
end
end
phi = abs(f_g-f_avg); #f_g: the optimized one
f_1 = Power_1;
w_1=abs((f_1-f_avg)/phi);
f_2 = Power_2;
w_2=abs((f_2-f_avg)/phi);
f_3 = Power_3;
w_3=abs((f_3-f_avg)/phi);
f_4 = Power_4;
w_4=abs((f_4-f_avg)/phi);
w_N = [w_1, w_2, w_3, w_4];
if Power_1_max < Power_1
P_1_max = Vo_1;
Power_1_max = Power_1;
end
if Power_2_max < Power_2
P_2_max = Vo_2;
Power_2_max = Power_2;
end
if Power_3_max < Power_3
P_3_max = Vo_3;
Power_3_max = Power_3;
end
if Power_4_max < Power_4
P_4_max = Vo_4;
Power_4_max = Power_4;
end
P_N_max = [P_1_max, P_2_max, P_3_max, P_4_max];
Power_N_max = [Power_1_max, Power_2_max, Power_3_max, Power_4_max];
v_1_next = w_1 * v_1 + c1*r1*(P_1_max - Vo_1) + c2*r2*(f_g_voltage - Vo_1);
v_2_next = w_2 * v_2 + c1*r1*(P_2_max - Vo_2) + c2*r2*(f_g_voltage - Vo_2);
v_3_next = w_3 * v_3 + c1*r1*(P_3_max - Vo_3) + c2*r2*(f_g_voltage - Vo_3);
v_4_next = w_4 * v_4 + c1*r1*(P_4_max - Vo_4) + c2*r2*(f_g_voltage - Vo_4);
if v_1_next > Uoc_module
v_1_next = Uoc_module;
end
if v_2_next > Uoc_module
v_2_next = Uoc_module;
end
if v_3_next > Uoc_module
v_3_next = Uoc_module;
end
if v_4_next > Uoc_module
v_4_next = Uoc_module;
end
v_N_next = [v_1_next, v_2_next, v_3_next, v_4_next];
Vo_1 = Vo_1 + v_1_next;
Vo_2 = Vo_2 + v_2_next;
Vo_3 = Vo_3 + v_3_next;
Vo_4 = Vo_4 + v_4_next;
if Vo_1 < 0
Vo_1 = 0;
end
if Vo_2 < 0
Vo_2 = 0;
end
if Vo_3 < 0
Vo_3 = 0;
end
if Vo_4 < 0
Vo_4 = 0;
end
if Vo_1 > 0.91*Uoc_array
Vo_1 = 0.91*Uoc_array;
end
if Vo_2 > 0.91*Uoc_array
Vo_2 = 0.91*Uoc_array;
end
if Vo_3 > 0.91*Uoc_array
Vo_3 = 0.91*Uoc_array;
end
if Vo_4 > 0.91*Uoc_array
Vo_4 = 0.91*Uoc_array;
end
#Power_1 = GetResult_U_to_P(Vo_1);
if round(Vo_1*100)==0
Power_1 = P_ss(1);
else
Power_1=P_ss(round(Vo_1*100));
end
#Power_2 = GetResult_U_to_P(Vo_2);
if round(Vo_2*100)==0
Power_2 = P_ss(1);
else
Power_2=P_ss(round(Vo_2*100));
end
#Power_3 = GetResult_U_to_P(Vo_3);
if round(Vo_3*100)==0
Power_3 = P_ss(1);
else
Power_3=P_ss(round(Vo_3*100));
end
#Power_4 = GetResult_U_to_P(Vo_4);
if round(Vo_4*100)==0
Power_4 = P_ss(1);
else
Power_4=P_ss(round(Vo_4*100));
end
Power_N = [Power_1,Power_2,Power_3,Power_4]
Vo_N = [Vo_1, Vo_2, Vo_3, Vo_4];
disp(currentN)
disp('P_N_max: individual output voltage of maximum power')
disp(P_N_max);
disp('Power_N_max: individual maximum power output')
disp(Power_N_max)
disp('f_g: global max')
disp(f_g)
disp('Vo_N: output voltage')
disp(Vo_N);
disp('power_N: output power')
disp(Power_N);
#disp('v_N_next')
#disp(v_N_next)
v_1 = v_1_next;
v_2 = v_2_next;
v_3 = v_3_next;
v_4 = v_4_next;
v_N = [v_1, v_2, v_3, v_4];
differ = (abs(Vo_1-sum( Vo_N)/N) + abs(Vo_2-sum( Vo_N)/N) + abs(Vo_2-sum( Vo_N)/N) + abs(Vo_4-sum( Vo_N)/N))
if gbest ==0
gbest = f_g;
else
gbest = [gbest, f_g];
endif
figure(2)
clf
plot(U_ss,P_ss)
xlabel('U/V')
ylabel('P/W')
title('U-W')
hold on
plot(Vo_N, Power_N, 'k+')
if differ< (0.1*Uoc_array)
break;
end
end
#v_i_next = w*v_i_k + c1*r1*(P_max - V_i_output) + c2*r2*(G_max - V_i_output);
#plot(f_g_voltage, f_g, 'k+')
#figure(2)
#hold on
#plot(v_N, Power_N, 'k+')2.1 绘制PV曲线（Matlab）
新建一个文件：mySolarPannel_4inseries.m。将以下代码复制进去。
clear
clc
%-----------------------------------------------
%-----------------------------------------------
%pannel in series
%first pannel
S_1=1000;
Tair_1=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_1 = Tair_1 + 0.028*S_1;
T_delta_1 = T_1 - Tref;
S_delta_1 = S_1/Sref - 1;
Isc_comp_1 = Isc*S_1/Sref*(1+a*T_delta_1);
Uoc_comp_1 = Uoc*(1-c*T_delta_1)*log(exp(1)+b*S_delta_1);
Im_comp_1 = Im*S_1/Sref*(1+a*T_delta_1);
Um_comp_1 = Um*(1-c*T_delta_1)*log(exp(1)+b*S_delta_1);
C2_1=(Um_comp_1/Uoc_comp_1-1)*(log(1-Im_comp_1/Isc_comp_1))^(-1);
C1_1=(1-Im_comp_1/Isc_comp_1)*exp(-Um_comp_1/(C2_1*Uoc_comp_1));
U_1=0:0.01:Uoc_comp_1;
Iph_1=Isc_comp_1*(1-C1_1*(exp(U_1/(C2_1*Uoc_comp_1))-1));
%-----------------------------------------------
%second pannel
S_2=800;
Tair_2=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_2 = Tair_2 + 0.028*S_2;
T_delta_2 = T_2 - Tref;
S_delta_2 = S_2/Sref - 1;
Isc_comp_2 = Isc*S_2/Sref*(1+a*T_delta_2);
Uoc_comp_2 = Uoc*(1-c*T_delta_2)*log(exp(1)+b*S_delta_2);
Im_comp_2 = Im*S_2/Sref*(1+a*T_delta_2);
Um_comp_2 = Um*(1-c*T_delta_2)*log(exp(1)+b*S_delta_2);
C2_2=(Um_comp_2/Uoc_comp_2-1)*(log(1-Im_comp_2/Isc_comp_2))^(-1);
C1_2=(1-Im_comp_2/Isc_comp_2)*exp(-Um_comp_2/(C2_2*Uoc_comp_2));
U_2=0:0.01:Uoc_comp_2;
Iph_2=Isc_comp_2*(1-C1_2*(exp(U_2/(C2_2*Uoc_comp_2))-1));
%-----------------------------------------------
%third pannel
S_3=600;
Tair_3=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_3 = Tair_3 + 0.028*S_3;
T_delta_3 = T_3 - Tref;
S_delta_3 = S_3/Sref - 1;
Isc_comp_3 = Isc*S_3/Sref*(1+a*T_delta_3);
Uoc_comp_3 = Uoc*(1-c*T_delta_3)*log(exp(1)+b*S_delta_3);
Im_comp_3 = Im*S_3/Sref*(1+a*T_delta_3);
Um_comp_3 = Um*(1-c*T_delta_3)*log(exp(1)+b*S_delta_3);
C2_3=(Um_comp_3/Uoc_comp_3-1)*(log(1-Im_comp_3/Isc_comp_3))^(-1);
C1_3=(1-Im_comp_3/Isc_comp_3)*exp(-Um_comp_3/(C2_3*Uoc_comp_3));
U_3=0:0.01:Uoc_comp_3;
Iph_3=Isc_comp_3*(1-C1_3*(exp(U_3/(C2_3*Uoc_comp_3))-1));
%-----------------------------------------------
%forth pannel
S_4=400;
Tair_4=25;
Sref=1000; %1000W/m^2
Tref=25; %25degree celcius
Uoc=44.2;
Um=35.4;
Isc=5.29;
Im=4.95;
a=0.00255;
b=0.55;
c=0.00285;
T_4 = Tair_4 + 0.028*S_4;
T_delta_4 = T_4 - Tref;
S_delta_4 = S_4/Sref - 1;
Isc_comp_4 = Isc*S_4/Sref*(1+a*T_delta_4);
Uoc_comp_4 = Uoc*(1-c*T_delta_4)*log(exp(1)+b*S_delta_4);
Im_comp_4 = Im*S_4/Sref*(1+a*T_delta_4);
Um_comp_4 = Um*(1-c*T_delta_4)*log(exp(1)+b*S_delta_4);
C2_4=(Um_comp_4/Uoc_comp_4-1)*(log(1-Im_comp_4/Isc_comp_4))^(-1);
C1_4=(1-Im_comp_4/Isc_comp_4)*exp(-Um_comp_4/(C2_4*Uoc_comp_4));
U_4=0:0.01:Uoc_comp_4;
Iph_4=Isc_comp_4*(1-C1_4*(exp(U_4/(C2_4*Uoc_comp_4))-1));
%{
plot(U_1,Iph_1)
hold on
plot(U_2,Iph_2)
hold on
plot(U_3,Iph_3)
hold on
plot(U_4,Iph_4)
%}
%-----------------------------------------------
% 4 in series
% U=C2*Uoc*log((Isc-I)/(C1*Isc)+1)
%Iph_1 > Iph_2 > Iph_3
U_s = 0:0.01:Uoc_comp_1*4;
Iph_s=Isc_comp_1*(1-C1_1*(exp(U_s/(C2_1*Uoc_comp_1*4))-1));
plot(U_s,Iph_s,'k')
U_ss = zeros(size(U_1)+size(U_2)+size(U_3)+size(U_4));
Iph_ss = U_ss;
%for i = 1 : size(U_1)(2)
% U_ss(i) = U_1(i);
% Iph_ss(i) = Iph_1(i);
%end
for i = 1 : length(U_1)
if Iph_1(i)>=Iph_2(1)
U_ss(i) = U_1(i);
Iph_ss(i) = Iph_1(i);
step1 = i;
else
break;
end
end
for i = 1 : length(U_2)
if Iph_2(i)>Iph_3(1)
U_ss(step1+i) = U_2(i) + U_ss(step1);
Iph_ss(step1+i) = Iph_2(i);
step2 = step1+i;
else
break;
end
end
for i = 1 : length(U_3)
if Iph_3(i)>Iph_4(1)
U_ss(step2+i) = U_3(i) + U_ss(step2);
Iph_ss(step2+i) = Iph_3(i);
step3 = step2+i;
else
break;
end
end
for i = 1 : length(U_4)
U_ss(step3+i) = U_ss(step3) + U_4(i);
Iph_ss(step3+i) = Iph_4(i);
end
%plot(U_1,Iph_1) I-U
%plot(U_2,Iph_2)
%plot(U_ss,Iph_ss,'+')
%plot(U_ss,Iph_ss,'+')
figure(1)
plot(U_ss,Iph_ss,'+')
xlabel('U/V')
ylabel('I/A')
title('U-I')
figure(2)
P_ss = U_ss .* Iph_ss;
plot(U_ss,P_ss)
xlabel('U/V')
ylabel('P/W')
title('U-W')
hold on
P_1 = U_1*4 .* Iph_1;
plot(U_1*4,P_1)
2.2 PSO.m（Matlab）
新建一个文件：PSO.m，将以下代码复制进去
gbest = 0;
% initialize
currentN = 1;
maxN = 100;
%2.
N=4;
Uoc_module = Uoc_comp_4;
Uoc_array = Uoc_comp_4 * N;
v_i = 0.25/(N-1)*Uoc_module;
x_max = 0.95 * Uoc_array;
v_max = Uoc_module;
i=1;
Vo_1 = 0.7* (i-1)*Uoc_module; %i<N
i=2;
Vo_2 = 0.7* (i-1)*Uoc_module; %i<N
i=3;
Vo_3 = 0.7* (i-1)*Uoc_module; %i<N
i=4;
Vo_4 = 0.7* Uoc_array; %i=N;
Uoc_array = Uoc * N;
v_i = 0.25/(N-1)*Uoc_module;
x_max = 0.95 * Uoc_array;
v_max = Uoc_module;
P_1_max = 0;
P_2_max = 0;
P_3_max = 0;
P_4_max = 0;
Power_1_max = 0;
Power_2_max = 0;
Power_3_max = 0;
Power_4_max = 0;
f_g = 0;
f_g_voltage =0;
c1_begin = 2.7;
c1_end = 1.2;
c2_begin = 0.5;
c2_end = 2.2;
v_1 = 0.25/(N-1)*Uoc_module;
v_2 = 0.25/(N-1)*Uoc_module;
v_3 = 0.25/(N-1)*Uoc_module;
v_4 = 0.25/(N-1)*Uoc_module;
%----------------------------------------------------------------------
%Power_1 = GetResult_U_to_P(Vo_1);
if round(Vo_1*100)==0
Power_1 = P_ss(1);
else
Power_1=P_ss(round(Vo_1*100));
end
%Power_2 = GetResult_U_to_P(Vo_2);
if round(Vo_2*100)==0
Power_2 = P_ss(1);
else
Power_2=P_ss(round(Vo_2*100));
end
%Power_3 = GetResult_U_to_P(Vo_3);
if round(Vo_3*100)==0
Power_3 = P_ss(1);
else
Power_3=P_ss(round(Vo_3*100));
end
%Power_4 = GetResult_U_to_P(Vo_4);
if round(Vo_4*100)==0
Power_4 = P_ss(1);
else
Power_4=P_ss(round(Vo_4*100));
end
Power_N = [Power_1,Power_2,Power_3,Power_4]
Vo_N = [Vo_1, Vo_2, Vo_3, Vo_4];
disp(currentN)
disp(Vo_N);
disp(Power_N);
for currentN = 1 : maxN-1
r1 = rand();
r2 = rand();
c1 = c1_begin + (c1_end - c1_begin)*(1-acos(-2*currentN/maxN+1)/pi);
c2 = c2_begin + (c2_end - c2_begin)*(1-acos(-2*currentN/maxN+1)/pi);
%f_avg = sum(f_i)/N; %paticle power average
%f_avg = (P_1+ P_2 + P_3 + P_4)/4;
f_avg = sum(Power_N)/N;
if f_g < max(Power_N)
f_g = max(Power_N);
for j=1:N
if Power_N(j) == max(Power_N)
f_g_voltage = Vo_N(j);
end
end
end
phi = abs(f_g-f_avg); %f_g: the optimized one
f_1 = Power_1;
w_1=abs((f_1-f_avg)/phi);
f_2 = Power_2;
w_2=abs((f_2-f_avg)/phi);
f_3 = Power_3;
w_3=abs((f_3-f_avg)/phi);
f_4 = Power_4;
w_4=abs((f_4-f_avg)/phi);
w_N = [w_1, w_2, w_3, w_4];
if Power_1_max < Power_1
P_1_max = Vo_1;
Power_1_max = Power_1;
end
if Power_2_max < Power_2
P_2_max = Vo_2;
Power_2_max = Power_2;
end
if Power_3_max < Power_3
P_3_max = Vo_3;
Power_3_max = Power_3;
end
if Power_4_max < Power_4
P_4_max = Vo_4;
Power_4_max = Power_4;
end
P_N_max = [P_1_max, P_2_max, P_3_max, P_4_max];
Power_N_max = [Power_1_max, Power_2_max, Power_3_max, Power_4_max];
v_1_next = w_1 * v_1 + c1*r1*(P_1_max - Vo_1) + c2*r2*(f_g_voltage - Vo_1);
v_2_next = w_2 * v_2 + c1*r1*(P_2_max - Vo_2) + c2*r2*(f_g_voltage - Vo_2);
v_3_next = w_3 * v_3 + c1*r1*(P_3_max - Vo_3) + c2*r2*(f_g_voltage - Vo_3);
v_4_next = w_4 * v_4 + c1*r1*(P_4_max - Vo_4) + c2*r2*(f_g_voltage - Vo_4);
if v_1_next > Uoc_module
v_1_next = Uoc_module;
end
if v_2_next > Uoc_module
v_2_next = Uoc_module;
end
if v_3_next > Uoc_module
v_3_next = Uoc_module;
end
if v_4_next > Uoc_module
v_4_next = Uoc_module;
end
v_N_next = [v_1_next, v_2_next, v_3_next, v_4_next];
Vo_1 = Vo_1 + v_1_next;
Vo_2 = Vo_2 + v_2_next;
Vo_3 = Vo_3 + v_3_next;
Vo_4 = Vo_4 + v_4_next;
if Vo_1 < 0
Vo_1 = 0;
end
if Vo_2 < 0
Vo_2 = 0;
end
if Vo_3 < 0
Vo_3 = 0;
end
if Vo_4 < 0
Vo_4 = 0;
end
if Vo_1 > 0.91*Uoc_array
Vo_1 = 0.91*Uoc_array;
end
if Vo_2 > 0.91*Uoc_array
Vo_2 = 0.91*Uoc_array;
end
if Vo_3 > 0.91*Uoc_array
Vo_3 = 0.91*Uoc_array;
end
if Vo_4 > 0.91*Uoc_array
Vo_4 = 0.91*Uoc_array;
end
%Power_1 = GetResult_U_to_P(Vo_1);
if round(Vo_1*100)==0
Power_1 = P_ss(1);
else
Power_1=P_ss(round(Vo_1*100));
end
%Power_2 = GetResult_U_to_P(Vo_2);
if round(Vo_2*100)==0
Power_2 = P_ss(1);
else
Power_2=P_ss(round(Vo_2*100));
end
%Power_3 = GetResult_U_to_P(Vo_3);
if round(Vo_3*100)==0
Power_3 = P_ss(1);
else
Power_3=P_ss(round(Vo_3*100));
end
%Power_4 = GetResult_U_to_P(Vo_4);
if round(Vo_4*100)==0
Power_4 = P_ss(1);
else
Power_4=P_ss(round(Vo_4*100));
end
Power_N = [Power_1,Power_2,Power_3,Power_4]
Vo_N = [Vo_1, Vo_2, Vo_3, Vo_4];
disp(currentN)
disp('P_N_max: individual output voltage of maximum power')
disp(P_N_max);
disp('Power_N_max: individual maximum power output')
disp(Power_N_max)
disp('f_g: global max')
disp(f_g)
disp('Vo_N: output voltage')
disp(Vo_N);
disp('power_N: output power')
disp(Power_N);
%disp('v_N_next')
%disp(v_N_next)
v_1 = v_1_next;
v_2 = v_2_next;
v_3 = v_3_next;
v_4 = v_4_next;
v_N = [v_1, v_2, v_3, v_4];
differ = (abs(Vo_1-sum( Vo_N)/N) + abs(Vo_2-sum( Vo_N)/N) + abs(Vo_2-sum( Vo_N)/N) + abs(Vo_4-sum( Vo_N)/N))
if gbest ==0
gbest = f_g;
else
gbest = [gbest, f_g];
end
figure(2)
clf
plot(U_ss,P_ss)
xlabel('U/V')
ylabel('P/W')
title('U-W')
hold on
plot(Vo_N, Power_N, 'k+')
if differ< (0.1*Uoc_array)
break;
end
end
%v_i_next = w*v_i_k + c1*r1*(P_max - V_i_output) + c2*r2*(G_max - V_i_output);
%plot(f_g_voltage, f_g, 'k+')
%figure(2)
%hold on
%plot(v_N, Power_N, 'k+')
需要注意的是，在Matlab中运行时候，需要把以上两份文件放在同一个文件夹内。先运行mySolarPannel_4inseries.m，然后运行pso.m。
3 仿真结果
初始化将4个个体均匀放置在曲线上。
第7次迭代计算：
第13次迭代计算：
第14次迭代计算：
我每次只允许个体改变2V电压，因此速度有点慢。
第33次迭代计算，个体的差异低于设定值1.768V，因此结束计算：
说明：
将上面提到的这二个文件放在同一个文件夹内。首先执行mySolarPannel_4.m，然后执行PSO.m。
另外在for循环内设置一个断点，那么每次迭代的结果都会在图形上显示。
204
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
