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


% Non-linear model with statcom and composite microgrid 
% Total no of states are 36 
% Statcom will be acting as centralized controller in the model 
% Author: Mohammad Ashraf Ali @4th_Nov_2012 
 
function xp=non_linear_36states_chip_controlled_pi_controller(t,x) 
 
 global values_MA  
 global values_PV 
 global values_Fc 
 global values_Wind 
 
% Assigning the important variables or parameters 
 
 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); 
  
%%%%% Parameters%%%%%%  
%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;  
 TdoP=5;  vto=1.032; Cdcs=2.0; 
% controller parameters  
 kpc=40;kic=300; 
  
% define constants 
x1=xdP+xt; 
x2=xq+xt; 
z1=rt^2+x1*x2; 
zb=rb^2+xb^2; 
A=g*zb*z1+rb*z1+rt*zb; 
B=b*zb*z1+x2*zb+xb*z1; 
C=x1*zb+b*zb*z1+xb*z1; 
D=rt*zb+g*zb*z1+rb*z1; 
Den=B*C+A*D; 
% vsq coefficients 
G=C*zb*z1/Den; 
H=A*zb*z1/Den; 
I=(C*x2*zb+A*rt*zb)/Den; 
K=((A*rb*z1+C*xb*z1)*vb)/Den; 
L=((C*rb*z1-A*xb*z1)*vb)/Den; 
 
 
vsq=G*(x(9)+x(17)+x(30)+x(34))+H*(x(10)+x(18)+x(31)+x(35))+I*x(3)+K*cos(x(1))+L*sin(x(1)); 
 
% vsd coefficients 
G1=(zb*z1-B*G)/A; 
H1=(-B*H)/A; 
I1=(x2*zb-B*I)/A; 
K1=(xb*z1*vb-B*K)/A; 
L1=(rb*z1*vb-B*L)/A; 
 
vsd=G1*(x(9)+x(17)+x(30)+x(34))+H1*(x(10)+x(18)+x(31)+x(35))+I1*x(3)+K1*cos(x(1))+L1*sin(x(1)); 
 
vs=sqrt(vsd*vsd+vsq*vsq); 
 
% itd and itq equations and power equations 
itd=(-rt*vsd+x2*(x(3)-vsq))/z1; 
itq=(x1*vsd+rt*(x(3)-vsq))/z1; 
 
vd=xq*itq; vq=x(3)-xdP*itd; 
Pe=vd*itd+vq*itq; 
vt=sqrt(vd*vd+vq*vq); 
 
ipf=sqrt(x(7)^2+x(8)^2);   
iff=sqrt(x(15)^2+x(16)^2); 
 
xp=zeros(36,1); 
% Microalternator non-linear equations 
phiP=kpc*(1-x(2))+x(37); % x(37) is the chiS State 
xp(1)=w0*(x(2)-1); 
 
if t>=1 && t<= 1.2 
      disturbance=0.1; 
     xp(2)=(1/M)*(Pm-Pe+disturbance); 
      else 
    xp(2)=(1/M)*(Pm-Pe); 
  end 
 
xp(3)=(1/TdoP)*(x(4)-x(3)-(xd-xdP)*itd); 
 
xp(4)=(Ka/Ta)*(vto-vt)-(1/Ta)*(x(4)-Efd0); 
 
% PV system Non-linear Equations 
  Vbase=40;G=1; TaC=25; 
  Vactualp = bp_sx150svoltcalc(x(5)*Ibasep,G,TaC); 
  vpv=Vactualp/Vbase; 
  dcp=1-0.861411783523035/1.000392156862745; 
   
xp(5)=(1/Ldc)*(vpv-(1-dcp)*x(6)); % ipv equation 
xp(6)=(1/Cdc)*((1-dcp)*x(5)-mp*x(7)*cos(phiP)-mp*x(8)*sin(phiP));%vdcp equation 
% ipfd and ipfq equations 
xp(7)=(-w0*Rpf/Lpf)*x(7)+w0*x(2)*x(8)+(w0*mp/Lpf)*x(6)*cos(phiP)+(-w0/Lpf)*x(11)+(-w0*Rpdr/Lpf)*(x(7)-x(9)); 
xp(8)=(-w0*Rpf/Lpf)*x(8)-w0*x(2)*x(7)+(w0*mp/Lpf)*x(6)*sin(phiP)+(-w0/Lpf)*x(12)+(-w0*Rpdr/Lpf)*(x(8)-x(10)); 
% ipd and ipq equations 
xp(9)=(-w0*Rp/Lp)*x(9)+w0*x(2)*x(10)+(w0/Lp)*(x(11)-vsd)+(w0*Rpdr/Lp)*(x(7)-x(9)); 
xp(10)=(-w0*Rp/Lp)*x(10)-w0*x(2)*x(9)+(w0/Lp)*(x(12)-vsq)+(w0*Rpdr/Lp)*(x(8)-x(10)); 
%vcpd and vcpq equations 
xp(11)=(w0/Cpf)*(x(7)-x(9))+w0*x(2)*x(12); 
xp(12)=(w0/Cpf)*(x(8)-x(10))-w0*x(2)*x(11); 
% Fuel cell model  
  Vbase=50; 
  Vactualf = fuelcell_volt_calc(x(13)*Ibasef); 
  vfc=Vactualf/Vbase; 
  dcf=1-(0.9/1.000392156862745); 
     
%Fuel cell generation system non-linear equations. 
xp(13)=(1/Lfc)*(vfc-(1-dcf)*x(14)); % ifc equation 
xp(14)=(1/Cfc)*((1-dcf)*x(13)-mf*x(15)*cos(phiF)-mf*x(16)*sin(phiF)); %vcfd equation 
% iffd and iffq equations 
xp(15)=(-w0*Rff/Lff)*x(15)+w0*x(2)*x(16)+(w0*mf/Lff)*x(14)*cos(phiF)+(-w0/Lff)*x(19)+(-w0*Rfdr/Lff)*(x(15)-x(17)); 
xp(16)=(-w0*Rff/Lff)*x(16)-w0*x(2)*x(15)+(w0*mf/Lff)*x(14)*sin(phiF)+(-w0/Lff)*x(20)+(-w0*Rfdr/Lff)*(x(16)-x(18)); 
%ifd and ifq equations 
xp(17)=(-w0*Rf/Lf)*x(17)+w0*x(2)*x(18)+(w0/Lf)*(x(19)-vsd)+(w0*Rfdr/Lf)*(x(15)-x(17)); 
xp(18)=(-w0*Rf/Lf)*x(18)-w0*x(2)*x(17)+(w0/Lf)*(x(20)-vsq)+(w0*Rfdr/Lf)*(x(16)-x(18)); 
% vcfd and vcfq equations 
xp(19)=(w0/Cff)*(x(15)-x(17))+w0*x(2)*x(20); 
xp(20)=(w0/Cff)*(x(16)-x(18))-w0*x(2)*x(19); 
 
% wind system 
vgd=mC*x(27)*sin(x(23)); 
vgq=mC*x(27)*cos(x(23)); 
vg=sqrt(vgd*vgd+vgq*vgq); 
Pew=(x(21)^2+x(22)^2)*ra+vgd*x(21)+vgq*x(22); 
 
%PMSG output currents(igd & igq) 
xp(21)=(-w0*ra/xdw)*x(21)+(w0*xqw/xdw)*x(22)*x(24)-(w0/xdw)*vgd; 
xp(22)=(-w0*ra/xqw)*x(22)-(w0*xdw/xqw)*x(24)*x(21)+(w0*Efdw/xqw)*x(24)-(w0/xqw)*vgq; 
 
% PMSG Swing equations 
xp(23)=w0*(x(24)-1); %delta w equations 
 
xp(24)=(1/(2*Hg))*(Ks*x(25)-Pew-0.7*(x(24)-1));% ww equation 
 
xp(25)=w0*(x(26)-x(24)); % theta s equation 
 
xp(26)=(1/(2*Ht))*(Pmw-Ks*x(25)-0.7*(x(26)-1));%wt equation 
 
% wind system currents and voltages 
% vdcw equation 
xp(27)=(1/Cw)*(mC*x(21)*sin(x(23))+mC*x(22)*cos(x(23))-mW*x(28)*cos(phiW)-mW*x(29)*sin(phiW)); 
 
% iwfd and iwfq equation 
xp(28)=(-w0*Rwf/Lwf)*x(28)+w0*x(24)*x(29)+(w0*mW*cos(phiW)/Lwf)*x(27)+(-w0/Lwf)*x(32)+(-w0*Rwdr/Lwf)*(x(28)-x(30)); 
 
xp(29)=(-w0*Rwf/Lwf)*x(29)-w0*x(24)*x(28)+(w0*mW*sin(phiW)/Lwf)*x(27)+(-w0/Lwf)*x(33)+(-w0*Rwdr/Lwf)*(x(29)-x(31)); 
% iwd and iwq equation 
xp(30)=(-w0*Rw/Lw)*x(30)+w0*x(24)*x(31)+(w0/Lw)*(x(32)-vsd)+(w0*Rwdr/Lw)*(x(28)-x(30)); 
 
xp(31)=(-w0*Rw/Lw)*x(31)-w0*x(24)*x(30)+(w0/Lw)*(x(33)-vsq)+(w0*Rwdr/Lw)*(x(29)-x(31)); 
% vcwd and vcwq equation 
xp(32)=(w0/Cwf)*(x(28)-x(30))+w0*x(24)*x(33); 
 
xp(33)=(w0/Cwf)*(x(29)-x(31))-w0*x(24)*x(32); 
% Battery parameters 
Vbatt=0.964557519964956; Rbatt=0.01; 
% Statcom dynamic equations 
xp(34)=(-w0*rs/Ls)*x(34)+w0*x(2)*x(35)+(w0*ms*cos(chiS)/Ls)*x(36)+(-w0/Ls)*vsd; 
xp(35)=(-w0*rs/Ls)*x(35)-w0*x(2)*x(34)+(w0*ms*sin(chiS)/Ls)*x(36)+(-w0/Ls)*vsq; 
xp(36)=(-1/Cdcs)*(x(34)*ms*cos(chiS)+x(35)*ms*sin(chiS))+(Vbatt-x(36))/(Rbatt*Cdcs); 
xp(37)=kic*(1-x(2)); 
end