牛顿拉夫法matlab,matlab实现牛顿拉夫逊潮流法

% Program for Newton-Raphson Load Flow Analysis..

nbus = 14; % IEEE-14, IEEE-30, IEEE-57..

Y = ybusppg(nbus); % Calling ybusppg.m to get Y-Bus Matrix..

busd = busdatas(nbus); % Calling busdatas..

BMva = 100; % Base MVA..

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..

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); % Conductance 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);

% 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..

if Iter <= 7 && Iter > 2 % Only checked up to 7th iterations..

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

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 Vector

% Jacobian

% J1 - Derivative of Real Power Injections 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 Injections 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 - Derivative of Reactive Power Injections 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 Reactive Power Injections 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 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..

end

loadflow(nbus,V,del,BMva); % Calling Loadflow.m..

获取完整代码请阅读全文

隐藏内容需要支付：¥10

查看隐藏内容需要支付：

¥10

支付宝

请使用

支付宝 扫码支付

本文地址：https://ai.52learn.online/2050