基于Matlab的PV模型仿真小实验

Name: Yang , Vorname: Ying Student ID: 359144 PV Assignment 1

Exercise 1 Code:

%Exercise 1 close all clear all

%define input values DOY=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 here

J=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->degree

for 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); end end

Figure 1:

. . .

. .

Exercise 2 Code:

% Exercise 2 close all; clear all;

%define V,G,T to call function V=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); end

subplot(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 Temp Isc_1=Isc*(1+Ko*(T-25)); %Isc for cell at any Temp

for 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 Tk

for j=1:1:10

if 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);

. . .

. .

end end

Figure 2:

Explain:

Best match for n, Rs and Rsh is that: n=1.6 Rs=0.002 Rsh=800

Exercise 3 Code:

% Exercise 3

. . .

. .

close 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: ŋmax=0.0015 FF=75.55%

Exercise 4 Code:

% Exercise 4 close all;

. . .

. .

clear all;

V=0:0.01:44.8;

%to get different values of G and T

G=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 P

for i=1:1:4481

P(i)=I(i)*V(i); end

%find the position of Pmax,then find its corresponding Vm and Im [r,c] = find(P == max(P(:))); Vm=V(r,c); Im=I(r,c); Rm=Vm/Im;

Result 4: Irradiance(W/m) o2(1) 200 0 43.9205 (2) 200 25 36.2839 (3) 600 0 14.5738 (4) 600 25 12.0328 (5) 1000 0 8.6985 (6) 1000 25 7.1796 Temperature(C) Resistive Load(Ω)

Exercise 5 Code:

. . .

. .

% Exercise 5 close all; clear all;

%to get different values of G, T and R

G=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 step for i=1:1:4481

P(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.08%

Exercise 6 Code:

. . .

. .

%Exercise 6 close all; clear all;

%define voltage step as 0.4V 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 here

V=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); end end

Figure 6(G=1000, T=25) would be given together with the figure 7

Exercise 7 Code:

% Exercise 7 close all; clear all;

%define the new voltage step as 0.1V 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):

. . .

. .

Compare 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 8 Code:

. . .

. .

% Exercise 8 close 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):

Compare 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 here

V=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 steps V(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/dV for 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)); end end

. . .