short-circuit fault newton-raphson power-flow

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)

Thanks in advance!

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
BD(faultbus,7:8)=0; %removing the load
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

Answer 1

[Sambit Senapati]

You might want to take a look in to the following link:

https://www.mathworks.com/matlabcentral/fileexchange/21059-newton-raphson-loadflow