Newton Raphson - Fluxo Potência

Newton Raphson - Fluxo Potência

% -ROTINA COMPUTACIONAL DESENVOLVIDA EM ATENDIMENTO A DISCIPLINA: SEP I.
% -SOLUÇÃO DE UM SISTEMA COM 4 BARRAS, PELO MÉTODO DE NEWTON RAPHSON COMPLETO.
%ALUNO: Lucas Alexander Mapa - 130900052
%ALUNO: Rafael Teixeira Gomes - 130900039
%------------------------------------------------------------------------------------
%LEITURA DOS DADOS DO SISTEMA
clc; clear all; close all; 
%Bloco responsável por aquisição dos dados da planilha excel
matriz_dados = xlsread('dados.xlsx');
aguarde  = waitbar(0, 'Carregando dados do sistema...');
waitbar(1,aguarde); close(aguarde);
v = xlsread('dados.xlsx',2);
teta = xlsread('dados.xlsx',3);
Pesp = xlsread('dados.xlsx',4);
Qesp = xlsread('dados.xlsx',5);
aguarde  = waitbar(0, 'Carregando dados do sistema...');
waitbar(1,aguarde); close(aguarde);
tipodebarra = [2 0 0 3];% Tipo de barra: 0-PQ; 3-SWING; 2-PV.
dimensao = length(tipodebarra);
% Artifício computacional para possibilitar cálculos subsequentes
for cont = 1:dimensao
    V(cont,1) = v(1,1);
    TETA(cont,1) = teta(1,1);
    for cont2 = 2:dimensao
        V(cont,cont2) = v(1,cont2);
        TETA(cont,cont2) = teta(1,cont2);
        
    end
     
end
%------------------------------------------------------------------------------------
% MONTAGEM DA MATRIZ DE ADMITÂNCIAS YBUS:
%Define a dimensão da matriz Ybus e  multiplica a parte real por (i)
[NL,NC] = size(matriz_dados);
Z_linha = matriz_dados(:,3) + (i*matriz_dados(:,4));
%Elementos fora da diagonal da Ybus
for k = 1:NL
    Ybij = -1/Z_linha(k);
    i = matriz_dados(k,1);
    j = matriz_dados(k,2);
    Ybus(i,j)= Ybij; 
    Ybus(j,i)= Ybij;
    
end
%Elementos da diagonal da Ybus
for i = 1:NC
    for j=1:NL
        if i== matriz_dados(j,1)
            Ybus(i,i) = Ybus(i,i) + (1/Z_linha(j));
        end
         if i== matriz_dados(j,2)
            Ybus(i,i) = Ybus(i,i) + (1/Z_linha(j));
         end
    end
       
end
% Fim da montagem da matriz de admitâncias Ybus
%------------------------------------------------------------------------------------
%VARIAVEIS AUXILIARES
 B=imag(Ybus);
 G=real(Ybus);
matrizP1 = ones(dimensao,1);
matrizQ1 = ones(dimensao,1);
matrizP2 = zeros(dimensao,1);
matrizQ2 = zeros(dimensao,1);
matrizP3 = zeros(dimensao,1);
matrizQ3 = zeros(dimensao,1);
 rodar = 1;
 tol=1e-3;
 itermax=5;
 iteracoes = 0;
%INICIO DO PROCESSO ITERATIVO
  while (rodar ==1 & iteracoes<itermax)
   
       
  % Loop para o calculo das potências ativas e reativas da diagonal da matriz
 for k = 1:(dimensao)
   
     Pcalc(k) =  G(k,k)*(V(1,k)^2); 
     Qcalc(k) = -B(k,k)*(V(1,k)^2);
 end    
    % Loop para o calculo das potências ativas e reativas forada diagonal
    % mais a potência da diagonal da matriz
        for a=1:(dimensao)
            for b = 1:(dimensao)
                 if a~=b       
    Pcalc(a) = Pcalc(a)+((V(a,a))*(V(a,b)))*(G(a,b)*cos(TETA(1,b)-TETA(1,a))- B(a,b)*sin(TETA(1,b)-TETA(1,a)));
    Qcalc(a) = Qcalc(a)-(((V(a,a))*(V(a,b)))*(G(a,b)*sin(TETA(1,b)-TETA(1,a))+B(a,b)*cos(TETA(1,b)-TETA(1,a))));
                end
            end  
        end
         
   % Variaveis auxiliares para o cáculo dos deltas P's e Q's do sistema
for k = 1:(dimensao)
    Piteracao(k) = (matrizP1(k)+ matrizP2(k)*V(k)+matrizP3(k)*(V(k)^2))*Pesp(1,k);
    Qiteracao(k) = (matrizQ1(k)+ matrizQ2(k)*V(k)+matrizQ3(k)*(V(k)^2))*Qesp(1,k);
   
       
    %Cáculo de deltas P's e Q's
    
        if tipodebarra(k) == 0 
        deltaP(k) = Piteracao(k) - Pcalc(k);
        deltaQ(k) = Qiteracao(k) - Qcalc(k);
       
        elseif tipodebarra(k) == 2
        deltaP(k) = Piteracao(k) - Pcalc(k);
       
        elseif tipodebarra(k) == 3
        deltaP(k) = 0;
        deltaQ(k) = 0;
                    
    end
end
%Montagem do vetor delta D
transp_deltaP = deltaP';
transp_deltaQ = deltaQ';
deltaD = [transp_deltaP;transp_deltaQ];
%---------------------------------------------------------------------------------
%Condição para cálculo da Jacobiana
contador = 0;
 for a = 1:dimensao
     if ((abs(deltaP(1,a))<tol) & (abs(deltaQ(1,a))<tol )) 
contador = contador +1;
     end
 end
 %---------------------------------------------------------------------------------
 if contador ~= dimensao
    %MONTAGEM DA MATRIZ JACOBIANA
        
        %Jacobiana =|   H    |    N    |
        %           |   J    |    L    |
        
       
         %Cálculo das diagonais da matriz jacobiana
        for k =1:dimensao
            H(k,k) = -Qcalc(k)-((V(1,k)^2)*B(k,k));
             if tipodebarra(k)==3
                 H(k,k) = 10^15; %Artifício do valor numerico grande (VNG) 
             end
            N(k,k) = (Pcalc(k)+(V(1,k)^2)*G(k,k)); 
            J(k,k) = Pcalc(k)-(V(1,k)^2)*G(k,k);
            L(k,k) = (Qcalc(k)-(V(1,k)^2)*B(k,k)); 
            if tipodebarra(k)>= 2 %Válido para barra de tensão e swing. 
                L(k,k)=10^15; % (VNG)
            end
        end
        
        % Calculo dos elementos fora da diagonal da matriz Jacobiana
       for a = 1:dimensao
           for b = 1:dimensao
           if a ~=b
               H(a,b) = (-V(1,a)*V(1,b))*(B(a,b)*cos(TETA(1,b)-TETA(1,a))+G(a,b)*sin(TETA(1,b)-TETA(1,a)));
               N(a,b) = (V(1,a)*V(1,b))*(G(a,b)*cos(TETA(1,b)-TETA(1,a))-B(a,b)*sin(TETA(1,b)-TETA(1,a)));
               J(a,b) = (V(1,a)*V(1,b))*(B(a,b)*sin(TETA(1,b)-TETA(1,a))-G(a,b)*cos(TETA(1,b)-TETA(1,a)));
               L(a,b) = (-V(1,a)*V(1,b))*(B(a,b)*cos(TETA(1,b)-TETA(1,a))+G(a,b)*sin(TETA(1,b)-TETA(1,a)));
           end
               
           end
       end
      JACOB = [H,N;J,L]; % Matriz jacobiana completa
 else
      rodar =0;
break
 end
       %Fim da montagem da matriz jacobiana
       
       %-----------------------------------------------------------------------------------------
       %MÉTODO LDU PARA CALCULO DA MATRIZ JACOBIANA INVERSA
       
 n=size(JACOB);
for i=1:n
    for j=1:n
        if i<j
            l(i,j)=0;
        end
        if i>j
            U(i,j)=0;
        end
    end
end
for i=1:n
    l(i,i)=1;
    for j=1:n
        U(j,j)=1;
        if i>j
            l(i,j)=JACOB(i,j);
        end
        if i<j
            U(i,j)=JACOB(i,j);
        end
        D(1,1)=JACOB(1,1);
        if i>j
            for k=j:-1:2
                if k==j
                    l(i,j)=JACOB(i,j);
                end
                l(i,j)=l(i,j)-(l(i,k-1)*D(k-1,k-1)*U(k-1,j));
            end
        end
        if i<j
            for k=i:-1:2
                if k==i
                    U(i,j)=JACOB(i,j);
                end
                U(i,j)=U(i,j)-(l(i,k-1)*D(k-1,k-1)*U(k-1,j));
            end
        end
        if i==j
            for k=i:-1:2
                if k==i
                    D(i,j)=JACOB(i,j);
                end
                D(i,j)=D(i,j)-(l(i,k-1)*D(k-1,k-1)*U(k-1,j));
            end
        end
        if i>j
            l(i,j)=l(i,j)/D(j,j);
        end
        if i<j
            U(i,j)=U(i,j)/D(i,i);
        end
    end
end 
  
%Fim da aquisição das matrizes L D e U para calculo implicito da inversa da matriz Jacobiana
%-------------------------------------------------------------------------------------
JACOB_INV = inv(U)*inv(D)*inv(l);
%Calculo do vetor de tetas e ângulos deltaX
deltaX = (JACOB_INV*deltaD);
deltaX(4,1) = deltaX(4,1)*(-1);
%Separa o vetor deltaX em dois outros vetores auxiliares compostos apenas por tetas e por ângulos
teta_deltaX = deltaX(1:dimensao,  1:1);
tensoes_deltaX = deltaX((dimensao+1):2*dimensao);
 % vetor das novas tensões e vetor dos novos tetas
NOVAS_TENSOES = (tensoes_deltaX + (V(1:1, 1:dimensao))')'; 
NOVOS_TETAS = (teta_deltaX +(TETA(1:1, 1:dimensao)'))'; 
%Artifício computacional para possíbilitar a rotina computacional
for cont4 = 1:dimensao
    V(cont4,1) = NOVAS_TENSOES(1,1);
    TETA(cont4,1) = NOVOS_TETAS(1,1);
    
    for cont3 = 2:dimensao
        V(cont4,cont3) = NOVAS_TENSOES(1,cont3);
        TETA(cont4,cont3) = NOVOS_TETAS(1,cont3);
       
    end
     
end
% incremento no número total de iterações para possibilitar a finalização
% do loop
  iteracoes = iteracoes +1;
   
  end
     if iteracoes < 5
  %	Resultados Obtidos: Tensões, ângulos e potências
fprintf('\nNúmero de iterações necessárias para convergência:\n');
iteracoes
fprintf('\nTensão nas barras (pu):\n');
V = V(1:1, 1:dimensao )
fprintf('\nÂngulos: \n');
TETA = TETA(1:1, 1:dimensao )
fprintf('\nPotência ativa calculada (pu)\n');
Pcalc
fprintf('\nPotência reativa calculada (pu)\n');
Qcalc
  
     else
       fprintf('\O método de Newton Raphson atingiu seu limite máximo permitido de iterações\n');
       fprintf('\nO sistema não converge\n');
     end
     
  % FIM DA ROTINA COMPUTACIONAL
     

Comentários