%  This program is to calculate power flow using Gass-Seidel iterative
%  method. Developed by Wang shisan, 2008-11-22 firstly, 2010-10-22 secondaly.

%  Part 1: Creating admittance matrix Y;
   clear;
   y13=1.17-j*4.71;   % Admittances linked 2-nodes;
   y12=5.88-j*23.5;
   y30=j*0.33;
   Y=[y12+y13  -y12  -y13;
       -y12     y12   0;
       -y13      0    y13+y30];

%  Part 2: Initialing voltage, power value;
   U(1)=1;   % Swing node;
   U3=1.1;   % Magnitude for PV node;
   P(3)=0.4;
   S(2)=-0.8-j*0.6;
   Nmax=10;
   Precision=1E-06;   % Tolerance level;
     
%  Part 3: Iterative process;
   U(2)=1.0;  % Setting initial voltage for 2-th node, need to solve.
   Sita3=0;   % Phase angle for 3th node, need to solve.
   MaxError=0;
   U(3)=U3*exp(j*Sita3);
   for k=1:Nmax
     U_old(2:3)=U(2:3);  %Save old voltage in order to determine congergence;
     % For PQ node, 2th node;
        Temp1=conj(S(2)/U(2))-Y(2,1)*U(1)-Y(2,3)*U(3);   %conjugated value;
        U(2)=Temp1/Y(2,2);
     % For PV nodes, 3th node;
        Temp3=conj(U(3))*(Y(3,3)*U(3)+Y(3,1)*U(1)+Y(3,2)*U(2));
        Q(3)=-imag(Temp3);
        S(3)=P(3)+j*Q(3);
        Temp3=conj(S(3)/U(3))-Y(3,1)*U(1)-Y(3,2)*U(2);
        U(3)=Temp3/Y(3,3);
        Sita3=angle(U(3));
        U(3)=U3*exp(j*Sita3);
     % Convergence criterion;
       for p=2:3
          Error(p)=U(p)-U_old(p);
       end
       MaxError=max(abs(Error(2)),abs(Error(3)));
       if MaxError<Precision
          disp('Iteration successful');
          disp('Iteration number=');
          disp(k);
          MaxError
          break
       end
       k
       U
       Sita3 
       MaxError
   end    
 
% Part 4: Output messages;
   if MaxError>Precision
     disp('Calculation is fail.');
     disp('Inceasing Nmax or deceasing presion!')
   end
%THE END, 2010-10-20, Wang Shishan.