nrlfppg.m

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

% % Program for Newton-Raphson Load Flow Analysis
% % 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 Number..
% type = busdata(:,2);        % Type of Bus 1-Slack, 2-PV, 3-PQ..
% V = busdata(:,3);           % Specified Voltage..
% del = 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..
% nbus = max(bus);            % To get no. of buses..
% P = Pg - Pl;                % Pi = PGi - PLi..
% Q = Qg - Ql;                % Qi = QGi - QLi..
% P = P/baseMVA;              % Converting 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..
% 
% pv = find(type == 2 | type == 1); % Index of PV Buses..
% pq = find(type == 3); % Index of PQ Buses..
% 
% npv = length(pv); % Number of PV buses..
% npq = length(pq); % Number of PQ buses..
% 
% while (Tol > 0.00001)   % 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
%     
%     % Calculate change from specified value
%     dPa = Psp-P;
%     dQa = Qsp-Q;
%     k = 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,nbus-npv);
%     for i = 1:(nbus-1)
%         m = i+1;
%         for k = 1:npv
%             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(nbus-npv,nbus-1);
%     for i = 1:(nbus-npv)
%         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(nbus-npv,nbus-npv);
%     for i = 1:(nbus-npv)
%         m = pq(i);
%         for k = 1:(nbus-npv)
%             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
%     
%     X = inv(J)*M;           % Correction Vector
%     dTh = X(1:nbus-1);
%     dV = X(nbus:length(X));
%     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
%     Iter = Iter + 1;
%     Tol = max(abs(M));
% end
%     
% Iter % Number of Iterations took..
% V;
% del = 180/pi*del;
% [V del] % Bus Voltages and angles..


% Program for Bus Power Injections, Line & Power flows (p.u)...
% clc

Y = ybusppg();          % Calling Ybus program..
nbus = length(Y); 
% Program for Newton-Raphson Load Flow Analysis..
% Assumption, Bus 1 is considered as Slack bus..

            % Calling ybusppg.m to get Bus Admittance Matrix..
busdata = busdata14();      % Calling busdata30.m to get busdatas..
baseMVA = 100;              % Base MVA..
bus = busdata(:,1);         % Bus Number..
type = busdata(:,2);        % Type of Bus 1-Slack, 2-PV, 3-PQ..
V = busdata(:,3);           % Specified Voltage..
del = busdata(:,4);         % Voltage Angle..
Pg = busdata(:,5);          % PGi..
Qg = busdata(:,6);          % QGi..
Pl = busdata(:,7);          % PLi..
Ql = busdata(:,8);
Qin=busdata(:,11);
% QLi..
Qmin = busdata(:,9);        % Minimum Reactive Power Limit..
Qmax = busdata(:,10);       % Maximum Reactive Power Limit..
nbus = max(bus);            % To get no. of buses..
P = Pg - Pl;
%random number
% 
% for i=1:30
%     Q(i)=i;
% end
% A=rand(5,5);
% Q;
% A=A.*10;
% A;
% A(1)=Q(1);
%----------------------------------------------
% 
%    Nh=randperm(5);
% Capsize=[17.5 20 18 16 15];
% for i=1:5
%     Ns(i)=Capsize(Nh(i));
% end
% select=Ns(1);
% M=randperm(14);
% for i=1:14
%     Qin(i)=0;
% end
% for i=1:5
%     N(i)=M(i);
% end
% for i=1:14
%     for j=1:5
%     if i==N(j)
%   Nh=randperm(5);
% sel=Nh(1);
%     Qin(i)=Ns(3);
%        else
%         j=j+1;        
%     end
%     end
% end
% Qin=sparse(Qin)
%-----------------------------------
%for i=1:2
% Pi = PGi - PLi..

Q = Qg - Ql+Qin;                % Qi = QGi - QLi..
P = P/baseMVA;              % Converting 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..

pv = find(type == 2 | type == 1); % Index of PV Buses..
pq = find(type == 3); % Index of PQ Buses..

npv = length(pv); % Number of PV buses..
npq = length(pq); % Number of PQ buses..

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
                if Q(n) < Qmin(n)
                    V(n) = V(n) + 0.01;
                elseif Q(n) > 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
    
    X = inv(J)*M;           % Correction Vector
    dTh = X(1:nbus-1);
    dV = X(nbus:end);
    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
    Iter = Iter + 1;
    Tol = max(abs(M));
end
    
Iter; % Number of Iterations took..
V;  % Bus Voltage Magnitudes in p.u ...
Del = 180/pi*del; % Bus Voltage Angles in Degree...
Vm = pol2rect(V,del);
% disp('------------------------------');
% disp('|  Bus  |    V    |  Angle   | ');
% disp('|  No   |   pu    |  Degree  | ');
% disp('------------------------------');
% for m = 1:nbus
%     fprintf('%4g', m); fprintf('    %8.4f', V(m)); fprintf('    %8.4f', Del(m)); fprintf('\n');
% end
% disp('-----------------------------');

%for
   % Converting polar to rectangular..
Iij = zeros(nbus,nbus);
% Iji = zeros(nbus,nbus);
Sij = zeros(nbus,nbus);
Sji = zeros(nbus,nbus);
Lij = zeros(nbus,nbus);


% Line Current Flows..
 for m = 1:nbus
     for n = 1:nbus
         if m ~= n
             Iij(m,n) = -(Vm(m) - Vm(n))*Y(m,n);%Vm(m)*b(m);
             Iji(n,m)=-(Vm(n)-Vm(m))*Y(n,m); % Y(m,n) = -y(m,n)..
         end
     end
 end
 Iij = sparse(Iij);
 Iji=sparse(Iji);
 Iijm = abs(Iij);
 Iija = angle(Iij);
sum=0;t=0;
 % Load Flows..
 for m = 1:nbus
     for n = 1:nbus
         if m ~= n
             Sij(m,n) = Vm(m)*conj(Iij(m,n));
             Tij(m,n)=100*real(Sij(m,n));
             sum=sum+Tij(m,n);
             Mij(m,n)=100*imag(Sij(m,n));
             t=t+Mij(m,n);
%              Sji(n,m)=Vm(n)*conj(Iji(n,m));
%              Lij(m,n)=(Sij(m,n)+Sji(n,m));
        end
     end
 end

  for m = 1:nbus
     for n = 1:nbus
         if m ~= n
             Sji(m,n) = Vm(m)*conj(Iji(m,n));
             
             
         
%              Sji(n,m)=Vm(n)*conj(Iji(n,m));
%              Lij(m,n)=(Sij(m,n)+Sji(n,m));
        end
     end
  end
 Sji=sparse(Sji);
 Sij = sparse(Sij);
Tij=sparse(Tij);
Mij=sparse(Mij);

 Pij = sparse(real(Sij));
 Qij = sparse(imag(Sij));
%  Lij=sparse(Lij);
%  realloss=100*sparse(real(Lij));
%  imajinaryloss=100*sparse(imag(Lij));
%  tloss=0;
%  for i=1:length(realloss)
%      tloss=tloss+realloss(i);
%  end
%  tloss
%  realloss;
%  imajinaryloss;

 % Bus Power Injections..
 Si = zeros(nbus,1);
 for i = 1:nbus
     for k = 1:nbus
         Si(i) = Si(i) + conj(Vm(i))* Vm(k)*Y(i,k);
     end
 end
 Pi = 100*real(Si);
 Qi = 100*(-imag(Si));


%M=100*real(S);
%N=100*imag(S);
Vm
real_powerloss=t

imag_powerloss=sum
%end

% n=0;
% for u=1:14
%     n=n+Ql(u,1)
% end