Relatorio de análise

Prévia do material em texto

FACULDADE INDEPENDENTE DO NORDESTE – FAINOR
CURSO DE BACHARELADO EM ENGENHARIA ELÉTRICA
DIEGO SANTOS PRADO
EDIMAR MARES FRANÇA
GRAZIELA BRAZ PITANGUEIRAS
SIMULAÇÃO DO SISTEMA TESTE DE 13 BARRAS IEEE NO MATLAB UTILIZANDO O MÉTODO DE NEWTON RAPHSON
VITÓRIA DA CONQUISTA – BA
2018
DIEGO SANTOS PRADO
EDIMAR MARES FRANÇA
GRAZIELA BRAZ PITANGUEIRAS
SIMULAÇÃO DO SISTEMA TESTE DE 13 BARRAS IEEE UTILIZANDO O MÉTODO DE NEWTON RAPHSON
Relatório apresentado a disciplina Análise de Sistemas Elétricos, do curso de Bacharelado em Engenharia Elétrica, da Faculdade Independente do Nordeste, como requisito para a obtenção de nota da III unidade.
Professora: Selma Alves de Oliveira.
VITÓRIA DA CONQUISTA – BA
2018
RESUMO
Os estudos para a obtenção de resultados quando se trata de Sistemas Elétricos de Potência (SEP) inclui diversos métodos que são baseados em simulações do fluxo de potência das redes tanto de distribuição quanto de transmissão, estão utilizando-se ferramentas computacionais e métodos de resoluções bem confiáveis para se obter informações fundamentais do sistema. 
O trabalho em questão propõe uma técnica composicional através do software Matlab, que encontra dados do fluxo de carga em um sistema real de 13 barras adquirido no Institute of Electrical and Electronics Engineers (IEEE). O sistema foi escolhido devido possuir todos os dados necessários para se aplicar o método de Newton Raphson. Este método é baseado no processo de convergência do sistema em questão, para se obter certos tipos de dados que são importantes para o entendimento do sistema como níveis de tesão. Com a aplicação desse método é possível identificar níveis indesejados de tensão e outros tipos de eventuais distúrbios assim ajustando onde for necessário. 
Palavras-chave: Carga. fluxo. IEEE. método. Newton.
LISTA DE FIGURAS
Figura 1 - Sistema de duas barras no modelo pi para análise do fluxo de carga	18
Figura 2 - Curva da função g com relação ao vetor x	20
Figura 3 - Interpretação do método de Newton Raphson	23
Figura 4 - Alimentador de teste de nó IEEE 13	24
Figura 5 - Diagrama unifilar do sistema IEEE de 13 barras	24
Figura 6 - Esquemático do alimentador de teste de nó IEEE 13	27
Figura 7 - Diagrama de linha única para rede de distribuição IEEE 13-barras	29
Figura 8 - Matriz admitância gerada no MATLAB	81
Figura 9 - Matriz admitância gerada no MATLAB	81
Figura 10 - Matriz admitância gerada no MATLAB	82
Figura 11 - Matriz admitância gerada no MATLAB	82
Figura 12 - Ângulos e tensões geradas da iteração no MATLAB	83
 
LISTA DE TABELAS
Tabela 1 - Parâmetros de equipamentos nas expressões de fluxo	19
Tabela 2 - Valores das tensões e ângulos calculados da iteração	33
Tabela 3 - Valores das tensões e ângulos calculados da iteração	37
Tabela 4 - Valores das tensões e ângulos calculados da iteração	41
Tabela 5 - Valores das tensões e ângulos calculados da iteração	44
Tabela 6 - Valores das tensões e ângulos calculados da iteração	48
Tabela 7 - Valores das tensões e ângulos calculados da iteração	52
Tabela 8 - Valores das tensões e ângulos calculados da iteração	56
Tabela 9 - Parâmetros de equipamentos nas expressões de fluxo	60
Tabela 10 - Valores das tensões e ângulos calculados da iteração	63
Tabela 11 - Valores das tensões e ângulos calculados da iteração	67
Tabela 12 - Valores das tensões e ângulos calculados da iteração	71
Tabela 13 - Valores das tensões e ângulos calculados da iteração	75
Tabela 14 - Valores das tensões e ângulos calculados da iteração	80
INTRODUÇÃO
O processo para a análise de sistemas elétricos tem se evoluído com o passar dos anos e com isso, existe sempre a necessidade de melhoria tais análises. Para que a distribuição da energia elétrica seja feita de forma satisfatória, é de suma importância que as tensões de entrega à carga esteja operando em condições de regime permanente, isto é, não pode haver muitas oscilações no sinal da tensão. Sendo de bastante relevância o estudo do fluxo de carga observando seu comportamento, planejamento e operação no sistema de energia elétrica, onde este imprime a análise em função do tempo. 
O estudo do fluxo de carga é realizado para simular o sistema operando em regime permanente, isto é, o sistema operando sem a interferência dos transitórios eletromagnéticos causados por efeitos externos a rede. Já á análise tem por necessidade encontrar três parâmetros comuns que são: tensão (em módulo e fase), potência (ativa, reativa e aparente) e a corrente. Existem métodos numéricos para essa em um SEP, podendo estes serem calculados a mão, porém quando se tem muitas barras no sistema fica impossível, em alguns métodos, se calcularem a mão, então, sendo o projetista da rede fazer esta análise através de softwares que realizam esta análise. 
Este relatório apresenta uma simulação de um sistema elétrico de 13 barras do IEEE, a qual é utilizada para encontrar as tensões nas barras através do método de Newton Raphson, observando suas potencias ativas e reativas, com a utilização o software MATLAB para a implementação do mesmo. Sendo de extrema importância os conhecimentos obtidos em sala de aula para a realização aprofundada dos estudos sobre a análise de fluxo de carga e acerca do método abordado, já que o objetivo central do projeto é observar o comportamento do fluxo durante toda a extensão do SEP. 
Dentre os vários métodos de análise de sistemas elétricos de potência, neste trabalho trataremos do método de Newton Raphson que entre as diversas funcionalidades tanto dentro do cálculo numérico e em outras aplicações matemáticas é também utilizado para o propósito de analise sistemas. Sabendo disso dentro do mesmo assunto esse método é bastante utilizado em sistemas que a quantidade de barras são maiores ou seja, esse método corresponde bem quando se está tratando sistemas grandes.
MÉTODO DE NEWTON RAPHSON
Quando se trata de SEP vêm logo em mente sistemas de enorme porte e de grande complexidade, sendo necessário a criação de métodos, os quais são essenciais para a compreensão do sistema tanto quando se está analisando uma parte do sistema quanto ao todo. O método de Newton Raphson entra com um diferencial nesse aspecto, já que é um método iterativo o qual aproxima um conjunto de equações não-lineares simultâneas por um conjunto de equações lineares usando expansão por séries de Taylor. Tem como vantagem a robustez já que converge quase sempre e com menos interações que outros métodos assim diminuindo bastante em ralação a memória e desenvolvimento de processo, além disso a convergência do sistema independe da dimensão do sistema. 
Ao decorrer do trabalho ira ser abordado tanto as técnicas do método que torna possível e satisfatória em relação aos resultados obtidos, quanto a programação que foi realizada no software Matlab, visto que se trata de uma técnica que é bastante utilizada em sistema de grande porte, então envolve sistemas que se calculado a mão, são bem trabalhosos e propício a erros humanos (VIENA, 2013). Em seguida é apresentado o problema prático, onde é demonstrado o sistema elétrico e seus parâmetros e também o cálculo da matriz de admitância feita a mão, para se analisar os resultados simulados no softwares. Posteriormente, é mostrado o processo de simulação do código do sistema de 13 barras no software e enfim, os resultados finais e dificuldades do projeto e conclusão.
O sistema que irá ser utilizado para se aplicar o método de Newton Raphson é um sistema de 13 Barras, o qual é um sistema real que foi adquirido no IEEE13, porem é importante salientar que um sistema com essa quantidade de barras pode ser definido tanto sendo um conjunto de subestações onde cada barra pode ser uma cidade quanto também podemos definir esses sistemas sendo um conjunto de transformadores em uma cidade, portanto pode-se notarque o método utilizado é bastante útil. 
ANÁLISE DE FLUXO DE CARGA EM SEP
Análise de fluxo de carga, ou conhecido também como fluxo de potência, é um estudo realizado em sistemas elétricos de potência onde consiste essencialmente na determinação das tensões complexas nas barras, das distribuições dos fluxos de carga que fluem pelas linhas e dentre outros diversas demandas. A modelagem matemática para análise de fluxo de carga baseia-se pelo modelo estático onde utiliza-se equações algébricas para solução de alguns parâmetros.
A análise de fluxo de carga possuem infinitas aplicações, uma delas são: permitir a determinação do estado operativo do sistema elétrico, verificar se o sistema em análise está operando adequadamente e indicar, indiretamente, o que deve ser feito para corrigir ou prevenir situações inadequadas de operação e análise de curto-circuito nas condições de pré-falta. Eles são basicamente utilizados tanto para o planejamento da rede como também para a operação. No planejamento, segundo Moura (2013), a análise de fluxo de carga permite a comparação de alternativas de expansão, bem como a avaliação do impacto no sistema da entrada de novas unidades geradoras, como é o caso da geração distribuída. Já na operação, segundo Moura (2013) o estudo de fluxo de potência visa, com informações mais confiáveis a respeito da carga, definir o melhor perfil de tensões para a operação do sistema bem como os ajustes de tapes dos transformadores, condições para chaveamento de banco de capacitores, etc.
Nas barras de referência, que comumente são as barras de geração e chamadas de barra 1, são as barras onde os valores de potência aparente ( e as tensões já são conhecidas e as outras barras do sistemas podem ser do tipo PQ (potência aparente da barra conhecida), barra do tipo PV (potência ativa e módulo da tensão conhecidas) e barra do tipo Vθ (onde o módulo e a fase da tensão são conhecidas). Na geração é necessário que se tenham os valores da potência ativa (PG), da potência reativa (QG) e do módulo da tensão controlada (V). Já nas cargas é necessário conhecer as potências que cada barra, onde estão dispostas as cargas, consomem, em tais barras é o ponto onde o objetivo da aplicação do fluxo de carga é realizado. Para a análise de fluxo de carga, devemos tomar como base um sistema elétrico simples a seguir:
Figura 1 - Sistema de duas barras no modelo pi para análise do fluxo de carga
Fonte: BORGES (2005).
O objetivo é encontrar o fluxo de potência ativa e reativa entre as barras k e m e calcular as perdas ativas e reativas da linha. Analisando o fluxo da corrente saindo da barra k para a barra m teremos a seguinte equação:
Já para o fluxo da corrente saindo da barra m para a barra k, a equação segue:
Para diferentes equipamentos, utilizamos a seguinte tabela de parâmetros nas expressões de fluxo:
Tabela 1 - Parâmetros de equipamentos nas expressões de fluxo
	
	
	
	
	Linha de transmissão
	1
	0
	
	Transformador em fase
	
	0
	0
	Transformador defasador puro
	1
	
	0
	Transformador defasador
	
	
	0
Fonte: OLIVEIRA (2017).
Temos que,
A partir da expressão final da dedução acima do fluxo da potência aparente de k para m, temos que as potências ativa e reativa da barra k para a barra m são:
De maneira análoga pode-se obter os fluxos de potência ativa e reativa saindo da barra m em direção à barra k, assim:
Observe que as expressões de Pmk e Qmk foram obtidas simplesmente trocando os índices k e m nos parâmetros. Logo, teremos então, que as perda ativa é dada pela seguinte expressão:
E a perda reativa é obtida pela seguinte expressão:
MÉTODOS DE SOLUÇÃO DE FLUXO DE CARGA
Anteriormente foram vistas expressões genéricas para o cálculo de fluxo de carga em um sistema contendo 2 barras. Porém, para os cálculos para fluxo de carga em sistemas elétricos onde se encontram muitas barras só podem ser realizadas usando os métodos matemáticos conhecidos como método de Newton Raphson e de Gauss-Seidel, existem outros métodos também, mas este dois são os principais (ARAÚJO, 2005). 
Os métodos de solução para análise de fluxo de carga possuem dois métodos interativos, um baseado na matriz Y (admitância), outro na matriz Z (impedância) ou o método de Newton. Os métodos baseados na matriz Y podem ser baseados no sistema da equação linear como por exemplo o método de Gauss e Gauss-Seidel e o método de Newton é onde pode ser utilizado tanto o método jacobiano como o método desacoplado para encontrar os resultados das tensões nas barras. Daremos enfoque neste capítulo sobre o Método de Newton Raphson (ARAÚJO, 2005). 
Método de solução através de Newton Raphson
Considere um sistema n-direcional.
Pretende-se determinar o valor de x onde a função não-linear g(x) se anula. Em termos geométricos corresponde ao ponto xs da figura a seguir:
Figura 2 - Curva da função g com relação ao vetor x
Fonte: OLIVEIRA (2017).
Onde temos sendo uma função vetorial (nx1) e x o vetor das incógnitas (nx1), ou seja, (1.1)
A linearização da função vetorial para é dada pelos dois primeiros termos da série de Taylor:
 (1.2)
Sendo J a matriz jacobiana dada por:
 (1.3)
O vetor de correção é calculado impondo-se que:
 (1.4)
Que é a maneira linearizada de se resolver o problema . No caso particular em que, por exemplo, n=2, a equação (1.4) assume a forma:
O algoritmo para resolução do sistema de equações:
Pelo método de Newton Raphson, temos a algoritmo genérico abaixo:
Fazer e escolher uma solução inicial 
Calcular 
Testar Convergência: para i = 1, 2, ...,, n o processo convergiu para caso contrário, realizar o próximo passo;
Calcular a matriz jacobiana ;
Determinar a nova solução:
Fazer e retornar ao passo II.
Fluxo de Potência pelo método de Newton Raphson
Agora serão feitas a aplicação do método de Newton Raphson para resolução do subsistema 1, o ponto central do processo de resolução consiste em determinar o vetor de correção o que exige a resolução do seguinte sistema linear:
Como e são constantes, o Jacobiano é simplificado para:
Assim, o sistema passa a ser:
As componentes das submatrizes jacobianas H, N, M e L são dadas por:
O método de Newton-Raphson aplicado ao problema do fluxo de potência possui então o seguinte algoritmo:
Fazer v = 0 e escolher os valores iniciais dos ângulos das tensões das barras PQ e PV, e os módulos das tensões das barras PQ(V = V0);
Calcular para as barras PQ e PV, e para as barras PQ e determinar os resíduos e;
Testar a convergência: se e , o processo interativo convergiu para ; caso contrário, desvia para a próxima etapa;
Calcular a matriz Jacobiana:
Determinar a nova solução: sendo e determinadas resolvendo o sistema linear:
Fazer e retornar ao passo II.
Temos a seguir a curva que interpreta o método de solução através de Newton-Raphson.
Figura 3 - Interpretação do método de Newton Raphson
Fonte: OLIVEIRA (2017).
SISTEMA TESTE IEEE DE 13 BARRAS
O IEEE disponibiliza sistemas testes de variadas quantidade de barras, neste relatório é abordado o sistema teste de 13 barras, e demostrado algumas das possibilidade de configurações desse sistema, demostrando que algumas configurações podem ser alteradas para fins de testes e estudos. A Figura 4 apresenta o sistema fornecido pelo IEEE, sistema base de onde outros podem se originar.
Figura 4 - Alimentador de teste de nó IEEE 13
Fonte: IEEE (2000-2010).
Esta rede permite testar várias conexões de transformadores, na forma de rebaixador, reguladores de tensão e unidades de geração distribuída com as cargas balanceadas ou desbalanceadas.
Figura 5 - Diagrama unifilar do sistema IEEE de 13 barras
Fonte: CARVALHO (2012).
Notam-se mudanças em relação ao diagrama apresentado no arquivo original representadona Figura 4, como a inserção de uma unidade de GD na barra 680, a representação dos bancos de capacitores e o ponto de regulação do regulador destacado na barra 671 (Figura 5).
 
 ALIMENTADOR DE TESTE DE NÓ IEEE 13
A seguir apresenta-se uma esquemático de uma simulação de uma adaptação do sistema base do IEEE para 13 barras (Figura 6), com a utilização de dados fornecido pelo IEEE (2000-2010). Este exemplo mostra o uso da ferramenta Load Flow do Powergui para inicializar o circuito do alimentador de teste de nó IEEE 13, utilizando a nova versão teste 2017.b do Matlab.
i) Da descrição do Teste:
O esquemático apresenta doze blocos de fluxo de carga usados para calcular um fluxo desequilibrado, no modelo representa o circuito do alimentador de teste de nó IEEE 13, originalmente publicado pelo relatório do Subcomitê de Análise do Sistema de Distribuição IEEE. Observa-se que o modelo não inclui o transformador de regulação entre os nós 650 e 632 do modelo de teste de referência.
Da Simulação:
Para sua realização são necessário alguns passos descritos pelo Mathworks (1994-2017), os quais serão reescritos abaixo:
Abra a ferramenta de fluxo de carga do Powergui e pressione o botão Compute. 
Pressione o botão 'Aplicar' para aplicar a solução de fluxo ao modelo para iniciar a simulação no estado estacionário. Observa-se que se pode visualizar a magnitude da voltagem do barramento individual e os valores de ângulo de fase na guia "Fluxo de carga" do bloco de fluxo de carga correspondente ou pode especificá-los como anotações de bloco por conveniência.
Pressione o botão "Relatório" para obter um relatório que mostra um resumo de fluxo de carga e resultados de fluxo de carga detalhados em cada ônibus. Os resultados do relatório de fluxo de carga parcial são reproduzidos no subsistema de Resultado do Fluxo de Energia no modelo (o bloco verde da Figura 6). 
Comece a simulação e verifique se ela começa no estado estacionário, com fluxo de energia esperado.
Figura 6 - Esquemático do alimentador de teste de nó IEEE 13
Fonte: MATHWORKS (1994-2017).
DESENVOLVIMENTO
Para o desenvolvimento deste projeto foram realizados alguns testes, a princípio com valores do sistema teste, em outros momentos com valores aleatórios, e por vim com auxílio dos dados fornecidos pelo artigo de Aref et al. (2012). Para a criação do sistema de 13 IEEE, como um todo, foi feita a simulação por passo a passo utilizando a interação entre barras e depois para a confirmação do processo foi executada a simulação com todo o sistema.
	A seguir apresenta o código desenvolvido no software Matlab, foram inseridos os dados mencionados acima, como inicialmente determinando o número de barras que nesta simulação é de 13, foi definido qual é a barra de referência, em seguida, os valores de tensão e fase, que assumimos como 1 e 0, respectivamente. Posteriormente, definiu-se os valores para a potência ativa e reativa gerada na barra de referência, e potência ativa e reativa na carga e definiu-se os valores de capacitores como nulos.
	Definiu-se também as interações a partir das barras e os valores de resistência, reatâncias, capacitores em serie que são nulos, e o ângulo de fase. O próximo passo foi realizar a determinação da dimensão das matrizes, a obtenção dos dados da barra, a obtenção dos trechos entre as barras. Posteriormente, a matriz de admitância, calcula-se as potências da carga e a que fluem entre as barras.
	Em seguida, é calculada a jacobiana, calcula-se as diagonais, cálculo das componentes, H, N, M e L, estruturando assim a matriz jacobiana. Após isso, definiu-se a condição de limite de interações para evitar o loop infinito. Por fim, o cálculo das potências e fluxo de potências. O sistema da Figura 7, apresenta o diagrama de 13 barras monofásico (unifilar).
Figura 7 - Diagrama de linha única para rede de distribuição IEEE 13-barras.
Fonte: AREF et al. (2012) (adaptada pelo autores).
 Implementação e Uso
Implementando o algoritmo em MATLAB em passo a passo com as interações entre as barras:
3.1.1 ITERAÇÃO ENTRE AS BARRAS 1 E 2
 clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 1 e 2. 
 
Y1 = 0.463+0.290i; 
% Convertendo para a matriz de admitãncia A.
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia 
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.386]; 
PotRG = [0.370; 0.380]; 
PotAD = [0; 0];
PotRD = [0; 0];
P = (PotAG - PotAD);
Q = (PotRG - PotRD);
%Estabelecendo as tensões iniciais V1 e V2
V = [1;1];
Vang = [0:0];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
Pc = zeros(Qtbarras,1);
Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
variacaoP=P-Pc;
variacaoQ=Q-Qc;
vP=variacaoP(2:Qtbarras);
x=1;
vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema 
for i=2:2
 	vQ(x,1)=variacaoQ(i);
 	x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % MudandoPara o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
Tabela 2: Valores das tensões e ângulos calculados da iteração.
	Tensão V1
	Ângulo θ
	Tensão V2
	Ângulo θ
	Inicial: 1.0000 pu
	Inicial: 0
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: -1.1251
	Final: 3.18
Fonte: AUTORIA PRÓPRIA (2018).
3.1.2 ITERAÇÃO ENTRE AS BARRAS 2 E 3
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 2 e 3. 
Y1 = 0.240+0.630i;
 
% Convertendo para a matriz de admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia
disp('MATRIZ ADMITÂNCIA')
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 2 e 3 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366];
PotRG = [0.370; 0.342];
PotAD = [0; 0];
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
%Estabelecendo as tensões iniciais da iteração
V = [1;-1.1251];
 
Vang = [0;3.18];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema 
for i=2:2
vQ(x,1)=variacaoQ(i);
 	x=x+1;
end
Mvariacao = [vP;vQ]; 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
Tabela 3: Valores das tensões e ângulos calculados da iteração.
	Tensão V2
	Ângulo θ
	Tensão V3
	Ângulo θ
	Inicial: ;-1.1251pu
	3.18
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: 1.2378 pu
	Final: -112.9771
Fonte: AUTORIA PRÓPRIA (2018).
3.1.3 ITERAÇÃO ENTRE AS BARRAS 3 E 4
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
 
%Impedância entre o barramento 3 e 4. 
Y1 = 0.212+0.440i;
 
% Convertendo para a matriz de admitãncia A.
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22];
disp(A)
 
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
g=real(A);b=imag(A);
 
% Parametros Iniciais das potências do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.386]; 
PotRG = [0.370; 0.380]; 
PotAD = [0; 0]; 
PotRD = [0; 0]; 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
%Estabelecendo as tensão iniciais da iteração 
V = [1;1.2378];
Vang = [0;-112.97];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema 
for i=2:2
vQ(x,1)=variacaoQ(i);
 	x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)-Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
Tabela 4: Valores das tensões e ângulos calculados da iteração.
	Tensão V3
	Ângulo θ
	Tensão V4
	Ângulo θ
	Inicial: 1.2378 pu
	-112.9771
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: 1.0815 pu
	Final: 0.1698
Fonte: AUTORIA PRÓPRIA (2018).
3.1.4 ITERAÇÃO ENTRE AS BARRAS 4 E 5
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 4 e 5. 
 Y1 = 0.312+0.441i;
 
% Convertendo para a matriz de admitãncia A.
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia
disp('MATRIZ ADMITÂNCIA')
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 4 e 5 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366];
PotRG = [0.370; 0.342];
PotAD = [0; 0];
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
%Estabelecendo as tensões iniciais da iteração
V = [1;1.0815];
%Estabelecendo os ângulos inicias da iteração
Vang = [0;0.1698];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
vQ(x,1)=variacaoQ(i);
 	x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 5: Valores das tensões e ângulos calculados da iteração.
	Tensão V4
	Ângulo θ
	Tensão V5
	Ângulo θ
	Inicial: ; 1.0815 pu
	0.1698
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: -1.0152pu
	Final: 0.0743
Fonte: AUTORIA PRÓPRIA (2018).
3.1.5 ITERAÇÃO ENTRE AS BARRAS 5 E 6
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 5 e 6. 
Y1 = 0.380+0.170i;
 
% Convertendo para a matriz de admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia 
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 5 e 6 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366]; 
PotRG = [0.370; 0.342]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
%Estabelecendo as tensões inicias 
V = [1;-1.0152];
%Estabelecendo os ângulos iniciais
Vang = [0;0.0743];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1]; 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:QtbarrasPc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
vQ(x,1)=variacaoQ(i);
 	x=x+1;
end
 Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
Tabela 6: Valores das tensões e ângulos calculados da iteração.
	Tensão V5
	Ângulo θ
	Tensão V6
	Ângulo θ
	Inicial: - 1.0152 pu
	0.0743
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: 1.0788 pu
	Final: 1.5021
Fonte: AUTORIA PRÓPRIA (2018).
3.1.6 ITERAÇÃO ENTRE AS BARRAS 5 E 7
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 5 e 7. 
Y1 = 0.360+0.181i; 
% Convertendo para a matriz de admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia
disp('MATRIZ ADMITÂNCIA')
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 5 e 7 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366]; 
PotRG = [0.370; 0.342]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
 
V = [1;-1.0152];
Vang = [0;0.0743];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
 vQ(x,1)=variacaoQ(i);
 x=x+1;
end
 Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 7: Valores das tensões e ângulos calculados da iteração.
	Tensão V5
	Ângulo θ
	Tensão V7
	Ângulo θ
	Inicial: - 1.0152 pu
	0.0743
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: 0.8164 pu
	Final: 5.8408
Fonte: AUTORIA PRÓPRIA (2018).
3.1.7 ITERAÇÃO ENTRE AS BARRAS 7 E 8
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 7 e 8. 
Y1 = 0.355+0.143i;
 
% Convertendo para a matrizde admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia 
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 7 e 8 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366]; 
PotRG = [0.370; 0.342]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD);
Q = (PotRG - PotRD);
 
V = [1;0.8164]; 
Vang = [0;5.8408];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
 vQ(x,1)=variacaoQ(i);
 x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 8: Valores das tensões e ângulos calculados da iteração.
	Tensão V7
	Ângulo θ
	Tensão V8
	Ângulo θ
	Inicial: 0.8164 pu
	5.8408
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: 2.8922 pu
	Final: -3.7554
Fonte: AUTORIA PRÓPRIA (2018).
3.1.8 ITERAÇÃO ENTRE AS BARRAS 8 E 9
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 8 e 9. 
Y1 = 0.450+0.190i;
 
% Convertendo para a matriz de admitãncia A.
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia 
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22]; 
disp(A) 
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 1 e 2 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366]; 
PotRG = [0.370; 0.342]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
 
V = [1;2.8922]; 
Vang = [0;-3.7554];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
 vQ(x,1)=variacaoQ(i);
 x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))-b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 9: Valores das tensões e ângulos calculados da iteração.
	Tensão V8
	Ângulo θ
	Tensão V9
	Ângulo θ
	Inicial: 2.8922 pu
	-3.7554
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: -1.2623 pu
	Final: 3.0489
Fonte: AUTORIA PRÓPRIA (2018).
3.1.9 ITERAÇÃO ENTRE AS BARRAS 8 E 10
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 2 e 3. 
Y1 = 0.3915+0.235i;
 
% Convertendo para a matriz de admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia 
disp('MATRIZ ADMITÂNCIA')
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 1 e 2 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366];
PotRG = [0.370; 0.342]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
 
V = [1;2.8922]; 
Vang = [0;-3.7554];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
vQ(x,1)=variacaoQ(i);
 	x=x+1;
end
Mvariacao = [vP;vQ]; 
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 10: Valores das tensões e ângulos calculados da iteração.
	Tensão V8
	Ângulo θ
	Tensão V10
	Ângulo θ
	Inicial: 2.8922 pu
	-3.7554
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: -2.1323 pu
	Final: -2.1858
Fonte: AUTORIA PRÓPRIA (2018).
3.1.10 ITERAÇÃO ENTRE AS BARRAS 7 E 11
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 7 e 11. 
Y1 = 0.222+0.150i;
 
% Convertendo para a matriz de admitãncia A.
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia
disp('MATRIZ ADMITÂNCIA')
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências nos barramentos 7 e 11 do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.366]; 
PotRG = [0.370; 0.342]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
 
V = [1;0.8164]; 
Vang = [0;5.8404];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
 vQ(x,1)=variacaoQ(i);
 x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componenteN
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 11: Valores das tensões e ângulos calculados da iteração.
	Tensão V7
	Ângulo θ
	Tensão V11
	Ângulo θ
	Inicial: 0.8164 pu
	5.8408
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: -1.1731 pu
	Final: -1.2536
Fonte: AUTORIA PRÓPRIA (2018).
3.1.11 ITERAÇÃO ENTRE AS BARRAS 11 E 12
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 11 e 12. 
Y1 = 0.1515+0.160i;
 
% Convertendo para a matriz de admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia 
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.386]; 
PotRG = [0.370; 0.380]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
 
V = [1;-1.1731]; 
Vang = [0;-1.2536];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema (entre 2 barramentos) 
for i=2:2
 vQ(x,1)=variacaoQ(i);
 x=x+1;
end
Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
Tabela 12: Valores das tensões e ângulos calculados da iteração.
	Tensão V11
	Ângulo θ
	Tensão V12
	Ângulo θ
	Inicial: -1.1731 pu
	-1.2536
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: 1.0696 pu
	Final: -1.9836
Fonte: AUTORIA PRÓPRIA (2018).
3.1.12 ITERAÇÃO ENTRE AS BARRAS 7 E 13
clc
clear all
 
% Método de Newton Raphson 13 BARRAS
%Iterações feitas entre barras 
%Impedância entre o barramento 7 e 13. 
Y1 = 0.410+0.270i;
 
% Convertendo para a matriz de admitãncia A. 
A12 = -1/Y1; 
A21 = -1/Y1;
A11 = (-A12);
A22 = (-A12);
 
%Matriz Admitancia
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A21 A22];
disp(A)
disp('-----------------------------------------------------')
%Determinação da equação Y=g+jb para a injeção de potência
 g=real(A);b=imag(A);
 
% Parametros Iniciais das potências do sistema
%Dados das potencias ativa e reativa: Gerador e carga 
PotAG = [0.890; 0.386]; 
PotRG = [0.370; 0.380]; 
PotAD = [0; 0]; 
PotRD = [0; 0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
 
V = [1;2.8922]; 
Vang = [0;5.8404];
%Estabelecendo o número de barras
difervt=1;
int=0;
Qtbarras = 2;
Qtbarras_PQ = 1;
Qtbarras_PV = 0;
barras_PQ = [1];
 
% Metodo de Newton-Raphson
%Estabelecendo a taxa de erro e o número de iteração para evitar o loop
%infinito: 8000 (limite de iterações)
%Método de convergência quadrática
while (difervt>0.0001 && int< 8000)
 
 Pc = zeros(Qtbarras,1);Qc = zeros(Qtbarras,1);
%Estabelecimento das potências nas barras P(km) e Q(kVAr) nas barras
%Considerando o ponto inicial a expansão de 1 ordem de Taylor ? x(i)=-(g(x0))/(g^' (x0))
for i=1:Qtbarras
for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
end
end
%Variações do fluxo entre as potências
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
%Determinando a característica da iteração do sistema 
for i=2:2
 vQ(x,1)=variacaoQ(i);
 x=x+1;
end
Mvariacao = [vP;vQ]; 
 
% Formação da Matriz Jacobiana, diagonais, componentes [H N; J L]
% J(x) = [(NPQ+NPV) (NPQ)]
% Formação da componente H
 H=zeros(Qtbarras-1,Qtbarras-1);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras-1;
 n=x+1;
if m==n
for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);
else
 H(i,x)=V(m)*V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente N
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (N)
 N=zeros(Qtbarras-1,Qtbarras_PQ);
for i=1:Qtbarras-1
 m=i+1;
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 N(i,x)=N(i,x)+V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 N(i,x)=N(i,x)+V(m)*g(m,m);
else
 N(i,x)=V(m)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente J
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (J)
 J=zeros(Qtbarras_PQ,Qtbarras-1);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras-1
 n=x+1;
if m==n
for n=1:Qtbarras
 J(i,x)=J(i,x)+V(m)*V(n)*(g(m,n)*cos(Vang(m)- Vang(n))+b(m,n)*sin(Vang(m)-Vang(n)));
end
 J(i,x)=J(i,x)-V(m)^2*g(m,m);
else
 J(i,x)=V(m)*V(n)*(-g(m,n)*cos(Vang(m)-Vang(n))- b(m,n)*sin(Vang(m)-Vang(n)));
end
end
end
 
% Formação da componente L
%Detrminação do número de barras P e Q através das devivadas parciais
%Jacobianas (L)
 L=zeros(Qtbarras_PQ,Qtbarras_PQ);
for i=1:Qtbarras_PQ
 m=barras_PQ(i);
for x=1:Qtbarras_PQ
 n=barras_PQ(x);
if m==n
for n=1:Qtbarras
 L(i,x)=L(i,x)+V(n)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
 L(i,x)=L(i,x)-V(m)*b(m,m);
else
 L(i,x)=V(m)*(g(m,n)*sin(Vang(m)-Vang(n))- b(m,n)*cos(Vang(m)-Vang(n)));
end
end
end
%Equação Jacobiana para cálculo dos valoes de tensão e ângulo nas iterações
 Jacob = [H N;J L] % Declaração da matriz Jacobiana
 X=inv(Jacob);
 L=X*Mvariacao;
 
 vTh=L(1:Qtbarras-1); % Mudando Para o angulo
 vV=L(Qtbarras:end); % Mudando para o modulo
 
 Vang(2:Qtbarras)=Vang(2:Qtbarras)+(vTh) %Atualização do angulo da tensao
 
 p=1;
for n=2:2
 V(n)=V(n)+vV(p); % Atualização da Tensao 
 p=p+1;
end
 variacao=max(abs(Mvariacao));
 int=int+1
end
disp('TENSÕES CALCULADAS')
disp(V)
disp('-----------------------------------------------------') 
 
Tabela 13: Valores das tensões e ângulos calculados da iteração.
	Tensão V7
	Ângulo θ
	Tensão V13
	Ângulo θ
	Inicial: 0.8164 pu
	5.8404
	Inicial: 1.0000 pu
	Inicial: 0
	
	
	Final: -1.1604 pu
	Final: -2.7004
Fonte: AUTORIA PRÓPRIA (2018).
3.1.13 PROGRAMAÇÃO COM TODO O SISTEMA (13 BARRAMENTOS)
Implementando o algoritmo em MATLAB como um todo entre as barras:
clc
clear all
 
% Método de Newton Raphson 13 BARRAS(TODO SISTEMA)
% Impedâncias Fornecidas 
Y1 = 0.463 + 0.290i;
Y2 = 0.240 + 0.630i;
Y3 = 0.212 + 0.440i;
Y4 = 0.312 + 0.441i;
Y5 = 0.380 + 0.170i;
Y6 = 0.360 + 0.181i;
Y7 = 0.355 + 0.143i;
Y8 = 0.450 + 0.190i;
Y9 = 0.312 + 0.185i;
Y10 = 0.222 + 0.150i;
Y11 = 0.1515 + 0.160i;
Y12 = 0.410 + 0.270i;
 
% Convertendo para admitância 
A12 = -1/Y1;
A23 = -1/Y2;
A34 = -1/Y3;
A45 = -1/Y4;
A56 = -1/Y5;
A57 = -1/Y6;
A78 = -1/Y7;
A89 = -1/Y8;
A810 = -1/Y9;
A711 = -1/Y10;
A1112 = -1/Y11;
A713 = -1/Y12;
 
% Definindo ida e volta 
A21 = A12;
A32 = A23;
A43 = A34;
A54 = A45;
A65 = A56;
A75 = A57;
A87 = A78;
A98 = A89;
A108 = A810;
A117 = A711;
A1211 = A1112;
A137 = A713;
 
% definindo quem é zero 
A13=0; A24=0; A31=0; A41=0; A51=0; A61=0; A71=0; A81=0; A91=0; A101=0; A111=0; A121=0; A131=0;
A14=0; A25=0; A35=0; A42=0; A52=0; A62=0; A72=0; A82=0; A92=0; A102=0; A112=0; A122=0; A132=0;
A15=0; A26=0; A36=0; A46=0; A53=0; A63=0; A73=0; A83=0; A93=0; A103=0; A113=0; A123=0; A133=0;
A16=0; A27=0; A37=0; A47=0; A58=0; A64=0; A74=0; A84=0; A94=0; A104=0; A114=0; A124=0; A134=0;
A17=0; A28=0; A38=0; A48=0; A59=0; A67=0; A76=0; A85=0; A95=0; A105=0; A115=0; A125=0; A135=0;
A18=0; A29=0; A39=0; A49=0; A510=0; A68=0; A79=0; A86=0; A96=0; A106=0; A116=0; A126=0; A136=0; 
A19=0; A210=0; A310=0; A410=0; A511=0; A69=0; A710=0; A811=0; A97=0; A107=0; A118=0; A127=0; A138=0;
A110=0; A211=0; A311=0; A411=0; A512=0; A610=0; A712=0; A812=0; A910=0; A109=0; A119=0; A128=0; A139=0;
A111=0; A212=0; A312=0; A412=0; A513=0; A611=0; A813=0; A911=0; A1011=0; A1110=0; A129=0; A1310=0;
A112=0; A213=0; A313=0; A413=0; A612=0; A912=0; A1012=0; A1113=0; A1210=0; A1311=0;
A113=0; A613=0; A912=0; A913=0; A1013=0; A1213=0; A1312=0;
 
% Diagonal principal 
A11 = -(A12);
A22 = -(A21 + A23);
A33 = -(A32 + A34);
A44 = -(A43 + 45);
A55 = -(A54 + A56 + A57);
A66 = -(A65);
A77 = -(A75 + A78 + A711 + A713);
A88 = -(A87 + A89 + A810);
A99 = -(A98);
A1010 = -(A108);
A1111 = -(A117 + A1112);
A1212 = -(A1211);
A1313 = -(A137);
 
%Matriz Admimitancia 
disp('MATRIZ ADMITÂNCIA') 
A = [A11 A12 A13 A14 A15 A16 A17 A18 A19 A110 A111 A112 A113;
 A21 A22 A23 A24 A25 A26 A27 A28 A29 A210 A211 A212 A213;
 A31 A32 A33 A34 A35 A36 A37 A38 A39 A310 A311 A312 A313;
 A41 A42 A43 A44 A45 A46 A47 A48 A49 A410 A411 A412 A413;
 A51 A52 A53 A54 A55 A56 A57 A58 A59 A510 A511 A512 A513;
 A61 A62 A63 A64 A65 A66 A67 A68 A69 A610 A611 A612 A613;
 A71 A72 A73 A74 A75 A76 A77 A78 A79 A710 A711 A712 A713;
 A81 A82 A83 A84 A85 A86 A87 A88 A89 A810 A811 A812 A813;
 A91 A92 A93 A94 A95 A96 A97 A98 A99 A910 A911 A912 A913;
 A101 A102 A103 A104 A105 A106 A107 A108 A109 A1010 A1011 A1012 A1013;
 A111 A112 A113 A114 A115 A116 A117 A118 A119 A1110 A1111 A1112 A1113;
 A121 A122 A123 A124 A125 A126 A127 A128 A129 A1210 A1211 A1212 A1213;
 A131 A132 A133 A134 A135 A136 A137 A138 A139 A1310 A1311 A1312 A1313];
 
 disp(A)
 disp('-----------------------------------------------------')
 g=real(A);b=imag(A);
 % Parametros Iniciais
PotAG = [0.890; 0.386;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0]; 
PotRG = [0.370; 0.380;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0;0.0]; 
PotAD = [0;0;0;0;0;0;0;0;0;0;0;0;0]; 
PotRD = [0;0;0;0;0;0;0;0;0;0;0;0;0];
 
P = (PotAG - PotAD); 
Q = (PotRG - PotRD);
V = [1;1.1;1;1;1;1;1;1;1;1;1;1;1]; 
Vang = [0;0;0;0;0;0;0;0;0;0;0;0;0];
 
difervt=1;
int=0;
Qtbarras = 13;
Qtbarras_PQ = 12;
Qtbarras_PV = 0;
barras_PQ = [2;3;4;5;6;7;8;9;10;11;12;13];
 
% Metodo de Newton-Raphson
while (difervt>0.0001 && int< 10000)
 Pc = zeros(Qtbarras,1);
 Qc = zeros(Qtbarras,1);
 for i=1:Qtbarras
 for x=1:Qtbarras
 Pc(i)=Pc(i)+V(i)*V(x)*(g(i,x)*cos(Vang(i)- Vang(x))+b(i,x)*sin(Vang(i)-Vang(x)));
 Qc(i)=Qc(i)+V(i)*V(x)*(g(i,x)*sin(Vang(i)-Vang(x))- b(i,x)*cos(Vang(i)-Vang(x)));
 end
 end
 
 variacaoP=P-Pc;
 variacaoQ=Q-Qc;
 vP=variacaoP(2:Qtbarras);
 x=1;
 vQ=zeros(Qtbarras_PQ,1);
 for i=2:13
 vQ(x,1)=variacaoQ(i);
 x=x+1;
 end
 Mvariacao = [vP;vQ];
 
% Formação da Matriz Jacobiana [H N; J L]
 % Formação de H
 H=zeros(Qtbarras-1,Qtbarras-1);
 for i=1:Qtbarras-1
 m=i+1;
 for x=1:Qtbarras-1;
 n=x+1;
 if m==n
 for n=1:Qtbarras
 H(i,x)=H(i,x)+V(m)*V(n)*(-g(m,n)*sin(Vang(m)- Vang(n))+b(m,n)*cos(Vang(m)-Vang(n)));
 end
 H(i,x)=H(i,x)-V(m)^2*b(m,m);