# short-circuit fault newton-raphson power-flow

4 views (last 30 days)
Can anyone explain to me, how to write a MATLAB code to solve for Newton-Raphson power-flow during a short-circuit fault? Like, as far as I know, I need to change the Y-bus matrix (remove the row & column corresponding to that bus) then modify the Jacobian and the delta P&Q matrix too and then calculate for the voltage and angle for the next iteration. If I program like this, the program isn't converging. Can you guys tell me if I'm missing/wrong about something?
If someone could point me to a resource, that would be okay too.
(P.S: I'm not considering the sequence networks. Just the normal NR powerflow)
Edit: I've included my code below.
%% Short-circuit Power Flow
%-----------------------------------------------------------------------
% B# code |V| Theta PG QG PL QL Qmin Qmax
BD = [ 1 1 1.04 0.0 0 0 0 0 0 0
2 2 1.025 0.0 1.63 0 0 0 0 0
3 2 1.025 0.0 0.85 0 0 0 0 0
4 3 1.0 0.0 0 0 0 0 0 0
5 3 1.0 0.0 0 0 1.25 0.5 0 0
6 3 1.0 0.0 0 0 0.9 0.3 0 0
7 3 1.0 0.0 0 0 0 0 0 0
8 3 1.0 0.0 0 0 1.0 0.35 0 0
9 3 1.0 0.0 0 0 0 0 0 0];
% ----------------------------------------------------------------------
% Line code
% = 1 for no transformer
% tb fb R X B/2 > 1 or < 1 for tr. tap at bus nl
LD = [ 1 4 0.0 0.0576 0.0 1
2 7 0.0 0.0625 0.0 1
3 9 0.0 0.0586 0.0 1
4 5 0.01 0.085 0.0880 1
4 6 0.017 0.092 0.0790 1
5 7 0.032 0.161 0.153 1
6 9 0.039 0.17 0.179 1
7 8 0.0085 0.072 0.0745 1
8 9 0.0119 0.1008 0.1045 1];
faultbus=5;
BD(2,5)=0; BD(3,5)=0; %bring the generation to zero
LD(4,3:4)=LD(4,3:4)*(10^-5); %Setting the line resistance & reactance very low
LD(6,3:4)=LD(6,3:4)*(10^-5);
Vm(faultbus)=0; Va(faultbus)=0; % setting Vmagnitude & Vangle to zero
N=max(BD(:,1)); % Find the no. of buses (N)
PG=BD(:,5); QG=BD(:,6); PL=BD(:,7); QL=BD(:,8); % Powers of generator & Loads
% Things related to Y-matrix
[Y]=Ybus(LD);
G=real(Y); B=imag(Y); Ym=abs(Y); Ya=angle(Y);
%initializing variables
it=0; maxit = 20; errmax=.0001; err=2*errmax; % Start the iterative process
P(N,1)=0; Q(N,1)=0; J1=zeros(N,N); J2=zeros(N,N); J3=zeros(N,N); J4=zeros(N,N);
% Iterative loop starts here:
while err>errmax & it < maxit; it=it+1;
% Forming the non-diag elements J11,J12,J21,&J22
for i=1:N; for j=1:N; if j~=i
J1(i,j)= +Vm(i)*Vm(j)*Ym(i,j)*sin(Va(i)-Va(j)-Ya(i,j));
J2(i,j)= +Vm(i)*Ym(i,j)*cos(Va(i)-Va(j)-Ya(i,j));
J3(i,j)= -Vm(i)*Vm(j)*Ym(i,j)*cos(Va(i)-Va(j)-Ya(i,j));
J4(i,j)= +Vm(i)*Ym(i,j)*sin(Va(i)-Va(j)-Ya(i,j)); end; end; end;
% Forming the diag elements J11,J12,J21,&J22
for i=1:N ; J1(i,i)=0; J2(i,i)= 2*Vm(i)*G(i,i); J3(i,i)=0; J4(i,i)=-2*Vm(i)*B(i,i);
for j=1:N; if j~=i
J1(i,i)=J1(i,i)-J1(i,j); J2(i,i)=J2(i,i)-J3(i,j)/Vm(i);
J3(i,i)=J3(i,i)-J3(i,j); J4(i,i)=J4(i,i)+J1(i,j)/Vm(i); end; end;
P(i)= Vm(i)^2 * G(i,i)+J3(i,i); Q(i)=-Vm(i)^2 * B(i,i)-J1(i,i); end
J=[J1 J2; J3 J4]; % Full Jacobian Matrix
% Define a vector consisting of the PV bus numbers ; Call this vector PV
k=1; for i=1:N; if BD(i,2)==2 ; PV(k)=BD(i,1); k=k+1; end ; end
J([faultbus N+faultbus],:)=[]; J(:,[faultbus N+faultbus])=[]; %modify the Jacobean; remove rows and columns corresponding to faultbus
% Finding Delta P & Q
dP=(PG-PL)-P;
dQ=(QG-QL)-Q;
dP(faultbus)=[]; dQ(faultbus)=[]; %remove P&Q for faultbus
dPQ=[dP;dQ];
% Update the values voltage magnitudes/angles
err=max(abs(dPQ)); dVmVa=(J)\dPQ;
Va=Va+[dVmVa(1:faultbus-1);0;dVmVa(faultbus:N-1)];
Va=Va+[dVmVa(N:N+faultbus-1);0;dVmVa(N+faultbus:2*N-2)];
end
it % Show the no of iteration for convergence
%
if it >= maxit; warning('Loadflow Solution Failed to Converge'); end
%printing the output
fprintf('\nBus\t |V|\t Angle(deg)\t P\t\t\t Q\t\n');
fprintf('----------------------------------------------------\r');
for i = 1:1:N
fprintf('%2d\t %6.3f\t %10.3f\t %10.3f\t %10.3f\n',i, Vm(i), Va(i)*180/pi, P(i), Q(i)');
end

Sambit Senapati on 14 Jun 2019
You might want to take a look in to the following link: