short-circuit fault newton-raphson power-flow

10 vues (au cours des 30 derniers jours)
Joseph Nikhil Raj Koppula
Joseph Nikhil Raj Koppula le 1 Juin 2019
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

Réponses (1)

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

Catégories

En savoir plus sur Mathematics dans Help Center et File Exchange

Produits


Version

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by