Matlab - flujo de potencia por Gauss Seidel en MATLAB

clc;clear;

ZDATA=[1 2 0.02 0.04;1 3 0.01 0.03;2 3 0.0125 0.0250];

na=ZDATA(:,1);nb=ZDATA(:,2);R=ZDATA(:,3);X=ZDATA(:,4);nbr=length(ZDATA(:,1));

nb1=length(ZDATA(:,1));nbus=max(max(na),max(nb));

Z=R+X*1i;

y=ones(nbr,1)./Z;

Y=zeros(nbus,nbus);

for k=1: nbr

    if na(k)>0 && nb(k)>0

        Y(na(k),nb(k))=Y(na(k),nb(k))-y(k);

        Y(nb(k),na(k))=Y(na(k),nb(k));

    end

end

    for n=1: nbus

        for k=1:nbr

        if na(k)==n || nb(k)==n

            Y(n,n)=Y(n,n)+y(k);

        end

        end

    end

   n1=na;nr=nb;

%inicia el ciclo de interaciones para determinar voltajes nodales

    %Datos el ejemplo 6.12 

    Edata=[1 0 0 1.05 0;2 -4 -2.5 1.0 0;3 2 0 1.04 0];

    ntn=length(Edata(:,1)); P=Edata(:,2); Q=Edata(:,3); mv=Edata(:,4);

    av=Edata(:,5).*pi/180;

    nosc= input('No. de nodo oscilante: ');

    nvc= input('No. de nodo de voltaje controlado: ');

    tol= input('Valor del error: ');

    Pbase=input('potencia base :');

    V= mv.*(cos(av)+ 1i*sin(av));

    iter=0;

    VA=V;

    while iter<50

    for k1=1:ntn

            sum=0;

      if k1 == nosc

            else

                for k2=1:ntn

                    if k2 == k1 % igual a k

                    else

                        %n diferente de k

                        sum=sum+Y(k1,k2)*V(k2);

                    end

                    if k1 == nvc %Se calcula la potencia reactiva

                        suma=0;

                        for k3=1:ntn

                            suma=suma+Y(k1,k3)*V(k3);

                        end

                        qima=conj(V(k1))*suma;

                        Q(k1)= -imag(qima);

                    else

                    end

                end

                %se determina voltaje calculado y se corrije

                d=0;

                while d<=1;

                    V(k1)=((1+0i)/(Y(k1,k1))*((P(k1)-1i*Q(k1))/conj(V(k1))-sum));

                    d=d+1;

                end

      end

    end

        % se calcula el error(partes real e imaginaria:a)

        for k4=1:ntn

        e(k4)=abs(real(V(k4))-real(VA(k4)));

        f(k4)=abs(imag(V(k4))-imag(VA(k4)));

        end

        %Se comparan las partes real e imaginaria con el error

        for k=1:ntn

        if e(k)<tol && f(k)<tol

            flag=1;

        else

            flag=0;

        end

        end

           if flag==1

        %Si la respuesta es si, se llego a la convergencia, se interrumpen

       %las interacciones y se mantiene la magnitud de voltaje constante en

        %el nodo de voltaje controlado

          anvc=angle(V(nvc));

          V(nvc)=Edata(nvc,4)*(cos(anvc) + 1i*sin(anvc));

          break

            else

                %Si la respuesta es no,se continua el ciclo de interaciones

            end

        iter=iter+1;

        %manteniendo la magnitud de voltaje constante en el nodo de voltaje

        %controlado

        anvc=angle(V(nvc));

        V(nvc)=Edata(nvc,4)*(cos(anvc) + 1i*sin(anvc));

        VA=V;

    end

    %Se terminas las iteraciones

    %Se calcula las potencias de los nodos

    for k7=1:ntn

        suma=0;

        for k8=1:ntn

            suma=suma+Y(k7,k8)*V(k8);

        end

        S(k7)=conj(V(k7))*suma*Pbase;

    end

    P=real(S);

    Q=-imag(S);

    % se calculas las corrientes en las lineas de transmision

    for k9=1:nb1

        C(n1(k9),nr(k9))=(V(n1(k9))-V(nr(k9)))/Z(k9);

    end

% Se calcula potencia trasmitida y la potencia recibida en c/linea

    for k11=1:nb1

        ST(n1(k11),nr(k11))=V(n1(k11))*conj(C(n1(k11),nr(k11)))*Pbase;

        SR(n1(k11),nr(k11))=V(n1(k11))*conj(C(n1(k11),nr(k11)))*Pbase;

    end

    %Se hace la preparacion de la informacion para la impresion en pantalla

    PG=zeros(ntn);QG=zeros(ntn);PC=zeros(ntn);QC=zeros(ntn);

    for k12=1:ntn

    if k12==nosc || k12==nvc

        PG(k12)=P(k12);QG(k12)=Q(k12);PC(k12)=0;QC(k12)=0;

    else

        PG(k12)=0;QG(k12)=0;PC(k12)=P(k12);QC(k12)=Q(k12);

    end

    end

        %Se imprime en pantalla el reporte final

        fprintf ('\nNodo  Mag. Voltaje  Ang.fase     Generacion(MW y MVAR)       Cargas(MW y MVAR)\n');

                for k5=1:ntn

         fprintf('%g       %3.3f     %3.3f     %5.2f    %5.2f           %5.2f  %5.2f \n\n\n',k5,abs(V(k5)),angle(V(k5))*180/pi,PG(k5),QG(k5),PC(k5),QC(k5));

                end

                 fprintf('\n Flujos de potencia \n');

                  fprintf('\nNodo - Nodo  P(MW)    Q(MVAR)    P(MW)  Q(MVAR)\n');

                  PT=real(ST);QT=imag(ST);PR=real(SR);QR=imag(SR);

                  for k13=1:nb1

                          fprintf('\n%g      %g    %5.2f  %5.2f   %5.2f    %5.2f\n',n1(k13),nr(k13),PT(n1(k13),nr(k13)),QT(n1(k13),nr(k13)),PR(n1(k13),nr(k13)),QR(n1(k13),nr(k13)));

                  end

                  fprintf('\n No. de interaciones : %g \n',iter);

                  fprintf('\n Dante omar cruz #11580499 \n');

                  fprintf('\n Sitemas Electricos De Potencia \n');

                  clear;

         %Fin del programa