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 celciusUoc=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 celciusUoc=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 celciusUoc=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 celciusUoc=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_3U_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);
%endfor 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;elsebreak;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;elsebreak;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;elsebreak;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)==0Power_1 = P_ss(1);
elsePower_1=P_ss(round(Vo_1*100));    
end#Power_2 = GetResult_U_to_P(Vo_2);
if  round(Vo_2*100)==0Power_2 = P_ss(1);
elsePower_2=P_ss(round(Vo_2*100));    
end
#Power_3 = GetResult_U_to_P(Vo_3);
if  round(Vo_3*100)==0Power_3 = P_ss(1);
elsePower_3=P_ss(round(Vo_3*100));    
end
#Power_4 = GetResult_U_to_P(Vo_4);
if  round(Vo_4*100)==0Power_4 = P_ss(1);
elsePower_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-1r1 = 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:Nif Power_N(j) == max(Power_N)f_g_voltage = Vo_N(j);endendendphi = abs(f_g-f_avg);   #f_g: the optimized onef_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_1P_1_max = Vo_1;Power_1_max = Power_1;endif Power_2_max < Power_2P_2_max = Vo_2;Power_2_max = Power_2;endif Power_3_max < Power_3P_3_max = Vo_3;Power_3_max = Power_3;endif Power_4_max < Power_4P_4_max = Vo_4;Power_4_max = Power_4;endP_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_modulev_1_next = Uoc_module;endif v_2_next > Uoc_modulev_2_next = Uoc_module;endif v_3_next > Uoc_modulev_3_next = Uoc_module;endif v_4_next > Uoc_modulev_4_next = Uoc_module;endv_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 < 0Vo_1 = 0;endif Vo_2 < 0Vo_2 = 0;endif Vo_3 < 0Vo_3 = 0;endif Vo_4 < 0Vo_4 = 0;endif Vo_1 > 0.91*Uoc_arrayVo_1 = 0.91*Uoc_array;endif Vo_2 > 0.91*Uoc_arrayVo_2 = 0.91*Uoc_array;endif Vo_3 > 0.91*Uoc_arrayVo_3 = 0.91*Uoc_array;endif Vo_4 > 0.91*Uoc_arrayVo_4 = 0.91*Uoc_array;end  #Power_1 = GetResult_U_to_P(Vo_1);if  round(Vo_1*100)==0Power_1 = P_ss(1);elsePower_1=P_ss(round(Vo_1*100));    end#Power_2 = GetResult_U_to_P(Vo_2);if  round(Vo_2*100)==0Power_2 = P_ss(1);elsePower_2=P_ss(round(Vo_2*100));    end#Power_3 = GetResult_U_to_P(Vo_3);if  round(Vo_3*100)==0Power_3 = P_ss(1);elsePower_3=P_ss(round(Vo_3*100));    end#Power_4 = GetResult_U_to_P(Vo_4);if  round(Vo_4*100)==0Power_4 = P_ss(1);elsePower_4=P_ss(round(Vo_4*100));    endPower_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 ==0gbest = f_g;elsegbest = [gbest, f_g];endiffigure(2)clfplot(U_ss,P_ss)xlabel('U/V')ylabel('P/W')title('U-W')hold onplot(Vo_N, Power_N, 'k+')if differ< (0.1*Uoc_array)break;endend
#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 celciusUoc=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 celciusUoc=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 celciusUoc=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 celciusUoc=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_3U_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);
%endfor 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;elsebreak;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;elsebreak;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;elsebreak;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)==0Power_1 = P_ss(1);
elsePower_1=P_ss(round(Vo_1*100));    
end%Power_2 = GetResult_U_to_P(Vo_2);
if  round(Vo_2*100)==0Power_2 = P_ss(1);
elsePower_2=P_ss(round(Vo_2*100));    
end
%Power_3 = GetResult_U_to_P(Vo_3);
if  round(Vo_3*100)==0Power_3 = P_ss(1);
elsePower_3=P_ss(round(Vo_3*100));    
end
%Power_4 = GetResult_U_to_P(Vo_4);
if  round(Vo_4*100)==0Power_4 = P_ss(1);
elsePower_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-1r1 = 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:Nif Power_N(j) == max(Power_N)f_g_voltage = Vo_N(j);endendendphi = abs(f_g-f_avg);   %f_g: the optimized onef_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_1P_1_max = Vo_1;Power_1_max = Power_1;endif Power_2_max < Power_2P_2_max = Vo_2;Power_2_max = Power_2;endif Power_3_max < Power_3P_3_max = Vo_3;Power_3_max = Power_3;endif Power_4_max < Power_4P_4_max = Vo_4;Power_4_max = Power_4;endP_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_modulev_1_next = Uoc_module;endif v_2_next > Uoc_modulev_2_next = Uoc_module;endif v_3_next > Uoc_modulev_3_next = Uoc_module;endif v_4_next > Uoc_modulev_4_next = Uoc_module;endv_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 < 0Vo_1 = 0;endif Vo_2 < 0Vo_2 = 0;endif Vo_3 < 0Vo_3 = 0;endif Vo_4 < 0Vo_4 = 0;endif Vo_1 > 0.91*Uoc_arrayVo_1 = 0.91*Uoc_array;endif Vo_2 > 0.91*Uoc_arrayVo_2 = 0.91*Uoc_array;endif Vo_3 > 0.91*Uoc_arrayVo_3 = 0.91*Uoc_array;endif Vo_4 > 0.91*Uoc_arrayVo_4 = 0.91*Uoc_array;end  %Power_1 = GetResult_U_to_P(Vo_1);if  round(Vo_1*100)==0Power_1 = P_ss(1);elsePower_1=P_ss(round(Vo_1*100));    end%Power_2 = GetResult_U_to_P(Vo_2);if  round(Vo_2*100)==0Power_2 = P_ss(1);elsePower_2=P_ss(round(Vo_2*100));    end%Power_3 = GetResult_U_to_P(Vo_3);if  round(Vo_3*100)==0Power_3 = P_ss(1);elsePower_3=P_ss(round(Vo_3*100));    end%Power_4 = GetResult_U_to_P(Vo_4);if  round(Vo_4*100)==0Power_4 = P_ss(1);elsePower_4=P_ss(round(Vo_4*100));    endPower_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 ==0gbest = f_g;elsegbest = [gbest, f_g];endfigure(2)clfplot(U_ss,P_ss)xlabel('U/V')ylabel('P/W')title('U-W')hold onplot(Vo_N, Power_N, 'k+')  if differ< (0.1*Uoc_array)break;endend
%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。