%  This program is to calculate power flow using Newton-Raphonson iterative
%  method. Developed by Wang shishan, 2010-10-27.

%  Part 1: Creating admittance matrix Y;
     clear;
     Nn=3;           % Number of nodes;
     y=[0               1.25-j*3.75  5-j*15
         1.25-j*3.75        0          1.667-j*5
         5-j*15          1.667-j*5       0];
     Y=[y(1,2)+y(1,3)  -y(1,2)  -y(1,3);
             -y(1,2)     y(1,2)+y(2,3)  -y(2,3);
             -y(1,3)      -y(2,3)    y(1,3)+y(2,3)];
      G=real(Y);      % Real part, conductive matrix;
      B=imag(Y);      % Image part, susceptance matrix;
  

%  Part 2: Initialing voltage, power value;
     SB=100 ;                   % MVA, base power;
     Pin=[0;20;0];              % Input active power at 2th node;
     Pout=[0;50;60];           % Output active power at 2th node;
     Qin=[0;0;0];                % Inpouring reactive power;
     Qout=[0; -20;25];         % Outing reactive power;
     Q2_range=[0 35];         % Range of min and max for Q2;
     U=[1.05;1.03;1.0];
     sita=angle(U);              % Phase angle of voltage;
     Um=abs(U);                % Amplititude of voltage;
     Nmax=10;
     TL=1E-06;                  % Tolerance level;

% Part 3: Iterative process;    
    for Ite=1:Nmax                   % Iterative vary;                   
       % Part 3.1: Solving error vectors of active and reactive powers;
       for L=2:3        % L: No. of node;
         Temp1=0;
         Temp2=0;
         for k=1:3
           delta=sita(L)-sita(k);
           Temp1=Temp1+Um(k)*(G(L,k)*cos(delta)+B(L,k)*sin(delta));  
           Temp2=Temp2+Um(k)*(G(L,k)*sin(delta)-B(L,k)*cos(delta)); 
         end
         P(L)=Um(L)*Temp1;
         Q(L)=Um(L)*Temp2;
       end
       deltaP=(Pin-Pout)/SB-P';
       deltaQ=(Qin-Qout)/SB-Q';
       deltaPQ=[deltaP(2) deltaP(3) deltaQ(3)]';
       if max(deltaPQ)<TL
            disp('Iteration is sucessful.')
            Ite
            Um
            sita
            break;   
        end
      
     
       % Part 3.2: Creating Jacobi elements coreponding matrix J;
       % H matrix;   
         for L=2:3
           H(L,L)=Q(L)+B(L,L)*Um(L)*Um(L);   
           for k=(L+1):3
             delta=sita(L)-sita(k);
             Temp1=G(L,k)*sin(delta)-B(L,k)*cos(delta);
             H(L,k)=-Um(L)*Um(k)*Temp1;
             H(k,L)=H(L,k);
           end
         end
       % N matrix;
         delta=sita(2)-sita(3);
         Temp1=G(2,3)*cos(delta)+B(2,3)*sin(delta);
         N23=-Um(2)*Um(3)*Temp1;
         N33=-P(3)-G(3,3)*Um(3)*Um(3);
      % J matrix;   
        J32=-N23;
        J33=-P(3)+G(3,3)*Um(3)^2;
     % L matrix;
       L33=-Q(3)+B(3,3)*Um(3)^2;
     % Creating Jacobi matrix;   
       Jacobi=[H(2,2)  H(2,3)  N23;
               H(3,2)  H(3,3)  N33;
               J32      J33       L33]
        
     % Part 3.3: Solving the correcting equation;    
       Correct=JacobideltaPQ;
    
     % Part 3.4: Correcting voltage for PQ nodes and angles for PV nodes;     
       for k=2:3
         sita(k)=sita(k)-Correct(k-1);
       end
       Um(3)=Um(3)-Correct(3)*Um(3);
     
     % Part 3.5: Checking Q2?
       Q2_G=SB*Q(2)-Qout(2);
       if Q2_G>Q2_range(2)
          Q(2)=Q2_range(2)/SB;
       elseif Q2_G<Q2_range(1)
          Q(2)=Q2_range(1)/SB;  
       else
          disp('Reactive power at 2th node in a norm range');
       end
    end       % END-iteration.
   
% Part 4: Postprocess----Power flowing in every branch;
    for k=1:Nn
       U(k)=Um(k)*exp(j*sita(k));  
    end 
    for I=1:Nn
    for J=1:Nn
         S(I,J)=U(I)*conj((U(I)-U(J))*y(I,J));   
    end   
    end       

%THE END, 2010-10-28, 2010-11-3, Wang Shishan.