Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@battulaharish2
battulaharish2 authored Sep 29, 2022
1 parent 145b5d9 commit b373370
Show file tree
Hide file tree
Showing 2 changed files with 377 additions and 0 deletions.
211 changes: 211 additions & 0 deletions NR_adjQ.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,211 @@
clc;
clear;
%% data
load('busd14');
load('ldata14');
%% ybus formation
l = ldata14(:,1);
nline=length(l);
nbus=length(busd14(:,1));
fb = ldata14(:,2);
tb = ldata14(:,3);
R = ldata14(:,4);
X = ldata14(:,5);
hc = ldata14(:,6);
t = ldata14(:,7);
Z = complex(R,X);
y=(1./Z);
sb1=max(fb);
sb2=max(tb);
bus = max(sb1,sb2);
ybus=zeros(bus,bus);
for k=1:nline
p=fb(k);
q=tb(k);
ybus(p,p)=ybus(p,p)+y(k)+1i*hc(k);
ybus(q,q)=ybus(q,q)+y(k)+1i*hc(k);
ybus(p,q)=ybus(p,q)-y(k);
ybus(q,p)=ybus(p,q);
end
Y=ybus;
%% reading bus data
BMVA=100;
bus=busd14(:,1);
type=busd14(:,2);
v=busd14(:,3);
dell=busd14(:,4);
del=dell*pi/180;
Pg=busd14(:,5)/BMVA;
Qg=busd14(:,6)/BMVA;
Pl=busd14(:,7)/BMVA;
Ql=busd14(:,8)/BMVA;
Qmin=busd14(:,9)/BMVA;
Qmax=busd14(:,10)/BMVA;
P=Pg-Pl;
Q=Qg-Ql;
Psp=P;
Qsp=Q;
G=real(Y);
B=imag(Y);
PQ=find(type==3);
npq=length(PQ);
Tol=1e-3;
err = 1;
Iter=0;
vv=0;
%% NR load flow
st=clock; % start the iteration time clock
while err>Tol
P=zeros(nbus,1);
Q=zeros(nbus,1);
%caluclate 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
for i=2:nbus
if type(i) == 2
if Q(i) < Qmin(i) % Whether violated the lower limit.
Q(i) = Qmin(i);
type(i) = 3;
vv = i;
end
if Q(i) > Qmax(i) % No, violated the upper limit
Q(i) = Qmax(i);
type(i) = 3; % If Violated, change PV bus to PQ bus.
vv = i;
end
end
type;
PQ=find(type==3);
npq=length(PQ);
end

%calculate MISMATCH: change 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];
%jacobian
%J1- dervative of real power injection with angles
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)^2*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 derviative of reactive of reactive power injection with
%angles
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 reative 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(m)*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)^2*B(m,m);
else
J4(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
J=[J1 J2; J3 J4];
x=J\M;
dth=x(1:nbus-1);
dv=x(nbus:end);
%updating state vectors
del(2:nbus)=dth+del(2:nbus);
k=1;
for i=2:nbus
if type(i)==3
v(i)=dv(k)+v(i);
k=k+1;
end
end
for i=1:nbus
E(i)=v(i)*(cos(del(i))+sin(del(i))*1j);
end
Iter=Iter+1;
err=max(abs(M));
end
ste=clock; % stop the iteration time clock
lf = [type abs(E)' rad2deg(angle(E))' P.*100 Q.*100];
%% printing results
disp('--------------------------------------------------------------------------')
disp(' Newton Rapson Load-Flow Study');
disp('--------------------------------------------------------------------------')
disp(date);
fprintf('Number of iterations : %d \n', Iter);
fprintf('Error tolerance considered : %d \n', Tol);
fprintf('Solution time : %g sec.\n',etime(ste,st));
fprintf('Maximum mismatch : %d\n\n',max(abs(M)));
% fprintf('Jacobian Matrix:\n');
% disp('==================');
% disp(J);
% fprintf('BUS Data:\n');
% disp('============');
% disp(' bus type v delta Pg Qg Pl Q1 Qmin Qmax');
% disp(' ----- ----- --- ----- ---- ---- ----- --- ------ ------')
% disp(busd);
% fprintf('Load Flow:\n');
% disp('=============');
% fprintf('Q-limit voilated at bus: %d\n\n', vv);
disp(' Bus Type Vpu angle P(MW) Q(MVAr)')
disp(' -------- -------- -------- -------- -------');
disp(lf)
power_loss = sum(lf(:,4));
fprintf('Power loss:%d MW\n \n',power_loss);
166 changes: 166 additions & 0 deletions fdlf.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,166 @@
clc
clear all
% sl sb rb R x b tap
linedata=[1 1 2 0.042 0.168 0.082 1;
2 1 5 0.031 0.126 0.062 1;
3 2 3 0.031 0.126 0.062 1;
4 3 4 0.024 0.136 0.164 1;
5 3 5 0.053 0.210 0.102 1;
6 4 5 0.063 0.252 0.122 1];
% bus type V d Pg Qg Pl Ql Qmin() Qmax
busdata = [ 1 1 1.00 0 0.0 0 0 0 0 0;
2 2 1.00 0 50 0 0 0 0 500;
3 2 1.00 0 100 0 0 0 -500 500;
4 3 1.00 0 0.0 0 115 60 -500 500;
5 3 1.00 0 0.0 0 85 40 -500 500;
];
nline = linedata(:,1);
sb = linedata(:,2); %Sending End Bus column
rb = linedata(:,3); %Receiving End Bus column
R = linedata(:,4); %Resistance(pu) in lines
X = linedata(:,5); %Reactance(pu) in lines
hc = 1i*linedata(:,6); %Capacitive Susceptance(pu)
t = linedata(:,7); %nominal tap ratio
Z = R + 1i*X;y=(1./Z);
sb1=max(sb);sb2=max(rb);
bus = max(sb1,sb2); %number of buses
line = length(nline); %number of lines
ybus=zeros(bus,bus);
for k=1:line %loop for tap ratio
if t(k)~=1
t1 = (1-1/t(k));
t2 = -t1/t(k);
hc(k) = t2*y(k);
y(k) = y(k)/t(k);
end
end
for k=1:line %loop for Ybus
p=sb(k);q=rb(k);
ybus(p,p)=ybus(p,p)+y(k)+hc(k)/2;
ybus(q,q)=ybus(q,q)+y(k)+hc(k)/2;
ybus(p,q)=ybus(p,q)-y(k);
ybus(q,p)=ybus(p,q);
end
busnumber=busdata(:,1); %Bus number
type = busdata(:,2);%bus type
Vmag = busdata(:,3);%Voltage magnitude
V = busdata(:,3).*exp(1i*busdata(:,4));
ANG = busdata(:,4);%bus angle
PG = busdata(:,5);%generated real power
QG = busdata(:,6);%generated reactive power
PL = busdata(:,7);%load real power
QL = busdata(:,8);%load reactive power
BaseMVA =100;
Load = PL + 1i*QL;
Qmin = busdata(:,9)/100; % Minimum Reactive Power Limit
Qmax = busdata(:,10)/100; % Maximum Reactive Power Limit
G = real(ybus); % Conductance
B = imag(ybus); % Susceptance
Psp=(PG-PL)/100;
Qsp=(QG-QL)/100;
Ssp= Psp +1i*Qsp;
PV= find(type==2);
PQ= find(type==3);
NS= [PV;PQ];
M=length(PV); %no. of PV bus
N=length(PQ);%no. of PQ bus
O=length(NS);
dP=zeros(O,1);
dQ=zeros(N,1);
MV = [dP; dQ];
iter=1;
tol=1;
st=clock; % stop the iteration time clock
while (tol>1e-8)
S = V.*conj(ybus*V);
S1 = S-Ssp;
Qc = imag(S);
dP=real(S1(NS));
s=0; %set pv violation
for i = 1:bus
if type(i) == 2
if (Qc(i) > Qmax(i)) || (Qc(i) < Qmin(i))
if Qc(i) < Qmin(i) % Whether violated the lower limit.
Qsp(i) = Qmin(i);
else % No, violated the upper limit.
Qsp(i) = Qmax(i);
end
type(i) = 3; % If Violated, change PV bus to PQ bus.
end
s=s+1;
end
end
PV= find(type==2);
PQ= find(type==3);
NS= [PV;PQ];
M=length(PV); %no. of PV bus
N=length(PQ);%no. of PQ bus
O=length(NS);
dQ=imag(S1(PQ));
B1= -imag(ybus(NS,NS));
B11= -imag(ybus(PQ,PQ));
invB1 = inv(B1);
invB11 = inv(B11);
dth=-invB1*dP; %correction vector : theta
dV=-invB11*dQ; %correction vector : Voltage
MV=[dP; dQ];
ANG(NS) = ANG(NS)+ dth;
Vmag(PQ) = Vmag(PQ)+ dV;
V=Vmag.*exp(1i*ANG);
iter = iter+1;
tol =max(abs(MV));
end
PVnew= find(type==2);
PQnew= find(type==3);
dell = ANG*180/pi;
for m=1:bus
for n=1:bus
I(m,n)=-(V(m)-V(n))*ybus(m,n);
end
I(m,m)=ybus(m,:)*V;
end
for m=1:bus
for n=1:bus
Pow_Flow(m,n)=V(m)*I(m,n)';
end
Pow_Flow(m,m)=V(m)*I(m,m)';
end
s1=find(type==1);
Slack_bus_power=Pow_Flow(s1,s1)*BaseMVA;
for m=1:bus
for n=1:bus
Loss(m,n)=Pow_Flow(m,n)+Pow_Flow(n,m);
end
Loss(m,m)=0;
end
ste=clock; % stop the iteration time clock
lf = [type Vmag dell real(S).*100 Qc.*100];
Losses = real(Loss);
Total_loop_loss=sum((sum(real(Loss))))/2*BaseMVA;
disp('--------------------------------------------------------------------------')
disp(' FDLF Load-Flow Study');
disp('--------------------------------------------------------------------------')
disp(date);
fprintf('Number of iterations : %d \n', iter);
fprintf('Error tolerance considered : %d \n', tol);
fprintf('Solution time : %g sec.\n\n',etime(ste,st));
fprintf(' Ybus Matrix:\n');
disp('==================');
disp(ybus);
disp(" B'");
disp('========');
disp(B1)
disp(' B"');
disp('========');
disp(B11);
fprintf(' BUS Data:\n');
disp('============');
disp(' bus type v delta Pg Qg Pl Q1 Qmin Qmax');
disp(' ----- ----- --- ----- ---- ---- ----- --- ------ ------');
disp(busdata);
fprintf('Load Flow:\n');
disp('=============');
disp(' Bus Type Vpu angle P(MW) Q(MVAr)');
disp(' -------- -------- -------- -------- -------');
disp(lf);
fprintf('Power loss:%d MW\n \n',Total_loop_loss);

0 comments on commit b373370

Please sign in to comment.