% Program For Newton Raphson Load Flow Analsis ..
nbus=5; %IEEE - 5
Y= ybusppg1s(nbus); %Calling ybusppg.m to get Y- Bus Matrix
%Type
%1- Slack Bus ..
%2- PV Bus ..
%3- PQ Bus ..
busdata=callbusdatas(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
loadflow1s(nbus,V,del,BMva);
评论1