nrmethod.rar_NR power flow

% Program for Newton-Raphson Load Flow Analysis.. nbus =33; busd= [ 1 0 1 0 0.32 0.1425 2 3 1 0 0 0 3 3 1 0 0.23 0.1425 4 3 1 0 0.23 0.1425 5 3 1 0 0 0 6 3 1 0 0 0 7 3 1 0 0.23 0.1425 8 3 1 0 0.32 0.1425 9 3 1 0 0 0 10 3 1 0 0.23 0.1425 11 3 1 0 0.137 0.084 12 3 1 0 0.072 0.045 13 3 1 0 0.072 0.045 14 3 1 0 0.072 0.045 15 3 1 0 0.0135 0.0075 16 3 1 0 0.23 0.1425 17 3 1 0 0.23 0.1425 18 3 1 0 0.23 0.1425 19 3 1 0 0.23 0.1425 20 3 1 0 0.23 0.1425 21 3 1 0 0.23 0.1425 22 3 1 0 0.23 0.1425 23 3 1 0 0.23 0.1425 24 3 1 0 0.23 0.1425 25 3 1 0 0.23 0.1425 26 3 1 0 0.137 0.085 27 3 1 0 0.075 0.048 28 3 1 0 0.075 0.048 29 3 1 0 0.075 0.048 30 3 1 0 0.057 0.0345 31 3 1 0 0.057 0.0345 32 3 1 0 0.057 0.0345 33 3 1 0 0.057 0.0345 ]; 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.. P = busd(:,5)/BMva; % PGi.. Q = 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; for i=10 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 % 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.. e=inv(J2-J1*inv(J3)*J4); dvP=inv(J2-J1*inv(J3)*J4)*dP; dvq=inv(J4-J3*inv(J1)*J2)*dQ; Vnew1=V(2:nbus)+dvP; Vnew2=V(2:nbus)+dvq; end % results goes worong what mistake i had made sir tell me