基于某Matlab的PV模型仿真小实验.doc
wordName: Yang , Vorname: Ying Student ID: 359144PV Assignment 1Exercise 1Code:%Exercise 1close allclear all%define input valuesDOY=172;TZ=2;lambda_g=13.2;fi_g=52.3;LT=6:1:20;%call function sundata to get sun enevation and azimuth Am, Ys, SAz = sundata(lambda_g,fi_g,TZ,DOY,LT);plot(SAz,Ys);xlabel('Azimuth in degrees As');ylabel('Altitude Position in degrees Ys');title('Berlin(13.2,52.3) on 21.06.2013(172/365)');grid on;Function sundata:function Am, Ys, SAz = sundata(lambda_g,fi_g,TZ,DOY,LT)%UNTITLED Summary of this function goes here% Detailed explanation goes hereJ=360*DOY/365;TEQ=0.0066+7.3525*cos(J+85.9)*pi/180)+9.9359*cos(2*J+108.9)*pi/180)+0.3387*cos(3*J+105.2)*pi/180);delta=0.3948-23.2559*cos(J+9.1)*pi/180)-0.3915*cos(2*J+5.4)*pi/180)-0.176*cos(3*J+26)*pi/180); % delta units->degreefor k=1:1:15 TLT(k)=LT(k)-TZ+ (4*lambda_g+TEQ)/60; W(k)=(12-TLT(k)*15; Ys(k)=asin( cos(W(k)*pi/180)*cos(fi_g*pi/180)*cos(delta*pi/180) +sin(fi_g*pi/180)*sin(delta*pi/180) )*180/pi; SunA(k)=acos( (sin(Ys(k)*pi/180)*sin(fi_g*pi/180)-sin(delta*pi/180) )/(cos(Ys(k)*pi/180)*cos(fi_g*pi/180) )*180/pi;if LT(k)<=12 SAz(k)=180-SunA(k);else SAz(k)=180+SunA(k);end Am(k)=1./sin(Ys(k)*pi/180);endendFigure 1:Exercise 2Code:% Exercise 2close all;clear all;%define V,G,T to call functionV=0:0.01:44.8;G=1000;T=25;P=zeros(1,4481);%call function PVmod to get I I =PVmod( V, G, T );subplot(2,1,1);plot(V,I);xlabel('Voltage(V)');ylabel('Current(A)');grid on;for i=1:1:4481 P(i)=I(i)*V(i);endsubplot(2,1,2);plot(V,P);xlabel('Voltage(V)');ylabel('Power(W)');grid on;Function PVmod:function I = PVmod( V,G,T )%UNTITLED3 Summary of this function goes here% Detailed explanation goes here%Definations of constants:k=1.38*10(-23);q=1.60*10(-19);Vg=1.12;Voc=44.8;Isc=5.5;Tstc=25+273;Tk=T+273;Ko=0.00065;Rs=0.002;Rsh=800;a,b=size(V);n=1.6;U=zeros();IL=zeros();Is1=zeros();Is=zeros();C=zeros();I=zeros();Voc_1=(Voc-(0.16*(T-25)/72; %Voc for each cell at any TempIsc_1=Isc*(1+Ko*(T-25); %Isc for cell at any Tempfor i=1:1:b U(i)=V(i)/72;if V(i)=0 IL(i)=Isc_1;else IL(i)=Isc_1*(G/1000)*(1+Ko*(Tk-Tstc);end Is1(i)=IL(i)/(exp(q*Voc_1/(n*k*Tstc)-1);%Is at STC Is(i)=Is1(i)*(Tk/Tstc)(3/n)*exp(-q*Vg*(1/Tk-1/Tstc)/(n*k);%Is for any Temp Tkfor j=1:1:10if j=1; C(j)=4.5;end%Newton's method: difI=-1-(q*Is*Rs/(n*k*Tk)*exp(q*(U(i)+C(j)*Rs)/(n*k*Tk)-Rs/Rsh; fI=IL-C(j)-Is*(exp(q*(U(i)+C(j)*Rs)/(n*k*Tk)-1)-(U(i)+C(j)*Rs)/Rsh; err=-fI/difI; C(j+1)=C(j)+err;end I(i)=C(j);endendFigure 2:Explain:Best match for n, Rs and Rsh is that:Rsh=800Exercise 3Code:% Exercise 3close all;clear all;G=1000;A=1593*790/10000;Vmpp=36.5;Impp=5.1;Voc=44.8;Isc=5.5;etamax=(Vmpp*Impp)/(A*G);FF=(Vmpp*Impp)/(Voc*Isc);Result 3:FF=75.55%Exercise 4Code:% Exercise 4close all;clear all;V=0:0.01:44.8;%to get different values of G and TG=input('Please input the value of Irrad:');T=input('Please input the value of Temperature:');P=zeros(1,4481);%call function PVmod to get I I =PVmod( V, G, T );%get Pfor i=1:1:4481 P(i)=I(i)*V(i);end%find the position of Pmax,then find its corresponding Vm and Imr,c = find(P = max(P(:);Vm=V(r,c);Im=I(r,c);Rm=Vm/Im;Result 4:(1)(2)(3)(4)(5)(6)Irradiance(W/m2)20020060060010001000Temperature(oC)025025025Resistive Load()Exercise 5Code:% Exercise 5close all;clear all;%to get different values of G, T and RG=input('Please input the value of Irrad:');T=input('Please input the value of Temperature:');R=input('Please input the value of R:');V=(0:0.01:44.8);P=zeros(1,4481);R1=zeros(1,4481);delta_R=zeros(1,4481);%call function to get I I =PVmod( V, G, T );%get P and R for each stepfor i=1:1:4481P(i)=V(i)*I(i);R1(i)=V(i)/I(i);delta_R(i)=abs(R-R1(i);end%find the corresponding V and I for the given R(approx)r,c=find(delta_R = min(delta_R(:); V1=V(r,c);I1=I(r,c);P1=V1*I1;Pm=max(P);Pd=(P1/Pm)*100;Result 5:Percentage power difference for G=1000, T=25, R=5: Pd=81.0864%Exercise 6Code:%Exercise 6close all;clear all;delta_V=0.4;Vo=30;G=input('Please input the value of Irrad:');T=input('Please input the value of Temperature:');%call function MPP to get the D and Vnew D,Vnew = MPP( G,T,delta_V,Vo );subplot(2,1,1);plot(D);xlabel('Time steps');ylabel('Duty cycle');grid on;subplot(2,1,2);plot(Vnew);xlabel('Time steps');ylabel('Module operating voltage(V)');grid on;Function MPP:function D,Vnew = MPP( G,T,delta_V ,Vo)%UNTITLED4 Summary of this function goes here% Detailed explanation goes hereV=zeros(1,400);Vnew=zeros(1,400);Vstep=zeros(1,400);P=zeros(1,400);D=zeros(1,400);V(1)=0;D(1)=Vo/(V(1)+Vo);I=PVmod(V(1),G,T);P(1)=I*V(1);Vstep(1)=delta_V;Vnew(1)=V(1)+Vstep(1);for i=1:1:400 V(i+1)=V(i)+Vstep(i); I=PVmod(V(i+1),G,T); P(i+1)=V(i+1)*I;if (P(i+1)-P(i)>0 Vstep(i+1)=Vstep(i);else Vstep(i+1)=-Vstep(i);end D(i+1)=Vo/(Vo+V(i+1); Vnew(i+1)=V(i+1)+Vstep(i+1);endendFigure 6(G=1000, T=25) would be given together with the figure 7Exercise 7Code:% Exercise 7close all;clear all;delta_V=0.1;Vo=30;G=input('Please input the value of Irrad:');T=input('Please input the value of Temperature:'); %call function to get D and Vnew D,Vnew = MPP( G,T,delta_V,Vo ); subplot(2,1,1);plot(D);xlabel('Time steps');ylabel('Duty cycle');grid on;subplot(2,1,2);plot(Vnew);xlabel('Time steps');ylabel('Module operating voltage(V)');grid on;Figure 6(G=1000, T=25): Figure 7(G=1000, T=25):pare with the figures of exercise 6 and 7, get the conclusion:The more the Voltage step is, the speed of tracking is faster but the accuracy of tracking is lower. Exercise 8Code:% Exercise 8close all;clear all;Vo=30;G=input('Please input the value of Irrad:');T=input('Please input the value of Temperature:');%call function to get D and Vnew D,Vnew,P = AHCA( G,T,Vo );subplot(2,1,1);plot(D);xlabel('Time steps');ylabel('Duty cycle');grid on;subplot(2,1,2);plot(Vnew);xlabel('Time steps');ylabel('Module operating voltage(V)');grid on;Figure 8(G=1000, T=25): Figure 7(G=1000, T=25):pare with the figures of exercise 8 and 7, get the conclusion that: the adaptive hill climbing can much faster track the maximum power point and the accuracy is also better.function D, Vnew, P = AHCA( G,T ,Vo)%UNTITLED3 Summary of this function goes here% Detailed explanation goes hereV=zeros(1,400);Vnew=zeros(1,400);Vstep=zeros(1,400);P=zeros(1,400);D=zeros(1,400);dP=zeros(1,400);dV=zeros(1,400);%define the first two stepsV(1)=0;D(1)=Vo/(V(1)+Vo);I=PVmod(V(1),G,T);P(1)=I*V(1);Vstep(1)=0.4;V(2)=V(1)+Vstep(1);I=PVmod(V(2),G,T);P(2)=V(2)*I;Vnew(1)=V(1)+Vstep(1);D(2)=Vo/(V(2)+Vo);%Voltage step is dP/dVfor i=2:1:400 dP(i)=P(i)-P(i-1); dV(i)=V(i)-V(i-1); Vstep(i)=dP(i)/dV(i); Vnew(i)=V(i)+Vstep(i); V(i+1)=V(i)+ Vstep(i); %call function PVmod to get corresponding I and P I=PVmod(V(i+1),G,T); P(i+1)=V(i+1)*I; D(i+1)=Vo/(Vo+V(i+1);endend13 / 13
