nrlfppg2.m

This is inportent program for energy sector engineers. This matlab code is based on Newton Rapshion

% Program for Newton-Raphson Load Flow Analysis, alternate code..
% Assumtions
%   1. Bus number 1 is assumed as slask bus..
%   2. No Tap changing transformers..

Y = ybusppg();  % Calling ybusppg.m to get Bus Admittance Matrix..
busdata = busdata14();  % Calling busdata30.m to get busdatas..
baseMVA = 100;          % Base MVA..
bus = busdata(:,1);         % Bus
type = busdata(:,2);        % Type of Bus 1-Slack, 2-PV, 3-PQ
V = busdata(:,3);           % Specified Voltage.
Th = busdata(:,4);          % Voltage Angle.
Pg = busdata(:,5);          % PGi
Qg = busdata(:,6);          % QGi
Pl = busdata(:,7);          % PLi
Ql = busdata(:,8);          % QLi
Qmin = busdata(:,9);        % Minimum Reactive Power Limit
Qmax = busdata(:,10);       % Maximum Reactive Power Limit
n = max(bus);               % To get no. of buses
P = Pg - Pl;                % Pi = PGi - PLi
Q = Qg - Ql;                % Qi = QGi - QLi
P = P/baseMVA;              % Convert to p.u
Q = Q/baseMVA;
Qmin = Qmin/baseMVA;
Qmax = Qmax/baseMVA;
Tol = 10;     % Tolerence kept at high value.
Iter = 1;  % iteration starting
Psp = P;
Qsp = Q;
G = real(Y);    % Conductance..
B = imag(Y);    % Susceptance..

while (Tol > 0.00001)
    P = zeros(n,1);
    Q = zeros(n,1);
    % Calculate P and Q
    for i = 1:n
        for k = 1:n
            P(i) = P(i) + V(i)*V(k)*abs(Y(i,k))*cos(Th(i)-Th(k)-angle(Y(i,k)));
            Q(i) = Q(i) + V(i)*V(k)*abs(Y(i,k))*sin(Th(i)-Th(k)-angle(Y(i,k)));
        end
    end
    % Calculate change from specified value
    dPa = Psp-P;
    dQa = Qsp-Q;
    k = 1;
    for i = 1:n
        if type(i) == 3
            dQ(k,1) = dQa(i);
            k = k+1;
        end
    end
    dP = dPa(2:n);
    M = [dP; dQ];       % Mismatch Vector
    
    % Jacobian
    for i = 2:n
        for k = 2:n
            ii = i-1;
            kk = k-1;
            if k ~= i
                J1(ii,kk) = V(i) * V(k) * abs(Y(i,k)) * sin(Th(i)-Th(k)-angle(Y(i,k)));
                J2(ii,kk) = V(i) * abs(Y(i,k)) * cos(Th(i)-Th(k)-angle(Y(i,k)));
                J3(ii,kk) = -V(i) * V(k) * abs(Y(i,k)) * cos(Th(i)-Th(k)-angle(Y(i,k)));
                J4(ii,kk) = V(i) * abs(Y(i,k)) * sin(Th(i)-Th(k)-angle(Y(i,k)));
            else
                J1(ii,ii) = -Q(i) - V(i)^2*B(i,i);
                J2(ii,ii) = P(i)/V(i) + V(i)*G(i,i);
                J3(ii,ii) = P(i) - V(i)^2*G(i,i);
                J4(ii,ii) = -Q(i)/V(i) - V(i)*B(i,i);
            end
        end
    end
    % To remove rows & columns corresponding to PV buses from Jacobian.
    r = find(type == 2);
    l = length(r);
    for i = l:-1:1
        J2(:,(r(i)-1)) = [];
        J4(:,(r(i)-1)) = [];
    end
    for i = l:-1:1
        J3((r(i)-1),:) = [];
        J4((r(i)-1),:) = [];
    end
    J = [J1 J2; J3 J4];     % Jacobian
    
    X = inv(J)*M;           % Correction Vector
    dTh = X(1:n-1);
    dV = X(n:length(X));
    Th(2:n) = dTh + Th(2:n);
    k = 1;
    for i = 2:n
        if type(i) == 3
            V(i) = dV(k) + V(i);
            k = k+1;
        end
    end
    Iter = Iter + 1;
    Tol = max(abs(X));    
end
    
Iter % Number of Iterations took..
V;
del = 180/pi*Th;
[V del] % Bus Voltages and angles..