upfc.rar_FACT_Transmission_UPFC_UPFC FACT_upfc matlab

% Program For Newton Raphson Load Flow Analsis .. nbus=5; %IEEE - 5 Y=ybusppg1u(nbus); %Calling ybusppg.m to get Y- Bus Matrix %Type %1- Slack Bus .. %2- PV Bus .. %3- PQ Bus .. busdata=callbusdatau(nbus); BMva=100; bus = busdata (:,1); type=busdata(:,2); V=busdata(:,3); del=busdata(:,4); Pg=busdata(:,5)/BMva; Qg=busdata(:,6)/BMva; P1=busdata(:,7)/BMva; Q1=busdata(:,8)/BMva; Qmin=busdata(:,9)/BMva; Qmax=busdata(:,10)/BMva; P=Pg-P1; Q=Qg-Q1; Psp=P; Qsp=Q; G=real(Y); B=imag(Y); pv=find(type==2|type==1); pq= find(type==3); npv=length(pv); npq=length(pq); %PQ Buses .. Tol=1; Iter=1; while (Tol>1e-5) P=zeros(nbus,1); Q=zeros(nbus,1); %Calculate P and Q for i=1:nbus for k=1:nbus P(i)=P(i)+V(i)*V(k)*(G(i,k)*cos(del(i)-del(k))+B(i,k)*sin(del(i)-del(k))); Q(i)=Q(i)+V(i)*V(k)*(G(i,k)*sin(del(i)-del(k))-B(i,k)*cos(del(i)-del(k))); end end %Checking Q-limit violation .. if Iter<=7&&Iter>2 for n= 2:nbus if type(n)==2 QG=Q(n)+Q1(n); if QG<Qmin(n) V(n)=V(n)+0.01; elseif QG>Qmax(n) V(n)=V(n)-0.01; end end end end % Calculate chanage from specified value dPa=Psp-P; dQa=Qsp-Q; k=1; dQ=zeros(npq,1); for i=1:nbus if type(i)==3 dQ(k,1)=dQa(i); k=k+1; end end dP=dPa(2:nbus); M=[dP;dQ]; % Mismatch Vector % Jacobian .. % J1 - Derivative of Real Power Injection with Angle .. J1=zeros(nbus-1,nbus-1); for i=1:(nbus-1) m=i+1; for k=1:(nbus-1) n=k+1; if n==m for n=1:nbus J1(i,k)=J1(i,k)+V(m)*V(n)*(-G(m,n)*sin(del(m)-del(n))+B(m,n)*cos(del(m)-del(n))); end J1(i,k)=J1(i,k)-V(m)^2*B(m,m); else J1(i,k)=V(m)*V(n)*(G(m,n)*sin(del(m)-del(n))-B(m,n)*cos(del(m)-del(n))); end end end %J2- Derivative of Real power Injection with V .. J2=zeros(nbus-1,npq); for i= 1:(nbus-1) m=i+1; for k=1:npq n=pq(k); if n==m for n=1:nbus J2(i,k)=J2(i,k)+V(n)*(G(m,n)*cos(del(m)-del(n))+B(m,n)*sin(del(m)-del(n))); end J2(i,k)=J2(i,k)+V(m)*G(m,m); else J2(i,k)=V(m)*(G(m,n)*cos(del(m)-del(n))+B(m,n)*sin(del(m)-del(n))); end end end % J3- Derivative of Reactive Power Injection with Angle .. J3=zeros(npq,nbus-1); for i=1:npq m=pq(i); for k=1:(nbus-1) n=k+1; if n==m for n=1:nbus J3(i,k)=J3(i,k)+V(m)*V(n)*(G(m,n)*cos(del(m)-del(n))+B(m,n)*sin(del(m)-del(n))); end J3(i,k)=J3(i,k)-V(m)^2*G(m,m); else J3(i,k)= V(m)*V(n)*(-G(m,n)*cos(del(m)-del(n))-B(m,n)*sin(del(m)-del(n))); end end end % J4- Derivative of Reactive Power Injection with V .. J4=zeros(npq,npq); for i=1:npq m=pq(i); for k=1:npq n=pq(k); if n==m for n=1:nbus J4(i,k)=J4(i,k)+V(n)*(G(m,n)*sin(del(m)-del(n))-B(m,n)*cos(del(m)-del(n))); end J4(i,k)=J4(i,k)-V(m)*B(m,m); else J4(i,k)=V(m)*(G(m,n)*sin(del(m)-del(n))-B(m,n)*cos(del(m)-del(n))); end end end J=[J1 J2; J3 J4]; X=inv(J)*M; dTh=X(1:nbus-1); dV=X(nbus:end); %Updating State Vector .. del(2:nbus)=dTh+del(2:nbus); k=1; for i=2:nbus if type(i)==3 V(i)=dV(k)+V(i); % Voltage Magnitude k=k+1; end end Iter=Iter+1; Tol=max(abs(M)); end loadflow1u(nbus,V,del,BMva);