www.pudn.com > STATCOM_BESS.rar > plot_of_microgrid_voltage.m, change:2012-11-08,size:11908b


 % Matlab script to plot the microgrid voltage and total real power 
 % Date:08th-Nov-12 
 % Matlab script to plot the final composite model uncontrolled, with chis 
% control and chip control 
tic; % start of the time counter; 
clc; clear all; close all;  
 
% Uncontrolled  
[Statevar1]=Initial_values_with_statcom;% calculates the initial conditions 
timerange=[0 6]; 
[t1,SPC1]=ode45('non_linear_36states_unstable_with_battery',timerange,Statevar1); 
 
% Chis controlled 
[Statevar2] = Initial_values_with_statcom_Qload_battery; % calculates initial values 
  
[t2,SPC2]=ode45('non_linear_36states_pi_battery_controlled',timerange,Statevar2); 
 
% Chip controlled 
[Statevar3]=Initial_values_with_statcom_Qload_chiP_controlled; 
 
[t3,SPC3]=ode45('non_linear_36states_chip_controlled_pi_controller',timerange,Statevar3); 
%%%%%%%%%PLOTS%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
toc % end of matlab counter 
 
figure(1);  
plot(t1,SPC1(:,1)*180/pi,'--r','Linewidth',2); 
grid on; hold on; 
plot(t2,SPC2(:,1)*180/pi,'b','Linewidth',2); 
hold on; 
plot(t3,SPC3(:,1)*180/pi,'k','Linewidth',2); 
xlabel('time(seconds)');ylabel('Delta(r)');title('Plot for Delta vs Time'); 
legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
  
figure(2); 
plot(t1,SPC1(:,2),'--r','Linewidth',2); 
grid on; hold on; 
plot(t2,SPC2(:,2),'b','Linewidth',2); 
hold on; 
plot(t3,SPC3(:,2),'k','Linewidth',2); 
hold on; 
xlabel('time(seconds)');ylabel('Speed(r/s)');title('Plot for Speed vs Time'); 
legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 
figure(3); 
plot(t1,SPC1(:,3),'--r','Linewidth',2); 
grid on; hold on; 
plot(t2,SPC2(:,3),'b','Linewidth',2); 
hold on; 
plot(t3,SPC3(:,3),'k','Linewidth',2); 
hold on; 
xlabel('time(seconds)');ylabel('eqp');title('internal voltage along q-axis'); 
legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 
figure(4); 
plot(t1,SPC1(:,4),'--r','Linewidth',2); 
grid on; hold on; 
plot(t2,SPC2(:,4),'b','Linewidth',2); 
hold on; 
plot(t3,SPC3(:,4),'k','Linewidth',2); 
hold on; 
xlabel('time(seconds)');ylabel('efd)');title('SG field voltage along d-axis'); 
legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 
% ipfd and ipfq combined 
ipf1=(sqrt(SPC1(:,7).^2+SPC1(:,8).^2)); 
ipf2=(sqrt(SPC2(:,7).^2+SPC2(:,8).^2)); 
ipf3=(sqrt(SPC3(:,7).^2+SPC3(:,8).^2)); 
figure(5); 
plot(t1,ipf1,'--r','Linewidth',2); 
grid on; hold on; 
plot(t2,ipf2,'b','Linewidth',2); 
hold on; 
plot(t3,ipf3,'k','Linewidth',2); 
xlabel('time(seconds)');ylabel('ipf'); title('PV current before the filter'); 
legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 
 % ipd and ipq combined 
 figure(6); 
 plot(t1,(sqrt(SPC1(:,9).^2+SPC1(:,10)).^2),'--r','Linewidth',2); 
 grid on; hold on; 
 plot(t2,(sqrt(SPC2(:,9).^2+SPC2(:,10)).^2),'b','Linewidth',2); 
 hold on; 
 plot(t3,(sqrt(SPC3(:,9).^2+SPC3(:,10)).^2),'k','Linewidth',2); 
 xlabel('time(seconds)');ylabel('ip');title('PV current after the filter'); 
 legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
  
 % vcpd and vcpq combined 
 figure(7); 
 plot(t1,(sqrt(SPC1(:,11).^2+SPC1(:,12).^2)),'--r','Linewidth',2); 
 grid on; hold on; 
 plot(t2,(sqrt(SPC2(:,11).^2+SPC2(:,12).^2)),'b','Linewidth',2); 
 hold on; 
 plot(t3,(sqrt(SPC3(:,11).^2+SPC3(:,12).^2)),'k','Linewidth',2); 
 xlabel('time(seconds)');ylabel('Vop'); title('pv filter capacitor voltage'); 
 legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 % vdcs  
  figure(8); 
  plot(t1,SPC1(:,36),'--r','Linewidth',2); 
  grid on; hold on; 
  plot(t2,SPC2(:,36),'b','Linewidth',2); 
  hold on; 
  plot(t3,SPC3(:,36),'k','Linewidth',2); 
  xlabel('time(seconds)');ylabel('Vdcs');title('statcom dc link voltage'); 
  legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 % ipv  
  figure(9); 
  plot(t1,SPC1(:,5),'--r','Linewidth',2); 
  grid on; hold on; 
  plot(t2,SPC2(:,5),'b','Linewidth',2); 
  hold on; 
  plot(t3,SPC3(:,5),'k','Linewidth',2); 
  xlabel('time(seconds)');ylabel('ipv');title('PV array current'); 
   legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
 % vdcp  
%   figure(10); 
%   plot(t,SPC(:,6));grid on 
%   xlabel('time(seconds)');ylabel('Vdcp');title('PV dc link voltage'); 
%  pause; 
%   % iffd and iffq combined 
% iff=(sqrt(SPC(:,15).^2+SPC(:,16).^2)); 
% figure(11);plot(t,iff);grid on 
% xlabel('time(seconds)');ylabel('iff'); title('FC current before the filter'); 
% pause; 
% % ifd and ifq combined 
% figure(12);plot(t,(sqrt(SPC(:,17).^2+SPC(:,18)).^2));grid on 
% xlabel('time(seconds)');ylabel('ip');title('FC current after the filter'); 
% pause; 
% % vcfd and vcfq combined 
% figure(13);plot(t,(sqrt(SPC(:,19).^2+SPC(:,20).^2)));grid on 
% xlabel('time(seconds)');ylabel('Vop'); title('pv filter capacitor voltage'); 
% pause; 
% % ifc  
%  figure(14);plot(t,SPC(:,13));grid on 
%  xlabel('time(seconds)');ylabel('ifc');title('FC output current'); 
%  pause; 
% % vdcf  
%  figure(15);plot(t,SPC(:,14));grid on 
%  xlabel('time(seconds)');ylabel('Vdcf');title('FC dc link voltage'); 
%  pause; 
%   % igd and igq combined 
%  figure(16);plot(t,(sqrt(SPC(:,21).^2+SPC(:,22).^2)));grid on 
%  xlabel('time(seconds)');ylabel('ig');title('PMSG current'); 
%  pause; 
%  % Statcom current isd and isq combined 
%  figure(17);plot(t,(sqrt(SPC(:,34).^2+SPC(:,35).^2)));grid on 
%  xlabel('time(seconds)');ylabel('is');title('Statcom current'); 
%  pause; 
%  % vdcw 
%  figure(18);plot(t,SPC(:,27));grid on 
%  xlabel('time(seconds)');ylabel('Vdcw');title('wind turbine dc link voltage'); 
%  pause; 
%  % iwfd and iwfq combined 
%  figure(19);plot(t,(sqrt(SPC(:,28).^2+SPC(:,29).^2)));grid on 
%  xlabel('time(seconds)');ylabel('iwf');title('wind current before filter'); 
%  pause; 
%  % iwd and iwq combined 
%  figure(20);plot(t,(sqrt(SPC(:,30).^2+SPC(:,31).^2)));grid on 
%  xlabel('time(seconds)');ylabel('iw');title('wind current after filter'); 
%  pause; 
%  % vcwd and vcwq combined 
%  figure(20);plot(t,(sqrt(SPC(:,32).^2+SPC(:,33).^2)));grid on 
%  xlabel('time(seconds)');ylabel('Vow');title('wind filter capacitor voltage'); 
%  pause; 
 global values_MA  
 global values_PV 
 global values_Fc 
 global values_Wind 
  
 Efd0=values_MA(1); rt=values_MA(2); xt=values_MA(3);rb=values_MA(4); 
 xb=values_MA(5); vb=values_MA(6);g=values_MA(7);b=values_MA(8); 
 xd=values_MA(9);xq=values_MA(10);xdP=values_MA(11);Pm=values_MA(16); 
 rs=values_MA(17);Ls=values_MA(18); ms=values_MA(19);  
 %chiS=values_MA(20); 
   
Ldc=values_PV(1); Cdc=values_PV(2);Rpf=values_PV(3);Lpf=values_PV(4); 
Cpf=values_PV(5); Rpdr=values_PV(6); Rp=values_PV(7);Lp=values_PV(8); 
Ibasep=values_PV(9);phiP=values_PV(10);mp=values_PV(11); 
 
 
Lfc=values_Fc(1); Cfc=values_Fc(2);Rff=values_Fc(3);Lff=values_Fc(4); 
Cff=values_Fc(5);Rfdr=values_Fc(6);Rf=values_Fc(7);Lf=values_Fc(8); 
Ibasef=values_Fc(9);phiF=values_Fc(10);mf=values_Fc(11); 
 
Ks=values_Wind(1);Efdw=values_Wind(2);phiW=values_Wind(3); 
mW=values_Wind(4);mC=values_Wind(5);Rw=values_Wind(6); 
ra=values_Wind(7);Lw=values_Wind(8);Rwf=values_Wind(9); 
Lwf=values_Wind(10);Cwf=values_Wind(11);Rwdr=values_Wind(12); 
Pmw=values_Wind(15); 
 
%wind 
Xqw=0.7; Xdw=1;Hg =0.3;Ht =2;Cw=1;                        
%generator 
Ka=10; Ta=0.01;   M=1.2; w0=2*pi*60; wo=w0; 
TdoP=5;  Vto=1.032; 
 
%  Auxillary constants 
 X1=xdP+xt; 
 x2=xq+xt; 
 z1=rt^2+X1*x2; 
 zb=rb^2+xb^2; 
  
 A=g*zb*z1+z1*rb+rt*zb; 
 B=b*zb*z1+xb*z1+x2*zb; 
 C=z1*(xb*x2+b*zb*x2+zb)-rt^2*zb; 
 Den=B*C+A^2*x2; 
 D=C*zb*z1; 
 E=C*zb*x2+A*x2*rt*zb; 
 F=A*x2*z1*zb; 
  
 G=C*z1*(rb*sin(SPC1(:,1))+xb*cos(SPC1(:,1)))+A*x2*z1*(rb*cos(SPC1(:,1))-xb*sin(SPC1(:,1))); 
 % Unstable response 
 Vsq=(D*(SPC1(:,9)+SPC1(:,17)+SPC1(:,30)+SPC1(:,34))+F*(SPC1(:,10)+SPC1(:,18)+SPC1(:,31)+SPC1(:,35))+E*SPC1(:,3)+G*vb)/Den; 
 Vsd=(z1*zb*(SPC1(:,9)+SPC1(:,17)+SPC1(:,30)+SPC1(:,34))+zb*x2*SPC1(:,3)+z1*vb*(rb*sin(SPC1(:,1))+xb*cos(SPC1(:,1)))-B*Vsq)/A; 
 Vs=sqrt(Vsd.^2+Vsq.^2); 
 
% Micro Alternator power 
 itd=((SPC1(:,3)-Vsq)*x2-rt*Vsd)/z1; 
 itq=(Vsd*z1-rt^2*Vsd+rt*(SPC1(:,3)-Vsq)*x2)/(x2*z1); 
 Vd=xq*itq; 
 Vq=SPC1(:,3)-xdP*itd; 
 Pe=(Vd.*itd+Vq.*itq); 
 Vt=(sqrt(Vd.*Vd+Vq.*Vq)); 
  
% Wind system power calculation 
  vgq=mC*SPC1(:,27).*cos(SPC1(:,23)); 
  vgd=mC*SPC1(:,27).*sin(SPC1(:,23)); 
  Vg=sqrt(vgd.^2+vgq.^2); 
   
  Pewc= vgd.*SPC1(:,21)+vgq.*SPC1(:,22); 
  Vbase=40;G=1;TaC=25; 
  if t1>=1 
      Vactual = bp_sx150svoltcalc(0*Ibasep,G,TaC); 
  else  
  Vactual = bp_sx150svoltcalc(SPC1(:,5)*Ibasep,G,TaC); 
  end 
  Vpv=Vactual/Vbase; 
  Ppv=Vpv.*SPC1(:,5); 
   
  % Fuel cell power calculation 
  Vbasef=50; 
  if t1>=1 
      Vactualf = fuelcell_volt_calc(0*Ibasef); 
  else  
  Vactualf = fuelcell_volt_calc(SPC1(:,13)*Ibasef); 
  end 
  Vfc=Vactualf/Vbasef; 
  Pfc=Vfc.*SPC1(:,13);   
   
 % Stable response chist controlled 
  G2=C*z1*(rb*sin(SPC2(:,1))+xb*cos(SPC2(:,1)))+A*x2*z1*(rb*cos(SPC2(:,1))-xb*sin(SPC2(:,1))); 
 % Unstable response 
 Vsq2=(D*(SPC2(:,9)+SPC2(:,17)+SPC2(:,30)+SPC2(:,34))+F*(SPC2(:,10)+SPC2(:,18)+SPC2(:,31)+SPC2(:,35))+E*SPC2(:,3)+G2*vb)/Den; 
 Vsd2=(z1*zb*(SPC2(:,9)+SPC2(:,17)+SPC2(:,30)+SPC2(:,34))+zb*x2*SPC2(:,3)+z1*vb*(rb*sin(SPC2(:,1))+xb*cos(SPC2(:,1)))-B*Vsq2)/A; 
 Vs2=sqrt(Vsd2.^2+Vsq2.^2); 
  
% Micro Alternator power 
 itd2=((SPC2(:,3)-Vsq2)*x2-rt*Vsd2)/z1; 
 itq2=(Vsd2*z1-rt^2*Vsd2+rt*(SPC2(:,3)-Vsq2)*x2)/(x2*z1); 
 Vd2=xq*itq2; 
 Vq2=SPC2(:,3)-xdP*itd2; 
 Pe2=(Vd2.*itd2+Vq2.*itq2); 
 Vt2=(sqrt(Vd2.*Vd2+Vq2.*Vq2)); 
  
% Wind system power calculation 
  vgq2=mC*SPC2(:,27).*cos(SPC2(:,23)); 
  vgd2=mC*SPC2(:,27).*sin(SPC2(:,23)); 
  Vg2=sqrt(vgd2.^2+vgq2.^2); 
   
  Pewc2= vgd2.*SPC2(:,21)+vgq2.*SPC2(:,22); 
  Vbase=40;G=1;TaC=25; 
  if t2>=1 
      Vactual = bp_sx150svoltcalc(0*Ibasep,G,TaC); 
  else  
  Vactual = bp_sx150svoltcalc(SPC2(:,5)*Ibasep,G,TaC); 
  end 
  Vpv2=Vactual/Vbase; 
  Ppv2=Vpv2.*SPC2(:,5); 
   
  % Fuel cell power calculation 
  Vbasef=50; 
  if t2>=1 
      Vactualf = fuelcell_volt_calc(0*Ibasef); 
  else  
  Vactualf = fuelcell_volt_calc(SPC2(:,13)*Ibasef); 
  end 
  Vfc2=Vactualf/Vbasef; 
  Pfc2=Vfc2.*SPC2(:,13);   
   
 % Chip controlled stable response 
 G3=C*z1*(rb*sin(SPC3(:,1))+xb*cos(SPC3(:,1)))+A*x2*z1*(rb*cos(SPC3(:,1))-xb*sin(SPC3(:,1))); 
  
 Vsq3=(D*(SPC3(:,9)+SPC3(:,17)+SPC3(:,30)+SPC3(:,34))+F*(SPC3(:,10)+SPC3(:,18)+SPC3(:,31)+SPC3(:,35))+E*SPC3(:,3)+G3*vb)/Den; 
 Vsd3=(z1*zb*(SPC3(:,9)+SPC3(:,17)+SPC3(:,30)+SPC3(:,34))+zb*x2*SPC3(:,3)+z1*vb*(rb*sin(SPC3(:,1))+xb*cos(SPC3(:,1)))-B*Vsq3)/A; 
 Vs3=sqrt(Vsd3.^2+Vsq3.^2); 
  
% Micro Alternator power 
 itd3=((SPC3(:,3)-Vsq3)*x2-rt*Vsd3)/z1; 
 itq3=(Vsd3*z1-rt^2*Vsd3+rt*(SPC3(:,3)-Vsq3)*x2)/(x2*z1); 
 Vd3=xq*itq3; 
 Vq3=SPC3(:,3)-xdP*itd3; 
 Pe3=(Vd3.*itd3+Vq3.*itq3); 
 Vt3=(sqrt(Vd3.*Vd3+Vq3.*Vq3)); 
  
% Wind system power calculation 
  vgq3=mC*SPC3(:,27).*cos(SPC3(:,23)); 
  vgd3=mC*SPC3(:,27).*sin(SPC3(:,23)); 
  Vg3=sqrt(vgd3.^2+vgq3.^2); 
   
  Pewc3= vgd3.*SPC3(:,21)+vgq3.*SPC3(:,22); 
  Vbase=40;G=1;TaC=25; 
  if t3>=1 
      Vactual = bp_sx150svoltcalc(0*Ibasep,G,TaC); 
  else  
  Vactual = bp_sx150svoltcalc(SPC3(:,5)*Ibasep,G,TaC); 
  end 
  Vpv3=Vactual/Vbase; 
  Ppv3=Vpv3.*SPC3(:,5); 
   
  % Fuel cell power calculation 
  Vbasef=50; 
  if t3>=1 
      Vactualf = fuelcell_volt_calc(0*Ibasef); 
  else  
  Vactualf = fuelcell_volt_calc(SPC3(:,13)*Ibasef); 
  end 
  Vfc3=Vactualf/Vbasef; 
  Pfc3=Vfc3.*SPC3(:,13);   
 
 % Plotting the corresponding waveforms 
 figure(21); 
 plot(t1,Vs,'--r'); 
 grid on;hold on; 
 plot(t2,Vs2,'b'); 
 hold on; 
 plot(t3,Vs3,'k'); 
xlabel('time(seconds)');ylabel('Vs');title('Microgrid voltage'); 
legend('Uncontrolled','\psi_{st}Controlled','\psi_{p}Controlled'); 
  
 figure(22); 
 plot(t1,Pewc+Pe+Ppv+Pfc,'--r'); 
 grid on; hold on; 
 plot(t2,Pewc2+Pe2+Ppv2+Pfc2,'b'); 
 hold on; 
 plot(t3,Pewc3+Pe3+Ppv3+Pfc3,'k'); 
  xlabel('time(seconds)');ylabel('P');title('Total real power output');