使用粒子群PSO演算法實現MPPT-M語言模擬
阿新 • • 發佈：2019-01-26
在Octave上，模擬了使用粒子群PSO實現MPPT的過程。
目錄
本文主要是程式碼。可以在Octave上執行。我的軟體環境是winxp(32bit)，Octave4.4.0。
1、先繪製出PV曲線
光伏電阻串聯的方法如下圖。每個光伏電池都有各自的反向並聯二極體。此處二極體的作用是保護光伏電池避免經過大電流，造成過度發熱。
程式碼中模擬得到了，由4個光伏電池串聯成的光伏電池組特性曲線，如下圖。
紅色線代表理想情況下的曲線。而藍色線是串聯起來的光伏電池組的PV特性曲線。沒達到理想曲線的原因是，每個光伏電池都受到不同程度的遮擋。以下是程式碼。
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)，還有好多引數在工作區內。
function result=GetResult_P_from_PUcurv(Voutput,P_ss)
if round(Voutput*100)==0
result = P_ss(1);
else
result=P_ss(round(Voutput*100));
end
endfunction把這函式存成GetResult_P_from_PUcurv.m，放在工作目錄內。
這函式功能是根據工作電壓，查表得到當前的輸出功率。在實際專案中，這個功率是經過電壓電流取樣計算得到。
2、PSO演算法
演算法的作用是，在PV曲線中均勻佈防4只個體。4只個體會共享資料，並能引導群體往最優值靠近。因此並不會陷入區域性最優，可以找到整體最優解。
% initialize currentN = 1; maxN = 100; %2. N=4; Uoc_module = Uoc; Uoc_array = Uoc * 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*Uoc_module; %i<N i=2; Vo_2 = 0.7* i*Uoc_module; %i<N i=3; Vo_3 = 0.7* i*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 = 3; 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_P_from_PUcurv(Vo_1,P_ss); Power_2 = GetResult_P_from_PUcurv(Vo_2,P_ss); Power_3 = GetResult_P_from_PUcurv(Vo_3,P_ss); Power_4 = GetResult_P_from_PUcurv(Vo_4,P_ss); 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); endif endfor endif 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 > v_max v_1_next = v_max; end if v_1_next < -v_max v_1_next = -v_max; end if v_2_next > v_max v_2_next = v_max; end if v_2_next < -v_max v_2_next = -v_max; end if v_3_next > v_max v_3_next = v_max; end if v_3_next < -v_max v_3_next = -v_max; end if v_4_next > v_max v_4_next = v_max; end if v_4_next < -v_max v_4_next = -v_max; 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.95*Uoc_array Vo_1 = 0.95*Uoc_array; end if Vo_2 > 0.95*Uoc_array Vo_2 = 0.95*Uoc_array; end if Vo_3 > 0.95*Uoc_array Vo_3 = 0.95*Uoc_array; end if Vo_4 > 0.95*Uoc_array Vo_4 = 0.95*Uoc_array; end Power_1 = GetResult_P_from_PUcurv(Vo_1,P_ss); Power_2 = GetResult_P_from_PUcurv(Vo_2,P_ss); Power_3 = GetResult_P_from_PUcurv(Vo_3,P_ss); Power_4 = GetResult_P_from_PUcurv(Vo_4,P_ss); 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('v_N_next') disp(v_N_next) disp('Vo_N: output voltage') disp(Vo_N); disp('power_N: output power') disp(Power_N); figure(3) clf plot(Vo_N, Power_N,'d') axis([0 200 0 400]) xlabel('U/V') ylabel('P/W') title('U-W') %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 differ< (0.01*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+')
3、模擬結果
初始化將4個個體均勻放置在曲線上。
第7次迭代計算：
第13次迭代計算：
第14次迭代計算：
我每次只允許個體改變2V電壓，因此速度有點慢。
第33次迭代計算，個體的差異低於設定值1.768V，因此結束計算：
