www.pudn.com > STATCOM_BESS.rar > Initial_values_with_statcom.m, change:2012-11-05,size:9561b


% Matlab script to find the initial values of the states and other imp 
% variables for the microgrid model consisting of wind syste, Fuel cell and  
  %microalternator connected to a power system grid and a load  
% Input: Peo-->Microalternator real power output 
%        Pg--> Wind generator system real power output 
%        Load--> load on the grid in per unit 
% Output: Statevar--> An array consisting of initial values of all states 
% Values--> values of the imp system variables which are to be used in the 
% further analysis 
% var--> values of the variables w.rt. operating conditions 
% Author: Mohammad Ashraf Ali @28th_oct_2012 
% Here the analysis is carried out by assuming the terminal voltage of the 
% microalternator and common bus voltage Vs and infinite bus voltage vb are 
% then computed. 
% The number of states are 33 
function [Statevar]=Initial_values_with_statcom(Peo,Pg,Ppv,Pfc,Load) 
 
 global values_MA  
 global values_PV 
 global values_Fc 
 global values_Wind 
QLoad=0.15; 
% Parameters of MicroAlternator 
Sload=Load+1i*QLoad; 
Qeo=0.15; 
w0=377; 
Vto=1.032; 
xdP=0.3;  
xq=0.47; 
 xd=1.3; 
Pm=Peo; 
rt=0.1;xt=0.2; 
  
%Parameters of the wind generation system 
 Rw=0.1;Lw=0.2; % coupling line parameters 
 Rwf=0;Lwf=0.2;Cwf=0.2;Rwdr=0.1; % filter parameters 
 Ra=0.01; xw=Lw; mw=1.02; 
% Parameters of the infinite bus system 
  rb=0.15;  xb=0.26522;  
   
% Parameters of the PV system   
 Rpf=0;Lpf=0.2;Cpf=0.2; Rpdr=0; % filter at output of inverter 
 Rp=0.1;  Lp=0.2; % resistance and inductance of coupling line 
 xp=Lp; % PV reactance coupling since w=1.0 p.u 
 mp=1.02; % modulation index of PV side inverter 
 Ldc=0.05; Cdc=1.1; % boost converter constants 
  %initial values calculation 
  % microalternator model 
pf=cos(atan(Qeo/Peo)); 
its=Peo/(Vto*pf); 
itso=its*(pf-sin(acos(pf))*1i); 
angitso=angle(itso); 
Vs=Vto-(rt+1i*xt)*itso;% voltage at the mid-point 
Vsmag=abs(Vs); 
delVs=angle(Vs); 
eq=Vto+1i*itso*xq; % Location of quadrature Axis 
deltaq=angle(eq); 
theta=pi/2-(deltaq-delVs); 
 
vsd=Vsmag*cos(theta); 
vsq=Vsmag*sin(theta); 
ep=Vto+1i*itso*xdP; 
deltaep=angle(ep); 
eqpo=abs(ep)*cos(deltaq-deltaep); 
vqo=Vto*cos(deltaq); 
vdo=Vto*sin(deltaq)	; 
itd=its*sin(deltaq-angitso); 
itq=its*cos(deltaq-angitso); 
efdo=eqpo+(xd-xdP)*itd; 
  
 Y=conj(Sload/(abs(Vs))^2); 
 g=real(Y); b=-imag(Y); 
iload=Vs*(Y); 
 
%statcom 
rs=0.02; Ls=0.2; ms=1.02; 
chiS=pi/2-(deltaq-delVs);  % angle of statcom bus 
VdcS0=Vsmag/ms;  
isd=0; 
isq=0; 
 
 %wind system 
 Ks=0.2; Xqw=0.7; Xdw=1;W = 1; Wt = 1; 
 Vg = 1.032; Qg = 0.15;    
pfangW=atan2(Qg,Pg);  
pfW =cos(pfangW); 
Sg=Pg+1i*Qg; 
igm=Pg/(Vg*pfW);Ig=igm*(pfW-1i*sin(pfangW)); 
eqW = Vg + Ig*(Ra+1i*Xqw);  
deltaqW = angle(eqW); 
igd = abs(Ig)*sin(deltaqW+pfangW); 
igq = abs(Ig)*cos(deltaqW+pfangW); 
vgd = Vg*sin(deltaqW); 
vgq = Vg*cos(deltaqW); 
Efdw=vgq+Ra*igq+Xdw*igd; 
Tew=Efdw*igq+(Xqw-Xdw)*igd*igq; 
Pew=W*Tew;Pmw = Pew; 
Thetas = Pmw/Ks; % in elect. radians 
Pgcheck =Ra*(igm)^2+(vgd*igd+vgq*igq);  % p.u Turbine Power 
  
%Converter Circuits 
Vwmag=1.0204;Vdcw=Vwmag/mw; 
mC=vgd/(Vdcw*sin(deltaqW)); % converter modulation index 
vgq1=mC*Vdcw*cos(deltaqW);%check 
chiC=0; 
   
%Inverter  
Pw=Pg; 
Z1w=Rwf+1i*Lwf; Yw=1/(Rwdr+1/(1i*Cwf));Z2w=Rw+1i*Lw; 
Zw=Z1w*(1+Yw*Z2w)+Z2w;Zwm=abs(Zw); Zangw=angle(Zw); 
Aw=1+Yw*Z2w; Awm=abs(Aw); Awang=angle(Aw); 
chiw=acos(Zwm/(Vwmag*Vsmag)*(Vwmag^2*Awm/Zwm*cos(Zangw-Awang)-Pw))-Zangw; 
phiWo=chiw+theta; 
angvw=pi/2-deltaq-phiWo; 
Vw=Vwmag*(cos(angvw)-1i*sin(angvw)); 
vwd=Vwmag*cos(pi/2-deltaq-angvw); 
vwq=Vwmag*sin(pi/2-deltaq-angvw); 
iwf=(Vw*Aw-Vs)/Zw;iwfm=abs(iwf);angiwf=angle(iwf); 
iwfd=iwfm*cos(pi/2-deltaq+angiwf); 
iwfq=iwfm*sin(pi/2-deltaq+angiwf); 
Vow=(Vs+iwf*Z2w)/Aw;Vowm=abs(Vow);angVow=angle(Vow); 
Vowd=Vowm*cos(pi/2-deltaq+angVow); 
Vowq=Vowm*sin(pi/2-deltaq+angVow); 
Iw=iwf-Vow*Yw;Iwm=abs(Iw);angiw=angle(Iw); 
iwd=Iwm*cos(pi/2-deltaq+angiw); 
iwq=Iwm*sin(pi/2-deltaq+angiw); 
Vcwd=Vowd-(iwfd-iwd)*Rwdr; 
Vcwq=Vowq-(iwfq-iwq)*Rwdr; 
 
 Pwfcheck=vwd*iwfd+vwq*iwfq; 
 Pwcheck=Vowd*iwd+Vowq*iwq; 
  
 % Fuell cell generation system 
  
 %Fuel cell system parameters 
 Rff=0;Lff=0.2;Cff=0.2; Rfdr=0; % filter at output of inverter 
 Rf=0.1;  Lf=0.2; % resistance and inductance of coupling line 
 xf=Lf; % PV reactance coupling since w=1.0 p.u 
 mf=1.02; % modulation index of PV side inverter 
 Lfc=0.05; Cfc=1.1; % Boost converter constants 
% PEM FC fuel cell modelling  
%test code to find maximu power point 
 Iaf = 0; 
 Pa_maxf = 0; 
% Start process P& O algorithm to find mppt  
while Iaf< 133.3 % Nominal operating point 
Vaf = fuelcell_volt_calc(Iaf); 
Pa_newf = Iaf * Vaf; 
if Pa_newf > Pa_maxf 
Pa_maxf = Pa_newf; 
Impf = Iaf; 
Vmpf = Vaf; 
end 
Iaf = Iaf +0.0005; 
end 
Vbasef=50; % base value for Fuel cell system defined 
Vcfm=1.0204; % have to check from where this has come...???? 
 
 Ibasef=Pa_maxf/(Pfc*Vbasef); % base current in Fuel cell system 
 Ifc=Impf/Ibasef; 
  
Z1f=Rff+1i*Lff; %filter inductor along with resistor Rpf 
Yf=1/(Rfdr+1/(1i*Cff)); % admittance of the parallel branch 
Z2f=Rf+1i*Lf; % impedance of the series coupling branch 
Zf=Z1f*(1+Yf*Z2f)+Z2f; %STAR TO DELTA CONVERSION 
Zfm=abs(Zf); Zangf=angle(Zf); 
Af=1+Yf*Z2f; Afm=abs(Af); Afang=angle(Af); 
 
chif=acos(Zfm/(Vcfm*Vsmag)*(Vcfm^2*Afm/Zfm*cos(Zangf-Afang)-Pfc))-Zangf; % the inverter phase angle psip 
phiof=chif+theta; % initial phase angle psip+theta 
 
angvcf=pi/2-deltaq-phiof; 
Vcf=Vcfm*(cos(angvcf)-1i*sin(angvcf)); % voltage across filter capacitor 
vfd=Vcfm*cos(phiof); % inverter output voltage along d-axis 
vfq=Vcfm*sin(phiof); % inverter output voltage along q-axis 
 
 
 iff=(Vcf*Af-Vs)/Zf;iffm=abs(iff);angiff=angle(iff); % filter current or inverter outpyt current 
 iffd=iffm*cos(pi/2-deltaq+angiff); %along d-axis 
 iffq=iffm*sin(pi/2-deltaq+angiff); %along q-axis 
  
Vof=(Vs+iff*Z2f)/Af;Vofm=abs(Vof);angVof=angle(Vof); %voltage at common node 
Vofd=Vofm*cos(pi/2-deltaq+angVof); 
Vofq=Vofm*sin(pi/2-deltaq+angVof); 
 
If=iff-Vof*Yf;Ifm=abs(If);angif=angle(If); % current to PCC 
ifd=Ifm*cos(pi/2-deltaq+angif); 
ifq=Ifm*sin(pi/2-deltaq+angif); 
Vcfd=Vofd-(iffd-ifd)*Rfdr; % across filter capacitor d-axis 
Vcfq=Vofq-(iffq-ifq)*Rfdr; % across filter capacitor q-axis 
Vdcfo=Vcfm/mf;           % dc-link capacitor or inverter input voltage 
 
% to check the analysis for correctness... 
 Pfcfcheck=vfd*iffd+vfq*iffq; 
 Pfccheck=Vofd*ifd+Vofq*ifq; 
  
 Qfcfcheck=vfq*iffd-vfd*iffq; 
 Qfccheck=Vofq*ifd-Vofd*ifq; 
  
  
 %PV system modelling 
% Pv modelling  
G=1.0; TaC=25; 
Ia = 0; 
Pa_max = 0; 
% Start process P& O algorithm to find mppt  
while Ia < 4.75*G 
Va = bp_sx150svoltcalc(Ia,G,TaC); 
Pa_new = Ia * Va; 
if Pa_new > Pa_max 
Pa_max = Pa_new; 
Imp = Ia; 
Vmp = Va; 
end 
Ia = Ia + .00005; 
end 
Vbase=40; % base value for PV system defined 
Vcpm=1.0204; % have to check from where this has come...???? 
 
 Ibasep=Pa_max/(Ppv*Vbase); % base current in pv system 
 Ipv=Imp/Ibasep; 
  
Z1p=Rpf+1i*Lpf; %filter inductor along with resistor Rpf 
Yp=1/(Rpdr+1/(1i*Cpf)); % admittance of the parallel branch 
Z2p=Rp+1i*Lp; % impedance of the series coupling branch 
Zp=Z1p*(1+Yp*Z2p)+Z2p; Zpm=abs(Zp); Zangp=angle(Zp); 
Ap=1+Yp*Z2p; Apm=abs(Ap); Apang=angle(Ap); 
 
chiP=acos(Zpm/(Vcpm*Vsmag)*(Vcpm^2*Apm/Zpm*cos(Zangp-Apang)-Ppv))-Zangp; % the inverter phase angle psip 
phip=chiP+theta; % initial phase angle psip+theta 
 
angvcp=pi/2-deltaq-phip; 
Vcp=Vcpm*(cos(angvcp)-1i*sin(angvcp)); % voltage across filter capacitor 
vpd=Vcpm*cos(phip); % inverter output voltage along d-axis 
vpq=Vcpm*sin(phip); % inverter output voltage along q-axis 
 
 
 ipf=(Vcp*Ap-Vs)/Zp;ipfm=abs(ipf);angipf=angle(ipf); % filter current or inverter outpyt current 
 ipfd=ipfm*cos(pi/2-deltaq+angipf); %along d-axis 
 ipfq=ipfm*sin(pi/2-deltaq+angipf); %along q-axis 
  
Vop=(Vs+ipf*Z2p)/Ap;Vopm=abs(Vop);angVop=angle(Vop); %voltage at common node 
Vopd=Vopm*cos(pi/2-deltaq+angVop); 
Vopq=Vopm*sin(pi/2-deltaq+angVop); 
 
Ip=ipf-Vop*Yp;Ipm=abs(Ip);angip=angle(Ip); % current to PCC 
ipd=Ipm*cos(pi/2-deltaq+angip); 
ipq=Ipm*sin(pi/2-deltaq+angip); 
Vcpd=Vopd-(ipfd-ipd)*Rpdr; % across filter capacitor d-axis 
Vcpq=Vopq-(ipfq-ipq)*Rpdr; % across filter capacitor q-axis 
Vdcpo=Vcpm/mp;           % dc-link capacitor or inverter input voltage 
 
% to check the analysis for correctness... 
 Ppvfcheck=vpd*ipfd+vpq*ipfq; 
 Ppvcheck=Vopd*ipd+Vopq*ipq; 
  
 Qpvfcheck=vpq*ipfd-vpd*ipfq; 
 Qpvcheck=Vopq*ipd-Vopd*ipq; 
 
 
 % Infinite bus voltage and current 
  
 ibus=itso+Iw+If+Ip-iload; 
 Vbus=Vto-((rt+1i*xt)*itso+(rb+1i*xb)*ibus); 
 Vb=abs(Vbus); 
 delVb=angle(Vbus); 
 delta=deltaq-delVb;    % delta1 is the rotor angle (delta) 
  
 fprintf('The power transfered to the main grid is-->\n'); 
 Sbus=Vbus*conj(ibus) 
 Pbus=real(Sbus) 
 
 Statevar=[delta,1,eqpo,efdo,Ipv,Vdcpo,ipfd,ipfq,ipd,ipq,Vcpd,Vcpq... 
     ,Ifc,Vdcfo,iffd,iffq,ifd,ifq,Vcfd,Vcfq... 
     ,igd,igq,deltaqW,1,Thetas,1,Vdcw,iwfd,iwfq,iwd,iwq,Vcwd,Vcwq,isd,isq,VdcS0]'; 
  
  
 Statevar_MC=[delta,1,eqpo,efdo]'; 
  
 Statevar_Wind=[igd,igq,deltaqW,1,Thetas,1,Vdcw,iwfd,iwfq,iwd,iwfq,Vcwd,Vcwq]'; 
  
 Statevar_fuel_cell=[iffd,iffq,ifd,ifq,Vcfd,Vcfq,Ifc,Vdcfo]'; 
  
 Statevar_PV=[ipfd,ipfq,ipd,ipq,Vcpd,Vcpq,Ipv,Vdcpo]'; 
  
values_MA=[efdo rt xt rb xb Vb g b xd xq xdP itd itq vdo vqo Pm  rs Ls ms chiS]'; 
 
values_PV=[Ldc Cdc Rpf Lpf Cpf Rpdr Rp Lp Ibasep phip mp]'; 
 
values_Fc=[Lfc Cfc Rff Lff Cff Rfdr Rf Lf Ibasef phiof mf]'; 
 
values_Wind=[Ks Efdw  phiWo  mw mC Rw Ra Lw  Rwf Lwf Cwf Rwdr vgd vgq Pmw ]'; 
return