% Program for Newton-Raphson Load flow analsis
function [Sl V] = NewtonRaphson(Np,Nd, Nt, xMin, xMax, vMin, vMax, R, k, Vel, Bus, p, V)
clear;
clc
nbus=30;
Y = ybusppg(nbus);
busd = busdatas(nbus);
BMva = 100;
bus = busd(:,1); % Bus number
type = busd(:,2); % Type of Bus 1-Slack, 2-PV, 3-PQ
V = busd(:, 3); % Specified Voltage
del = busd(:,4); % Voltage angle
%Ubusdt = UpdatedBusdatas(Np,Nd, Nt, xMin, xMax, vMin, vMax, R,k,Vel, Bus,p,V);
Pg = busd(:,5)/BMva; % PGi...
Qg = busd(:,6)/BMva; % QGi...
Pl = busd(:,7)/BMva; % PLi...
Ql = busd(:,8)/BMva; % QLi..
Qmin = busd(:,9)/BMva; % Minimum Reactive Power Limit
Qmax = busd(:,10)/BMva; % Maximum Reactive Power Limit
P = Pg-Pl; % Pi = PGi -PLi
Q = Qg-Ql; % Qi = QGi-QLi
Psp = P; % P Specified
Qsp = Q; % Q Specified
G = real(Y); % Condactance Matrix...
B = imag(Y); % Susceptance Matrix..
pv = find(type == 2 |type == 1); % PV Buses..
pq = find(type == 3); % PQ Buses..
npv = length(pv); % No. of PV Buses..
npq = length(pq); % No. of PQ BUses..
Tol = 1;
Iter = 1;
while (Tol > 1e-5) % Iteration Starting...
P = zeros(nbus,1);
Q = zeros(nbus,1);
end
% 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 violations..
for n= 2:nbus
if type(n) == 2
QG = Q(n)+Ql(n);
if QG < Qmin(n)
V(n) = V(n)+0.01;
elseif QG > Qmax(n)
V(n) = V(n) -0.01;
end
end
end
% Calculate 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]; % Mismatch Vectors
% Jacobian
% J1 - Derivatives of Real Power injectons 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- Derivatives 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- Derivatives 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- Derivatives 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]; % Jacobian Matrix...
X = inv(J)*M; % Correction vector
dTh = X(1:nbus-1); % Change in Voltage angle
dV = X(nbus:end); % Change in Voltage Magnitude
% Updating state State Vectors..
del(2:nbus) = dTh + del(2:nbus); % Voltage Angle
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)); % Tolerance
PLQL = loadflow(nbus,V,del,BMva); % Calling Loadflow
V;
SL = PLQL;
end
% Proogram to for Admitance and Impedance Bus formation
function Y = ybusppg(num) % Return Y
linedata = linedatas(num); % Calling Linedatas
fb = linedata(:,1); % From Bus Number
tb = linedata(:,2); % To bus number
r = linedata(:,3); % Resistance, R
x = linedata(:,4); % Reactance, X
b = linedata(:,5); % Ground Admittance, B/2
a = linedata(:,6); % Tap Setting value
z = r + i*x; % z matrix
y = 1./z; % To get iinverse of each element
b = i*b; % Make B imaginary
nb = max(max(fb),max(tb)); % No. of buses
nl = length(fb); % No. of branches
Y = zeros(nb,nb); % Initialise Ybus
% Formulation of the off Diagnal Elements
for k = 1:nl
Y(fb(k),tb(k)) = Y(fb(k),tb(k)) - y(k)/a(k);
Y(tb(k),fb(k)) = Y(fb(k),tb(k));
end
% Formulation of Diagonal Elements
for m = 1:nb
for n = 1:nl
if fb(n) == m
Y(m,m) = Y(m,m) + y(n)/(a(n)^2) + b(n);
elseif tb(n) == m
Y(m,m) = Y(m,m) + y(n) + b(n);
end
end
end
end
评论0