使用粒子群PSO算法实现MPPT-M语言仿真

电力电子 同时被 2 个专栏收录
12 篇文章 41 订阅
3 篇文章 2 订阅

在Octave以及Matlab上,仿真了使用粒子群PSO实现MPPT的过程。粒子数为4。太阳能电池为4个串联。

2019年4月24日更新matlab代码。

目录

1.1 先绘制出PV曲线(Octave)

1.2 PSO算法(Octave)

2.1 绘制PV曲线(Matlab)

2.2 PSO.m(Matlab)

3 仿真结果


 


本文主要是代码。

我的软件环境是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循环内设置一个断点,那么每次迭代的结果都会在图形上显示。

 

评论34
请先登录 后发表评论~
©️2021 CSDN 皮肤主题: 编程工作室 设计师:CSDN官方博客 返回首页

打赏作者

qq_27158179

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值