Gauss Seidel Matalb

Prévia do material em texto

UNIVERSIDADE FEDERAL DA BAHIA
ESCOLA POLITÉCNICA
DEPARTAMENTO DE ENGENHARIA ELÉTRICA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA ELÉTRICA
ENGA80: ANÁLISE DE SISTEMAS DE ENERGIA ELÉTRICA
1º TRABALHO – EQUIPE 8
SOLUÇÃO DE UM SISTEMA DE QUATRO BARRAS USANDO O METÓDO DE GAUSS SEIDEL IMPLEMENTADO NO MATLAB E MODELAGEM NO POWER WORLD
ALUNOS: 
	 ANDRÉ JUAREZ
 	 DANIEL COELHO VASCONCELOS
 	 HUMBERTO ROCHA
QUESTÃO 1
IMPLEMENTAÇÃO MATLAB:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programa de Pós Graduação em Engenharia Elétrica
% Disciplina Análise de Sistemas de Energia 2016.1
% Método Gauss-Seidel
% Questão 1
% Prof.: Niraldo
% Alunos: 1 - André Duarte
% 2 - Daniel Vasconcelos
% 3 - Humberto Rocha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc;
clear;
 
tic;
nb = 4; %quantidade de barras
nvc = 4; %número da barra de tensão controlada;
e = 1; %tolerância
 
% Bases adotadas
Sb = 100; % 100 MVA
Vb = 230; % 230 kV
 
% Montando a matriz admitância de barra
Yb = zeros(nb,nb);
Yb(1,1) = 8.985190 - 44.835953*i;
Yb(2,1) = -3.815629 + 19.078144*i;
Yb(3,1) = -5.169561 + 25.847809*i;
Yb(1,2) = Yb(2,1);
Yb(2,2) = Yb(1,1);
Yb(4,2) = Yb(3,1);
Yb(1,3) = Yb(4,2);
Yb(3,3) = 8.193267 - 40.863838*i;
Yb(4,3) = -3.023705 + 15.118528*i;
Yb(2,4) = Yb(3,1);
Yb(3,4) = Yb(4,3);
Yb(4,4) = Yb(3,3);
 
% Vetor de tensões iniciais nas barras
Vo = [1; 1; 1; 1.02];
 
% Vetor de Potência Ativa e Reativa
Sc = zeros(nb,1);
Sg = zeros(nb,1);
 
% Barra 1 é a de referência
Sc(1) = (-50-30.99*i)/Sb;
% Barra 2 e 3 de carga
Sc(2) = (-170-105.35*i)/Sb;
Sc(3) = (-200-123.94*i)/Sb;
% Barra 4 é a de tensão controlada
Sc(4) = (-80-49.58*i)/Sb;
Sg(4) = 318/Sb;
 
P = real(Sc + Sg);
Q = imag(Sc + Sg);
 
% Metódo Gauss-Seidel
V = Vo;
k= 0;
Sg1o = zeros(nb,1);
while (e >= 10e-5)
 k = k+1;
 Icalc = zeros(nb,1);
 for m=2:nb
 soma = 0;
 for j=1:nb;
 if (m ~= j) 
 soma = soma + Yb(m,j)*V(j);
 end
 Icalc(m) = Yb(m,j)*V(j)+Icalc(m);
 end 
 if (m == nvc) % Calculo de Q e V para a barra PV
 Qcalc(k) = -imag(conj(V(m))*Icalc(m));
 V(m) = ((P(m)-i*Qcalc(k))/conj(V(m)) - soma)/ Yb(m,m);
 V(m) = abs(Vo(4))*V(m)/abs(V(m)); % Corrigindo o módulo da tensão da barra PV
 else
 V(m) = ((P(m)-i*Q(m))/conj(V(m)) - soma)/ Yb(m,m); 
 end 
 end
 
 % Calculando o valor da corrente e potência gerada na Barra de Referência - Barra 1
 Icalc(1) = 0;
 for j=1:nb;
 Icalc(1) = Yb(1,j)*V(j)+Icalc(1);
 end
 Sg(1) = V(1)*conj(Icalc(1))-Sc(1);
 
 %Potência gerada na Barra 4
 Sg(4) = V(4)*conj(Icalc(4))-Sc(4); 
 
 Sg = [Sg(1); 0; 0; Sg(4)]; %Vetor com potência gerada em cada barra
 er = max(abs(real((Sg - Sg1o)))); % Cálculo do erro entre as Potências Calculadas na iteração
 ei = max(abs(imag((Sg - Sg1o))));
 e = max(er,ei);
 Sg1o = Sg;
end
 
% Cálculo dos Fluxos de Potência nas Linhas
Yshunt = zeros(nb);
Yshunt(1,2) = 0.05125*i;
Yshunt(2,1) = Yshunt(1,2);
Yshunt(1,3) = 0.03875*i;
Yshunt(3,1) = Yshunt(1,3);
Yshunt(2,4) = Yshunt(1,3);
Yshunt(4,2) = Yshunt(2,4);
Yshunt(3,4) = 0.06375*i;
Yshunt(4,3) = Yshunt(3,4);
for m=1:nb
 for j=1:nb
 if (m ~= j)
 SL(j,m) = V(j)*(conj((V(j)-V(m))*(-Yb(j,m)))); % Potencia na linha
 S(j,m) = SL(j,m)- abs(V(j))*(abs(V(j))*Yshunt(j,m)); % Potencia na saida e entrada da barra
 I(j,m) = (S(j,m)/V(j))*Sb*10^3/(sqrt(3)*Vb); % Corrente na saída e entrada da barra
 end
 end
end
 
Perdas = sum((Sg+Sc))*Sb;
SL = SL*Sb; % Potência nas Linhas em MVA
S = S*Sb;
SLp = [SL(1,2)+SL(2,1);SL(1,3)+SL(3,1);SL(2,4)+SL(4,2);SL(3,4)+SL(4,3)]; %Perdas nas linhas (Não considera o shunt)
Sp = [S(1,2)+S(2,1);S(1,3)+S(3,1);S(2,4)+S(4,2);S(3,4)+S(4,3)]; %Perdas na Linhas (Considerando o Shunt)
SLpr = real(SLp); %Perdas Ativas em MW na linha
SLpi = imag(SLp); %Perdas Reativas em MVar na linha sem o Yshunt
Spi = imag(Sp); % Perdas Reativas em MVar considerando o Yshunt
Sg = Sg*Sb; % Potência gerada por barra em MVA
Stg = sum(Sg)*Sb; % Potência total gerada no sistema em MVA
Stc = sum(Sc)*Sb; % Potência total das cargas do sistema em MVA
t1=toc;
 
fprintf(1,' \n'); 
fprintf(1,'**********************************************\n');
fprintf(1,'*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***\n');
fprintf(1,'**********************************************\n');
fprintf(1,' \n'); 
fprintf(1,'Quantidade de iterações: %d \n\n',k); 
fprintf(1,'Tensões Finais nas barras em PU:\n V1 = %f/_%f°\n V2 = %f/_%f°\n V3 = %f/_%f°\n V4 = %f/_%f°\n\n',abs(V(1)),rad2deg(angle(V(1))),abs(V(2)),rad2deg(angle(V(2))),abs(V(3)),rad2deg(angle(V(3))),abs(V(4)),rad2deg(angle(V(4)))); 
fprintf(1,'Correntes entrando e saindo em cada Barra em A:\n I12 = %f/_%f° I21 = %f/_%f°\n I13 = %f/_%f° I31 = %f/_%f°\n I24 = %f/_%f° I42 = %f/_%f°\n I34 = %f/_%f° I43 = %f/_%f°\n\n',abs(I(1,2)),rad2deg(angle(I(1,2))),abs(I(2,1)),rad2deg(angle(I(2,1))),abs(I(1,3)),rad2deg(angle(I(1,3))),abs(I(3,1)),rad2deg(angle(I(3,1))),abs(I(2,4)),rad2deg(angle(I(2,4))),abs(I(4,2)),rad2deg(angle(I(4,2))),abs(I(3,4)),rad2deg(angle(I(3,4))),abs(I(4,3)),rad2deg(angle(I(4,3)))); 
fprintf(1,'Perdas Ativa e Reativa em cada Linha em MW e MVAr:\n PLinha12 = %f QLinha12(sem Yshunt) = %f QLinha12(com Yshunt) = %f\n PLinha13 = %f QLinha13(sem Yshunt) = %f QLinha13(com Yshunt) = %f\n PLinha24 = %f QLinha24(sem Yshunt) = %f QLinha24(com Yshunt) = %f\n PLinha34 = %f QLinha34(sem Yshunt) = %f QLinha34(com Yshunt) = %f\n\n',SLpr(1),SLpi(1),Spi(1),SLpr(2),SLpi(2),Spi(2),SLpr(3),SLpi(3),Spi(3),SLpr(4),SLpi(4),Spi(4));
fprintf(1,'Perda total no Sistema (Linhas) em MW e MVAr:\n Ptotal = %f\n Qtotal(sem Yshunt)= %f Qtotal(com Yshunt)= %f\n\n', sum(SLpr), sum(SLpi), sum(Spi)); 
fprintf(1,'Potência Ativa Gerada nas barras em MW:\n P1 = %f\n P2 = %f\n P3 = %f\n P4 = %f\n\n',real(Sg(1)),real(Sg(2)),real(Sg(3)),real(Sg(4)));
fprintf(1,'Potência Reativa Gerada nas barras em MVAr:\n Q1 = %f\n Q2 = %f\n Q3 = %f\n Q4 = %f\n\n',imag(Sg(1)),imag(Sg(2)),imag(Sg(3)),imag(Sg(4)));
fprintf(1,'Tempo total de execução: %f s\n', t1);
**********************************************
*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***
**********************************************
 
Quantidade de iterações: 15 
Tensões Finais nas barras em PU:
 V1 = 1.000000/_0.000000°
 V2 = 0.982421/_-0.976235°
 V3 = 0.969005/_-1.872250°
 V4 = 1.020000/_1.522949°
Correntes entrando e saindo em cada Barra em A:
 I12 = 112.106359/_29.952375° I21 = 126.614300/_-139.947930°
 I13 = 290.303258/_31.957700° I31 = 300.624852/_-144.913443°
 I24 = 385.769726/_-149.624664° I42 = 376.209021/_27.823683°
 I34 = 309.084272/_-147.730699° I43 = 293.399127/_27.000920°
Perdas Ativa e Reativa em cada Linha em MW e MVAr:
 PLinha12 = 0.226733 QLinha12(sem Yshunt) = 1.133664 QLinha12(com Yshunt) = -8.937735
 PLinha13 = 1.031482 QLinha13(sem Yshunt) = 5.157410 QLinha13(com Yshunt) = -2.356099
 PLinha24 = 1.715483 QLinha24(sem Yshunt) = 8.577417 QLinha24(com Yshunt) = 0.805906
 PLinha34 = 1.835435 QLinha34(sem Yshunt) = 9.177178 QLinha34(com Yshunt) = -3.441307
Perda total no Sistema (Linhas) em MW e MVAr:
 Ptotal = 4.809133
 Qtotal(sem Yshunt)= 24.045669 Qtotal(com Yshunt)= -13.929235
Potência Ativa Gerada nas barras em MW:
 P1 = 186.816033
 P2 = 0.000000
 P3 = 0.000000
 P4 = 317.992645
Potência Reativa Gerada nas barras em MVAr:
 Q1 = 114.499840
 Q2 = 0.000000
 Q3 = 0.000000
 Q4 = 181.431380
Tempo total de execução: 0.021296 s
Comentários: O código acima foi de fácil de implementação e foi otimizado ao máximo, desde o método em si (GAUSS SEIDEL) até o cálculo dos fluxos depotências entre as barras como pode ser observado. Foram necessárias 15 iterações para que o erro calculado ficasse abaixo de 10e-5 (0,0001) assim como solicitado no enunciado do problema. Observar que para o cálculo deste erro assim como solicitado está sendo calculado a cada iteração a diferença entre as potências (ativa e reativa), ou seja, o programa irá convergir quando o maior erro (diferença entre a potência calculada na iteração e a potência da iteração anterior) for menor que 10e-5.
�
MODELAGEM NO POWERWORLD:
Figura 1 – Modelagem do Sistema no PowerWorld
Resultados:
Tabela 1 – Resultados POWERWORLD
	Bus
	 
	 
	 
	 
	 
	 
	 
	 
	Number
	Nom kV
	PU Volt
	Volt (kV)
	Angle (Deg)
	Load MW
	Load Mvar
	Gen MW
	Gen Mvar
	1
	230
	1
	230
	0
	50
	30,99
	186,81
	114,5
	2
	230
	0,98242
	225,957
	-0,98
	170
	105,35
	 
	 
	3
	230
	0,96901
	222,871
	-1,87
	200
	123,94
	 
	 
	4
	230
	1,02
	234,6
	1,52
	80
	49,58
	318
	181,43
�
RESULTADOS:
Tabela 2 – Resultados Matlab vs Resultado PowerWolrd
	SOFTWARE
	TENSÕES FINAIS NAS BARRAS (PU)
	
	V1
	V2
	V3
	V4
	MATLAB
	1.000000/_0.000000°
	0.982421/_-0.976235°
	0.969005/_-1.872250°
	1.020000/_1.522949°
	POWERWOLRD
	1,00/_0°
	0,98242/_-0,98°
	0,96901/_-1,87°
	1,02/_1,52°
	SOFTWARE
	CORRENTES ENTRANDO E SAINDO EM CADA BARRA (A)
	
	I12/ I21
	I13/ I31
	I24/ I42
	I34/ 43
	MATLAB
	112.106359/_29.952375°
126.614300/_-139.947930°
	290.303258/_31.957700°
300.624852/_-144.913443°
	385.769726/_-149.624664°
376.209021/_27.823683°
	309.084272/_-147.730699°
293.399127/_27.000920°
	POWERWORLD
	112,10
126,61
	290,29
300,62
	385,77
376,21
	309,08
293,40
	SOFTAWARE
	PERDAS ATIVAS E REATIVAS NA LINHA (MW E MVAr) – Considerando Yshunt
	
	P12/Q12
	P13/Q13
	P24/Q24
	P34/Q34
	MATLAB
	PLinha12 = 0.226733 QLinha12 = -8.937735
	PLinha13 = 1.031482 QLinha13 = -2.356099
	PLinha24 = 1.715483 QLinha24 = 0.805906
	PLinha34 = 1.835435 QLinha34 = -3.441307
	
	Ptotal = 4.809133 Qtotal= -13.929235
	POWERWOLRD
	PLinha12 = 0,23
QLinha12 = -8,94
	PLinha13=1,03
QLinha13=-2,36
	PLinha24=1,72
QLinha24=0,81
	PLinha34 = 1,84
QLinha34=-3,44 
	
	Ptotal = 4,82 Qtotal= -13,93
	SOFTWARE
	POTENCIA GERADA EM CADA BARRA (MW E MVAr)
	
	BARRA 1
	BARRA 2
	BARRA 3
	BARRA 4
	MATLAB
	P1 = 186.816033
Q1 = 114.499840
	P2 = 0.000000
Q2 = 0.000000
	P3 = 0.000000
Q3 = 0.000000
	P4 = 317.992645
Q4 = 181.43138
	POWERWORLD
	P1=186,81
Q1=114,5
	P2 = 0
Q2 = 0
	P3 = 0
Q3 = 0
	P4=318
Q4=181,43
Comentários: Os resultados da implementação do método no MATLAB e da modelagem no POWERWORLD podem ser facilmente comparados na tabela 2, conforme pode ser observado todos os valores estão muito próximos, dessa forma pode ser afirmado que os resultados são satisfatórios e evidenciam a correta implementação do método GAUSS SEIDEL no MATLAB e modelagem no software POWERWORLD. Como observação pode ser visto que o balanço de reativo nas linhas fica negativo quando considerado o Yshunt que significa a geração de reativo nas linhas, enquanto que quando não é considerado o Yshunt fica positivo, ou seja, há consumo de energia reativa, esses valores podem ser vistos no resultado do MATLAB. 
Questão 2- 
IMPLEMENTAÇÃO MATLAB:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programa de Pós Graduação em Engenharia Elétrica
% Disciplina Análise de Sistemas de Energia 2016.1
% Método Gauss-Seidel
% Questão 2
% Prof.: Niraldo
% Alunos: 1 - André Duarte
% 2 - Daniel Vasconcelos
% 3 - Humberto Rocha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc;
clear;
 
tic;
nb = 4; %quantidade de barras
nvc = 4; %número da barra de tensão controlada;
e = 1; %tolerância
 
% Bases adotadas
Sb = 100; % 100 MVA
Vb = 230; % 230 kV
 
% Montando a matriz admitância de barra
Yb = zeros(nb,nb);
Yb(1,1) = 8.985190 - 44.835953*i;
Yb(2,1) = -3.815629 + 19.078144*i;
Yb(3,1) = -5.169561 + 25.847809*i;
Yb(1,2) = Yb(2,1);
Yb(2,2) = Yb(1,1);
Yb(4,2) = Yb(3,1);
Yb(1,3) = Yb(4,2);
Yb(3,3) = 8.193267 - 40.863838*i;
Yb(4,3) = -3.023705 + 15.118528*i;
Yb(2,4) = Yb(3,1);
Yb(3,4) = Yb(4,3);
Yb(4,4) = Yb(3,3);
 
% Vetor de tensões iniciais nas barras
Vo = [1; 1; 1; 1.02];
 
% Vetor de Potência Ativa e Reativa
Sc = zeros(nb,1);
Sg = zeros(nb,1);
 
% Barra 1 é a de referência
Sc(1) = (-50-30.99*i)/Sb;
% Barra 2 e 3 de carga
Sc(2) = (-170-105.35*i)/Sb;
Sc(3) = (-200-123.94*i)/Sb;
% Barra 4 é a de tensão controlada
Sc(4) = (-80-49.58*i)/Sb;
Sg(4) = 318/Sb;
 
P = real(Sc + Sg);
Q = imag(Sc + Sg);
 
% Metódo Gauss-Seidel
V = Vo;
k= 0;
Sg1o = zeros(nb,1);
while (e >= 10e-5)
 k = k+1;
 Icalc = zeros(nb,1);
 for m=2:nb
 soma = 0;
 for j=1:nb;
 if (m ~= j) 
 soma = soma + Yb(m,j)*V(j);
 end
 Icalc(m) = Yb(m,j)*V(j)+Icalc(m);
 end 
 if (m == nvc) % Calculo de Q e V para a barra PV
 Qcalc(k) = -imag(conj(V(m))*Icalc(m));
 if Qcalc(k)> (1.25+imag(Sc(4))) % Barra passa a ser PQ devido a limitação de reativo
 Qcalc(k) = (1.25+imag(Sc(4)));
 Icalc(m) = (P(m)-i*Qcalc(k))/conj(V(m));
 end
 V(m) = ((P(m)-i*Qcalc(k))/conj(V(m)) - soma)/ Yb(m,m);
 if abs(V(m))>= abs(Vo(m)) %Verifica se a Barra pode voltar a ser PV
 V(m) = abs(Vo(4))*V(m)/abs(V(m)); % Corrigindo o módulo da tensão da barra PV 
 end
 else
 V(m) = ((P(m)-i*Q(m))/conj(V(m)) - soma)/ Yb(m,m); 
 end 
 end
 
 % Calculando o valor da corrente e potência gerada na Barra de Referência - Barra 1
 Icalc(1) = 0;
 for j=1:nb;
 Icalc(1) = Yb(1,j)*V(j)+Icalc(1);
 end
 Sg(1) = V(1)*conj(Icalc(1))-Sc(1);
 
 %Potência gerada na Barra 4
 Sg(4) = V(4)*conj(Icalc(4))-Sc(4); 
 
 Sg = [Sg(1); 0; 0; Sg(4)]; %Vetor com potência gerada em cada barra
 er = max(abs(real((Sg - Sg1o)))); % Cálculo do erro entre as Potências Calculadas na iteração
 ei = max(abs(imag((Sg - Sg1o))));
 e = max(er,ei);
 Sg1o = Sg;
end
 
% Cálculo dos Fluxos de Potência nas Linhas
Yshunt = zeros(nb);
Yshunt(1,2) = 0.05125*i;
Yshunt(2,1) = Yshunt(1,2);
Yshunt(1,3) = 0.03875*i;
Yshunt(3,1) = Yshunt(1,3);
Yshunt(2,4) = Yshunt(1,3);
Yshunt(4,2) = Yshunt(2,4);
Yshunt(3,4) = 0.06375*i;
Yshunt(4,3) = Yshunt(3,4);
for m=1:nb
 for j=1:nb
 if (m ~= j)
 SL(j,m) = V(j)*(conj((V(j)-V(m))*(-Yb(j,m)))); % Potencia na linha
 S(j,m) = SL(j,m)- abs(V(j))*(abs(V(j))*Yshunt(j,m)); % Potencia na saida e entrada da barra
 I(j,m) = (S(j,m)/V(j))*Sb*10^3/(sqrt(3)*Vb); % Corrente na saída e entrada da barra
 end
 end
end
 
Perdas = sum((Sg+Sc))*Sb;
SL = SL*Sb; % Potência nas Linhas em MVA
S = S*Sb;
SLp = [SL(1,2)+SL(2,1);SL(1,3)+SL(3,1);SL(2,4)+SL(4,2);SL(3,4)+SL(4,3)]; %Perdas nas linhas
Sp = [S(1,2)+S(2,1);S(1,3)+S(3,1);S(2,4)+S(4,2);S(3,4)+S(4,3)]; %Perdas na Linhas (Considerando o Shunt)
SLpr = real(SLp); %Perdas Ativas em MW na linha
SLpi = imag(SLp); %Perdas Reativas em MVar na linha
Spi = imag(Sp); % Perdas Reativas em MVar considerando o Yshunt
Sg = Sg*Sb; % Potência gerada por barra em MVA
Stg = sum(Sg)*Sb; % Potência total gerada no sistema em MVA
Stc = sum(Sc)*Sb; % Potência total das cargas do sistema em MVA
t1=toc;
 
fprintf(1,' \n'); 
fprintf(1,'**********************************************\n');
fprintf(1,'*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***\n');
fprintf(1,'**********************************************\n');
fprintf(1,' \n'); 
fprintf(1,'Quantidade de iterações: %d \n\n',k); 
fprintf(1,'Tensões Finais nas barras em PU:\n V1 = %f/_%f°\n V2 = %f/_%f°\n V3 = %f/_%f°\n V4 =%f/_%f°\n\n',abs(V(1)),rad2deg(angle(V(1))),abs(V(2)),rad2deg(angle(V(2))),abs(V(3)),rad2deg(angle(V(3))),abs(V(4)),rad2deg(angle(V(4)))); 
fprintf(1,'Correntes entrando e saindo em cada Barra em A:\n I12 = %f/_%f° I21 = %f/_%f°\n I13 = %f/_%f° I31 = %f/_%f°\n I24 = %f/_%f° I42 = %f/_%f°\n I34 = %f/_%f° I43 = %f/_%f°\n\n',abs(I(1,2)),rad2deg(angle(I(1,2))),abs(I(2,1)),rad2deg(angle(I(2,1))),abs(I(1,3)),rad2deg(angle(I(1,3))),abs(I(3,1)),rad2deg(angle(I(3,1))),abs(I(2,4)),rad2deg(angle(I(2,4))),abs(I(4,2)),rad2deg(angle(I(4,2))),abs(I(3,4)),rad2deg(angle(I(3,4))),abs(I(4,3)),rad2deg(angle(I(4,3)))); 
fprintf(1,'Perdas Ativa e Reativa em cada Linha em MW e MVAr:\n PLinha12 = %f QLinha12(sem Yshunt) = %f QLinha12(com Yshunt) = %f\n PLinha13 = %f QLinha13(sem Yshunt) = %f QLinha13(com Yshunt) = %f\n PLinha24 = %f QLinha24(sem Yshunt) = %f QLinha24(com Yshunt) = %f\n PLinha34 = %f QLinha34(sem Yshunt) = %f QLinha34(com Yshunt) = %f\n\n',SLpr(1),SLpi(1),Spi(1),SLpr(2),SLpi(2),Spi(2),SLpr(3),SLpi(3),Spi(3),SLpr(4),SLpi(4),Spi(4));
fprintf(1,'Perda total no Sistema (Linhas) em MW e MVAr:\n Ptotal = %f\n Qtotal(sem Yshunt)= %f Qtotal(com Yshunt)= %f\n\n', sum(SLpr), sum(SLpi), sum(Spi)); 
fprintf(1,'Potência Ativa Gerada nas barras em MW:\n P1 = %f\n P2 = %f\n P3 = %f\n P4 = %f\n\n',real(Sg(1)),real(Sg(2)),real(Sg(3)),real(Sg(4)));
fprintf(1,'Potência Reativa Gerada nas barras em MVAr:\n Q1 = %f\n Q2 = %f\n Q3 = %f\n Q4 = %f\n\n',imag(Sg(1)),imag(Sg(2)),imag(Sg(3)),imag(Sg(4)));
fprintf(1,'Tempo total de execução: %f s\n', t1);
 
**********************************************
*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***
**********************************************
 
Quantidade de iterações: 15 
Tensões Finais nas barras em PU:
 V1 = 1.000000/_0.000000°
 V2 = 0.966770/_-0.803606°
 V3 = 0.958827/_-1.776892°
 V4 = 0.993838/_1.941819°
Correntes entrando e saindo em cada Barra em A:
 I12 = 165.127792/_54.088471° I21 = 186.121288/_-121.294698°
 I13 = 331.941851/_41.934864° I31 = 344.756390/_-135.677250°
 I24 = 361.575196/_-160.508333° I42 = 355.480761/_16.604278°
 I34 = 284.666621/_-159.505777° I43 = 275.229821/_14.397272°
Perdas Ativa e Reativa em cada Linha em MW e MVAr:
 PLinha12 = 0.493887 QLinha12(sem Yshunt) = 2.469436 QLinha12(com Yshunt) = -7.445621
 PLinha13 = 1.353062 QLinha13(sem Yshunt) = 6.765310 QLinha13(com Yshunt) = -0.672165
 PLinha24 = 1.518960 QLinha24(sem Yshunt) = 7.594802 QLinha24(com Yshunt) = 0.145660
 PLinha34 = 1.583990 QLinha34(sem Yshunt) = 7.919949 QLinha34(com Yshunt) = -4.237575
Perda total no Sistema (Linhas) em MW e MVAr:
 Ptotal = 4.949899
 Qtotal(sem Yshunt)= 24.749497 Qtotal(com Yshunt)= -12.209701
Potência Ativa Gerada nas barras em MW:
 P1 = 186.954808
 P2 = 0.000000
 P3 = 0.000000
 P4 = 317.999352
Potência Reativa Gerada nas barras em MVAr:
 Q1 = 172.640125
 Q2 = 0.000000
 Q3 = 0.000000
 Q4 = 125.000285
Tempo total de execução: 0.018615 s
Comentários: A diferença entre a implementação da questão 1 e a questão está nas linhas de código a seguir que foram inseridas com o objetivo de limitar o reativo gerado na barra 4, quando o Qcalc chega nesse limitador a barra PV passa a ser uma barra PQ, o Icalc é calculado novamente para essa limitação de reativo e em seguida a tensão para a barra 4 é calculada novamente. A barra 4 somente voltará a ser uma barra PV se a tensão calculada em alguma iteração for superior a tensão inicial, ou seja, estaria com reativo sobrando, mas não foi o caso desse sistema implementado.
if Qcalc(k)> (1.25+imag(Sc(4))) % Barra passa a ser PQ devido a limitação de reativo
 Qcalc(k) = (1.25+imag(Sc(4)));
 Icalc(m) = (P(m)-i*Qcalc(k))/conj(V(m));
 end
 V(m) = ((P(m)-i*Qcalc(k))/conj(V(m)) - soma)/ Yb(m,m);
 if abs(V(m))>= abs(Vo(m)) %Verifica se a Barra pode voltar a ser PV
 V(m) = abs(Vo(4))*V(m)/abs(V(m)); % Corrigindo o módulo da tensão da barra PV 
 end
�
MODELAGEM NO POWERWORLD:
Figura 2 – Modelagem do Sistema no PowerWorld
Resultados:
Tabela 3 – Resultados POWERWORLD
	Bus
	 
	 
	 
	 
	 
	 
	 
	 
	Number
	Nom kV
	PU Volt
	Volt (kV)
	Angle (Deg)
	Load MW
	Load Mvar
	Gen MW
	Gen Mvar
	1
	230
	1
	230
	0
	50
	30,99
	186,95
	172,65
	2
	230
	0,96677
	222,357
	-0,8
	170
	105,35
	 
	 
	3
	230
	0,95882
	220,53
	-1,78
	200
	123,94
	 
	 
	4
	230
	0,99384
	228,582
	1,94
	80
	49,58
	318
	125
�
RESULTADOS:
Tabela 4 – Resultados Matlab vs Resultado Power Wolrd
	SOFTWARE
	TENSÕES FINAIS NAS BARRAS (PU)
	
	V1
	V2
	V3
	V4
	MATLAB
	1.000000/_0.000000°
	0.966770/_-0.803606°
	0.958827/_-1.776892°
	0.993838/_1.941819°
	POWERWOLRD
	1,00/_0°
	0,96677/_-0,8°
	0,95882/_-1,78°
	0,99384/_1,94
	SOFTWARE
	CORRENTES ENTRANDO E SAINDO EM CADA BARRA (A)
	
	I12/ I21
	I13/ I31
	I24/ I42
	I34/ 43
	MATLAB
	165.127792/_54.088471°
186.121288/_-121.294698°
	331.941851/_41.934864°
344.756390/_-135.677250°
	361.575196/_-160.508333°
355.480761/_16.604278°
	284.666621/_-159.505777°
275.229821/_14.397272°
	POWERWORLD
	165,14
186,13
	331,95
344,76
	361,58
355,48
	284,67
275,23
	SOFTAWARE
	PERDAS ATIVAS E REATIVAS NA LINHA (MW E MVAr) – Considerando Yshunt
	
	P12/Q12
	P13/Q13
	P24/Q24
	P34/Q34
	MATLAB
	PLinha12 = 0.493887
QLinha12 = -7.445621
	PLinha13 = 1.353062
QLinha13 = -0.672165
	PLinha24 = 1.518960
QLinha24 = 0.145660
	PLinha34 = 1.583990
QLinha34 = -4.237575
	
	Ptotal = 4.949899 Qtotal= -12.209701
	POWERWOLRD
	PLinha12 = 0,49
QLinha12 = -7,45
	PLinha13=1,35
QLinha13=-0,67
	PLinha24=1,52
QLinha24=0,15
	PLinha34 = 1,58
QLinha34=-4,24 
	
	Ptotal = 4,94 Qtotal= -12,21
	SOFTWARE
	POTENCIA GERADA EM CADA BARRA (MW E MVAr)
	
	BARRA 1
	BARRA 2
	BARRA 3
	BARRA 4
	MATLAB
	P1 = 186.954808
Q1 = 172.640125
	P2 = 0.000000
Q2 = 0.000000
	P3 = 0.000000
Q3 = 0.000000
	P4 = 317.999352
Q4 = 125.000285
	POWERWORLD
	P1=186,95
Q1=172,65
	P2 = 0
Q2 = 0
	P3 = 0
Q3 = 0
	P4=318
Q4=125
Comentários: Nesta questão o reativo da barra PV (barra 4) ficou limitado em 125 MVAr, neste caso como o reativo necessário para manter a tensão da barra PV é superior ao limitado, essa barra passou a ser uma barra PQ, ela somente voltaria a ser uma barra PV caso o reativo necessário tivesse superado em algum momento a limitação de 125 MVAr que não foi o caso, por esse motivo a tensão final na barra 4 ficou abaixo da tensão inicial 1,02. Como pode ser observado na tabela 4, os resultados nos dois softwares estão bem próximos e dessa forma são satisfatórios e evidenciam a correta implementação do método no MATLAB e modelagem no POWERWOLRD.
QUESTÃO 3
IMPLEMENTAÇÃO MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programa de Pós Graduação em Engenharia Elétrica
% Disciplina Análise de Sistemas de Energia 2016.1
% Metódo Gauss-Seidel
% Questão 3
% Profº: Niraldo
% Alunos: 1 - André Duarte
% 2 - Daniel Vasconcelos
% 3 - Humberto Rocha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc;
clear;
 
tic;
nb = 4; %quantidade de barras
 
e = 1; %tolerância
 
% Bases adotadas
Sb = 100; % 100 MVA
Vb = 230; % 230 kV
 
% Montando a matriz admitância de barra
Yb = zeros(nb,nb);
Yb(1,1) = 8.985190 - 44.835953*i;
Yb(2,1) = -3.815629 + 19.078144*i;
Yb(3,1) = -5.169561 + 25.847809*i;
Yb(1,2) = Yb(2,1);
Yb(2,2) = Yb(1,1);
Yb(4,2) = Yb(3,1);
Yb(1,3) = Yb(4,2);
Yb(3,3) = 8.193267 - 40.863838*i;
Yb(4,3) = -3.023705 + 15.118528*i;
Yb(2,4) = Yb(3,1);
Yb(3,4) = Yb(4,3);
Yb(4,4) = Yb(3,3);
 
% Vetor de tensões iniciais nas barras
Vo = [1; 0.99; 1; 1.02];
 
% Vetor de Potência Ativa e Reativa
Sc = zeros(nb,1);
Sg = zeros(nb,1);
 
% Barra 1 é ade referência
Sc(1) = (-50-30.99*i)/Sb;
% Barra 2 e 3 de carga
Sc(2) = (-170-105.35*i)/Sb;
Sc(3) = (-200-123.94*i)/Sb;
% Barra 4 é a de tensão controlada
Sc(4) = (-80-49.58*i)/Sb;
Sg(4) = 318/Sb;
 
P = real(Sc + Sg);
Q = imag(Sc + Sg);
 
% Metódo Gauss-Seidel
V = Vo;
k= 0;
Sg1o = zeros(nb,1);
while (e >= 10e-5)
 k = k+1;
 Icalc = zeros(nb,1);
 for m=2:nb
 soma = 0;
 for j=1:nb;
 if (m ~= j) 
 soma = soma + Yb(m,j)*V(j);
 end
 Icalc(m) = Yb(m,j)*V(j)+Icalc(m);
 end 
 if (m == 2) || (m==4) % Calculo de Q e V para a barra PV
 Qcalc(k,m) = -imag(conj(V(m))*Icalc(m));
 V(m) = ((P(m)-i*Qcalc(k,m))/conj(V(m)) - soma)/ Yb(m,m);
 V(m) = abs(Vo(m))*V(m)/abs(V(m)); % Corrigindo o módulo da tensão da barra PV
 else
 V(m) = ((P(m)-i*Q(m))/conj(V(m)) - soma)/ Yb(m,m); 
 end 
 end
 
 % Calculando o valor da corrente e potência gerada na Barra de Referência - Barra 1
 Icalc(1) = 0;
 for j=1:nb;
 Icalc(1) = Yb(1,j)*V(j)+Icalc(1);
 end
 Sg(1) = V(1)*conj(Icalc(1))-Sc(1);
 
 %Potência gerada na Barra 2 e 4
 Sg(2) = i*imag(V(2)*conj(Icalc(2))-Sc(2)); 
 Sg(4) = V(4)*conj(Icalc(4))-Sc(4); 
 
 Sg = [Sg(1); Sg(2); Sg(3); Sg(4)]; %Vetor com potência gerada em cada barra
 er = max(abs(real((Sg - Sg1o)))); % Cálculo do erro entre as Potências Calculadas na iteração
 ei = max(abs(imag((Sg - Sg1o))));
 e = max(er,ei);
 Sg1o = Sg;
end
 
% Cálculo dos Fluxos de Potência nas Linhas
Yshunt = zeros(nb);
Yshunt(1,2) = 0.05125*i;
Yshunt(2,1) = Yshunt(1,2);
Yshunt(1,3) = 0.03875*i;
Yshunt(3,1) = Yshunt(1,3);
Yshunt(2,4) = Yshunt(1,3);
Yshunt(4,2) = Yshunt(2,4);
Yshunt(3,4) = 0.06375*i;
Yshunt(4,3) = Yshunt(3,4);
for m=1:nb
 for j=1:nb
 if (m ~= j)
 SL(j,m) = V(j)*(conj((V(j)-V(m))*(-Yb(j,m)))); % Potencia na linha
 S(j,m) = SL(j,m)- abs(V(j))*(abs(V(j))*Yshunt(j,m)); % Potencia na saida e entrada da barra
 I(j,m) = (S(j,m)/V(j))*Sb*10^3/(sqrt(3)*Vb); % Corrente na saída e entrada da barra
 end
 end
end
 
Perdas = sum((Sg+Sc))*Sb;
SL = SL*Sb; % Potência nas Linhas em MVA
S = S*Sb;
SLp = [SL(1,2)+SL(2,1);SL(1,3)+SL(3,1);SL(2,4)+SL(4,2);SL(3,4)+SL(4,3)]; %Perdas nas linhas (Não considera o shunt)
Sp = [S(1,2)+S(2,1);S(1,3)+S(3,1);S(2,4)+S(4,2);S(3,4)+S(4,3)]; %Perdas na Linhas (Considerando o Shunt)
SLpr = real(SLp); %Perdas Ativas em MW na linha
SLpi = imag(SLp); %Perdas Reativas em MVar na linha sem o Yshunt
Spi = imag(Sp); % Perdas Reativas em MVar considerando o Yshunt
Sg = Sg*Sb; % Potência gerada por barra em MVA
Stg = sum(Sg)*Sb; % Potência total gerada no sistema em MVA
Stc = sum(Sc)*Sb; % Potência total das cargas do sistema em MVA
t1=toc;
 
fprintf(1,' \n'); 
fprintf(1,'**********************************************\n');
fprintf(1,'*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***\n');
fprintf(1,'**********************************************\n');
fprintf(1,' \n'); 
fprintf(1,'Quantidade de iterações: %d \n\n',k); 
fprintf(1,'Tensões Finais nas barras em PU:\n V1 = %f/_%f°\n V2 = %f/_%f°\n V3 = %f/_%f°\n V4 = %f/_%f°\n\n',abs(V(1)),rad2deg(angle(V(1))),abs(V(2)),rad2deg(angle(V(2))),abs(V(3)),rad2deg(angle(V(3))),abs(V(4)),rad2deg(angle(V(4)))); 
fprintf(1,'Correntes entrando e saindo em cada Barra em A:\n I12 = %f/_%f° I21 = %f/_%f°\n I13 = %f/_%f° I31 = %f/_%f°\n I24 = %f/_%f° I42 = %f/_%f°\n I34 = %f/_%f° I43 = %f/_%f°\n\n',abs(I(1,2)),rad2deg(angle(I(1,2))),abs(I(2,1)),rad2deg(angle(I(2,1))),abs(I(1,3)),rad2deg(angle(I(1,3))),abs(I(3,1)),rad2deg(angle(I(3,1))),abs(I(2,4)),rad2deg(angle(I(2,4))),abs(I(4,2)),rad2deg(angle(I(4,2))),abs(I(3,4)),rad2deg(angle(I(3,4))),abs(I(4,3)),rad2deg(angle(I(4,3)))); 
fprintf(1,'Perdas Ativa e Reativa em cada Linha em MW e MVAr:\n PLinha12 = %f QLinha12(sem Yshunt) = %f QLinha12(com Yshunt) = %f\n PLinha13 = %f QLinha13(sem Yshunt) = %f QLinha13(com Yshunt) = %f\n PLinha24 = %f QLinha24(sem Yshunt) = %f QLinha24(com Yshunt) = %f\n PLinha34 = %f QLinha34(sem Yshunt) = %f QLinha34(com Yshunt) = %f\n\n',SLpr(1),SLpi(1),Spi(1),SLpr(2),SLpi(2),Spi(2),SLpr(3),SLpi(3),Spi(3),SLpr(4),SLpi(4),Spi(4));
fprintf(1,'Perda total no Sistema (Linhas) em MW e MVAr:\n Ptotal = %f\n Qtotal(sem Yshunt)= %f Qtotal(com Yshunt)= %f\n\n', sum(SLpr), sum(SLpi), sum(Spi)); 
fprintf(1,'Potência Ativa Gerada nas barras em MW:\n P1 = %f\n P2 = %f\n P3 = %f\n P4 = %f\n\n',real(Sg(1)),real(Sg(2)),real(Sg(3)),real(Sg(4)));
fprintf(1,'Potência Reativa Gerada nas barras em MVAr:\n Q1 = %f\n Q2 = %f\n Q3 = %f\n Q4 = %f\n\n',imag(Sg(1)),imag(Sg(2)),imag(Sg(3)),imag(Sg(4)));
fprintf(1,'Tempo total de execução: %f s\n', t1);
**********************************************
*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***
**********************************************
 
Quantidade de iterações: 15 
Tensões Finais nas barras em PU:
 V1 = 1.000000/_0.000000°
 V2 = 0.990000/_-1.047184°
 V3 = 0.969005/_-1.873075°
 V4 = 1.020000/_1.520761°
Correntes entrando e saindo em cada Barra em A:
 I12 = 98.141580/_10.858114° I21 = 105.777084/_-155.368025°
 I13 = 290.371492/_31.945413° I31 = 300.689474/_-144.925985°
 I24 = 361.623976/_-156.462520° I42 = 354.189959/_20.643035°
 I34 = 309.017425/_-147.717678° I43 = 293.326861/_27.013303°
Perdas Ativa e Reativa em cada Linha em MW e MVAr:
 PLinha12 = 0.164336 QLinha12(sem Yshunt) = 0.821680 QLinha12(com Yshunt) = -9.326332
 PLinha13 = 1.031946 QLinha13(sem Yshunt) = 5.159730 QLinha13(com Yshunt) = -2.353781
 PLinha24 = 1.513698 QLinha24(sem Yshunt) = 7.568492 QLinha24(com Yshunt) = -0.260946
 PLinha34 = 1.834587 QLinha34(sem Yshunt) = 9.172936 QLinha34(com Yshunt) = -3.445552
Perda total no Sistema (Linhas) em MW e MVAr:
 Ptotal = 4.544567
 Qtotal(sem Yshunt)= 22.722838 Qtotal(com Yshunt)= -15.386611
Potência Ativa Gerada nas barras em MW:
 P1 = 186.553846
 P2 = 0.000000
 P3 = 0.000000
 P4 = 317.990590
Potência Reativa Gerada nas barras em MVAr:
 Q1 = 99.560300
 Q2 = 34.103332
 Q3 = 0.000000
 Q4 = 160.811881
Tempo total de execução: 0.017724 s
COMENTÁRIOS: Em relação ao código da questão 1 neste foi alterado as linhas de código a fim de que a barra 2 seja uma barra PV e os cálculos de reativos sejam feitos a cada iteração a fim de manter a tensão na barra constante em 0,99 pu. Neste caso temos duas barras PV (Barras 2 e 4). Linha de código alterada:
 if (m == 2) || (m==4) % Calculo de Q e V para a barra PV
�
MODELAGEM NO POWERWORLD:
Figura 3 – Modelagem do Sistema no PowerWorld
Resultados:
Tabela 5 – Resultados POWERWORLD
	Bus
	 
	 
	 
	 
	 
	 
	 
	 
	Number
	Nom kV
	PU Volt
	Volt (kV)
	Angle (Deg)
	Load MW
	Load Mvar
	Gen MW
	Gen Mvar
	1
	230
	1
	230
	0
	50
	30,99
	186,54
	99,56
	2
	230
	0,99
	227,7
	-1,05
	170
	105,35
	0
	34,11
	3
	230
	0,96901
	222,871
	-1,87
	200
	123,94
	 
	 
	4
	230
	1,02
	234,6
	1,52
	80
	49,58
	318
	160,81
�
RESULTADOS:
Tabela 6 – Resultados Matlab vs Resultado Power Wolrd
	SOFTWARE
	TENSÕES FINAIS NAS BARRAS (PU)
	
	V1
	V2
	V3
	V4
	MATLAB
	1.000000/_0.000000°
	0.990000/_-1.047184°
	0.969005/_-1.873075°
	1.020000/_1.520761°
	POWERWOLRD
	1,00/_0°
	0,99/_-1,05°
	0,96901/_-1,87°
	1,02/_1,52°
	SOFTWARE
	CORRENTES ENTRANDO E SAINDO EM CADA BARRA (A)
	
	I12/ I21
	I13/ I31
	I24/ I42
	I34/ 43
	MATLAB
	98.141580/_10.858114°
105.777084/_-155.368025°
	290.371492/_31.945413°
300.689474/_-144.925985°
	361.623976/_-156.462520°
354.189959/_20.643035°
	309.017425/_-147.717678°
293.326861/_27.013303°
	POWERWORLD
	98,13
105,76
	290,36
300,68
	361,62
354,19
	309,02
293,33
	SOFTAWARE
	PERDAS ATIVAS E REATIVAS NA LINHA (MW E MVAr) – ConsiderandoYshunt
	
	P12/Q12
	P13/Q13
	P24/Q24
	P34/Q34
	MATLAB
	PLinha12 = 0.164336
QLinha12 = -9.326332
	PLinha13 = 1.031946
QLinha13 = -2.353781
	PLinha24 = 1.513698
QLinha24 = -0.260946
	PLinha34 = 1.834587
QLinha34 = -3.445552
	
	Ptotal = 4.544567 Qtotal= -15.386611
	POWERWOLRD
	PLinha12 = 0,16
QLinha12 = -9,33
	PLinha13=1,03
QLinha13=-2,35
	PLinha24=1,51
QLinha24=-0,26
	PLinha34 = 1,83
QLinha34=-3,45 
	
	Ptotal = 4,53 Qtotal= -15,39
	SOFTWARE
	POTENCIA GERADA EM CADA BARRA (MW E MVAr)
	
	BARRA 1
	BARRA 2
	BARRA 3
	BARRA 4
	MATLAB
	P1 = 186.553846
Q1 = 99.560300
	P2 = 0.000000
Q2 = 34.103332
	P3 = 0.000000
Q3 = 0.000000
	P4 = 317.990590
Q4 = 160.811881
	POWERWORLD
	P1=186,54
Q1=99,56
	P2 = 0
Q2 = 34,11
	P3 = 0
Q3 = 0
	P4=318
Q4=160,81
Comentários: Nesta questão a barra 2 também passa a ser uma barra PV e há um gerador injetando reativo na barra a fim de manter a tensão constante em 0,99 pu, o resultado pode ser observado na tabela 6 acima, onde os resultados nos dois softwares estão bem próximos e dessa forma são satisfatórios e evidenciam a correta implementação do método no MATLAB e modelagem no POWERWOLRD. Pode ser observado também que as perdas ativas total no sistema são menores em relação as questões 1 e 2, isso se deve pelo fato do menor carregamento nas linhas.
QUESTÃO 4
A Figura 4-1 a seguir ilustra o transformador ideal com tap “t” do lado do nó “i” com as correntes Ii e Ij entrando nas duas barras e as tensões Vi e Vj referentes ao nó de referência. 
Figura 4 – Transformador Ideal com tap no lado do nó “i”
As expressões da potência complexa das barras i e j são dadas pela equações 4.1 e 4.2 
 	(4.1)
 	(4.2)
Como o transformador é ideal, então podemos considerar que a potência complexa no lado do primário do transformador é igual à potência complexa do lado secundário. Assim:
 (4.3)
Das equações (4.1), (4.2) e (4.3) obtemos:
 (4.4)
A equação (4.4) pode ser reescrita como:
 (4.5)
A corrente 
 (primário do transformador) pode ser expressa pela equação (4.6).
 
 (4.6)
Multiplicando por 
 a equação (4.6), temos a equação (4.7), que representa a corrente no primário do transformador em função das tensões do primário e secundário:
	(4.7)
Sabendo-se que a equação matricial que representa o quadripolo de admitâncias é dado pela equação (4.8):
 (4.8)
Definindo
e rearranjando as equações (4.6) e (4.7), temos:
 (4.9)
O circuito equivalente π somente poderá ser determinado caso “t” nas equações anteriores seja real. Isto implica em Yij = Yji e o circuito π com admitâncias encontra-se representado na Figura 4-2. 
Figura 5 – Circuito π com admitâncias conforme eq. (4.9), para “t” real. 
�
QUESTÃO 5:
IMPLEMENTAÇÃO MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programa de Pós Graduação em Engenharia Elétrica
% Disciplina Análise de Sistemas de Energia 2016.1
% Metódo Gauss-Seidel
% Questão 5
% Profº: Niraldo
% Alunos: 1 - André Duarte
% 2 - Daniel Vasconcelos
% 3 - Humberto Rocha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc;
clear;
 
tic;
nb = 5; %quantidade de barras
 
e = 1; %tolerância
 
% Bases adotadas
Sb = 100; % 100 MVA
Vb = 230; % 230 kV
t = 1; % Valor inicial do TAP
Yt = 1/(0.02*i);
 
% Montando a matriz admitância de barra
Yb = zeros(nb,nb);
Yb(1,1) = 8.985190 - 44.835953*i;
Yb(2,1) = -3.815629 + 19.078144*i;
Yb(3,1) = -5.169561 + 25.847809*i;
Yb(1,2) = Yb(2,1);
Yb(2,2) = Yb(1,1);
Yb(4,2) = Yb(3,1);
Yb(1,3) = Yb(4,2);
Yb(3,3) = (8.193267 - 40.863838*i)+(t^2)*Yt;
Yb(4,3) = -3.023705 + 15.118528*i;
Yb(2,4) = Yb(3,1);
Yb(3,4) = Yb(4,3);
Yb(4,4) = (8.193267 - 40.863838*i);
Yb(3,5) = -t*Yt;
Yb(5,3) = -t*Yt;
Yb(5,5) = Yt;
 
% Vetor de tensões iniciais nas barras
Vo = [1; 1; 1; 1.02; 0.97];
 
% Vetor de Potência Ativa e Reativa
Sc = zeros(nb,1);
Sg = zeros(nb,1);
 
% Barra 1 é a de referência
Sc(1) = (-50-30.99*i)/Sb;
% Barra 2, 3 e 5 de carga
Sc(2) = (-170-105.35*i)/Sb;
Sc(3) = 0;
Sc(5) = (-200-123.94*i)/Sb;
% Barra 4 é a de tensão controlada
Sc(4) = (-80-49.58*i)/Sb;
Sg(4) = 318/Sb;
 
P = real(Sc + Sg);
Q = imag(Sc + Sg);
 
% Metódo Gauss-Seidel
V = Vo;
k= 0;
Sg1o = zeros(nb,1);
while (e >= 10e-5)
 % Recalculando os valores da matriz Ybus que dependem de t
 
 Yb(3,3) = (8.193267 - 40.863838*i)+(t^2)*Yt;
 Yb(3,5) = -t*Yt;
 Yb(5,3) = -t*Yt;
 
 k = k+1;
 Icalc = zeros(nb,1);
 for m=2:nb
 soma = 0;
 for j=1:nb;
 if (m ~= j) 
 soma = soma + Yb(m,j)*V(j);
 end
 Icalc(m) = Yb(m,j)*V(j)+Icalc(m);
 end 
 if (m == 4) % Calculo de Q e V para as barras PV
 Qcalc(k) = -imag(conj(V(m))*Icalc(m));
 V(m) = ((P(m)-i*Qcalc(k))/conj(V(m)) - soma)/ Yb(m,m);
 V(m) = abs(Vo(m))*V(m)/(abs(V(m))); % Corrigindo o módulo da tensão da barra PV
 else
 V(m) = ((P(m)-i*Q(m))/conj(V(m)) - soma)/ Yb(m,m);
 end 
 end
 
 deltat = Vo(5)-abs(V(5)); % Calculo da variação do TAP proporcional a variação de tensão na barra
 t = t + deltat; % Valor do TAP do Transformador 
 
 % Calculando o valor da corrente e potência gerada na Barra de Referência - Barra 1
 Icalc(1) = 0;
 for j=1:nb;
 Icalc(1) = Yb(1,j)*V(j)+Icalc(1);
 end
 Sg(1) = V(1)*conj(Icalc(1))-Sc(1);
 
 %Potência gerada na Barra 4
 Sg(4) = V(4)*conj(Icalc(4))-Sc(4); 
 
 Sg = [Sg(1); 0; 0; Sg(4); 0]; %Vetor com potência gerada em cada barra
 er = max(abs(real((Sg - Sg1o)))); % Cálculo do erro entre as Potências Calculadas na iteração
 ei = max(abs(imag((Sg - Sg1o))));
 e = max(er,ei);
 Sg1o = Sg;
end
 
% Cálculo dos Fluxos de Potência nas Linhas
Yshunt = zeros(nb);
Yshunt(1,2) = 0.05125*i;
Yshunt(2,1) = Yshunt(1,2);
Yshunt(1,3) = 0.03875*i;
Yshunt(3,1) = Yshunt(1,3);
Yshunt(2,4) = Yshunt(1,3);
Yshunt(4,2) = Yshunt(2,4);
Yshunt(3,4) = 0.06375*i;
Yshunt(4,3) = Yshunt(3,4);
for m=1:nb
 for j=1:nb
 if (m ~= j)
 SL(j,m) = V(j)*(conj((V(j)-V(m))*(-Yb(j,m)))); % Potencia na linha
 S(j,m) = SL(j,m)- abs(V(j))*(abs(V(j))*Yshunt(j,m)); % Potencia na saida e entrada da barra
 I(j,m) = (S(j,m)/V(j))*Sb*10^3/(sqrt(3)*Vb); % Corrente na saída e entrada da barra
 end
 end
end
% Cálculo das correntes entre as barras 3 e 5 devido à presença do transformador regulador
Iload5to3 = (-t*Yt*V(3,1)+Yt*V(5,1))*(Sb*10^3/(sqrt(3)*Vb)); % Corrente no secundário do transformador.
Iload3to5 = (t*t*Yt*V(3,1)-t*Yt*V(5,1))*(Sb*10^3/(sqrt(3)*Vb)); % Corrente no primário do transformador.
 
Perdas = sum((Sg+Sc))*Sb;
SL = SL*Sb; % Potência nas Linhas em MVA
S = S*Sb;
Sg = Sg*Sb; % Potência gerada por barra em MVA
Stg = sum(Sg)*Sb; % Potência total gerada no sistema em MVA
Stc = sum(Sc)*Sb; % Potência total das cargas do sistema em MVA
t1=toc;
 
fprintf(1,' \n'); 
fprintf(1,'**********************************************\n');
fprintf(1,'*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***\n');
fprintf(1,'**********************************************\n');
fprintf(1,' \n'); 
fprintf(1,'Quantidade de iterações: %d \n\n',k); 
fprintf(1,'Tensões Finais nas barras em PU:\n V1 = %f/_%f°\n V2 = %f/_%f°\n V3 = %f/_%f°\n V4 = %f/_%f°\n V5 = %f/_%f°\n\n',abs(V(1)),rad2deg(angle(V(1))),abs(V(2)),rad2deg(angle(V(2))),abs(V(3)),rad2deg(angle(V(3))),abs(V(4)),rad2deg(angle(V(4))),abs(V(5)),rad2deg(angle(V(5)))); 
fprintf(1,'Correntes entrando e saindo em cada Barra em A:\n\n I12 = %f/_%f° I21 = %f/_%f°\n I13 = %f/_%f° I31 = %f/_%f°\n I24 = %f/_%f° I42 = %f/_%f°\n I34 = %f/_%f° I43 = %f/_%f°\nI35 = %f/_%f° I53 = %f/_%f°\n\n',abs(I(1,2)),rad2deg(angle(I(1,2))),abs(I(2,1)),rad2deg(angle(I(2,1))),abs(I(1,3)),rad2deg(angle(I(1,3))),abs(I(3,1)),rad2deg(angle(I(3,1))),abs(I(2,4)),rad2deg(angle(I(2,4))),abs(I(4,2)),rad2deg(angle(I(4,2))),abs(I(3,4)),rad2deg(angle(I(3,4))),abs(I(4,3)),rad2deg(angle(I(4,3))), abs(Iload3to5),rad2deg(angle(Iload3to5)),abs(Iload5to3),rad2deg(angle(Iload5to3))); 
fprintf(1,'Potência Ativa Gerada nas barras em MW:\n P1 = %f\n P2 = %f\n P3 = %f\n P4 = %f\n P5 = %f\n\n',real(Sg(1)),real(Sg(2)),real(Sg(3)),real(Sg(4)),real(Sg(5)));
fprintf(1,'Potência Reativa Gerada nas barras em MVAr:\n Q1 = %f\n Q2 = %f\n Q3 = %f\n Q4 = %f\n Q5 = %f\n\n',imag(Sg(1)),imag(Sg(2)),imag(Sg(3)),imag(Sg(4)),imag(Sg(5)));
fprintf(1,'TAP FINAL: %f\n\n', t);
fprintf(1,'Tempo total de execução: %f s\n', t1);
RESULTADO:
**********************************************
*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***
**********************************************
 
Quantidade de iterações: 28 
Tensões Finais nas barras em PU:
 V1 = 1.000000/_0.000000°
 V2 = 0.982421/_-0.975774°
 V3 = 0.966053/_-1.846793°
 V4 = 1.020000/_1.523655°
 V5 = 0.970000/_-4.218671°
Correntes entrando e saindo em cada Barra em A:
 I12 = 112.077087/_29.965204° I21 = 126.589961/_-139.934454°
 I13 = 301.543612/_35.100552° I31 = 312.689848/_-142.000106°
 I24 = 385.791013/_-149.628036° I42 = 376.231122/_27.820389°
 I34 = 315.607122/_-145.995497° I43 = 299.058688/_28.945472°
 I35 = 628.019513/_-36.005138° I53 = 608.887054/_143.994862°
Potência Ativa Gerada nas barras em MW:
 P1 = 186.960890
 P2 = 0.000000
 P3 = 0.000000
 P4 = 317.999982
 P5 = 0.000000
Potência Reativa Gerada nas barras em MVAr:
 Q1 = 122.364961
 Q2 = 0.000000
 Q3 = 0.000000
 Q4 = 186.117348
 Q5 = 0.000000
TAP FINAL: 1.031422
Tempo total de execução: 0.003059 s
COMENTÁRIOS: Neste programa foi inserido mais uma barra (barra 5) e a admitância do transformador na matriz Ybarra (Yb) entre as barras 3 e 5, na questão é informado que o valor da reatância deve ser de 0,2 pu, no entanto este valor não é um valor usual para transformadores mesmo assim foi simulado e encontrado um tap final de 1.377936, o que evidencia a inexequibilidade com este valor de reatância, dessa forma foi considerado uma reatância de 0,02 pu para a reatância deste transformador, tanto na questão 5 como na questão 6 deste trabalho. Como o problema estima o valor inicial do tap em 1 e pre-especifica o valor de |V5| em 0,97 o programa calcula em cada iteração o valor do tap necessário para manter a tensão na barra 5 fixa (0,97 pu), conforme literatura consultada (A. Monticelli, Fluxo de Carga em Redes de Energia Elétrica, Edgar Blucher, 1983) a variação do tap é proporcional a variação de tensão na barra 5 (Vesp –Vcalc), onde o Vesp é 0,97, ou seja, valor especificado para a barra 5 e o Vcalc é a tensão calculada na iteração, dessa forma os valores da matriz Ybarra que dependem do tap são recalculados e uma nova iteração é iniciada.
 deltat = Vo(5)-abs(V(5)); % Calculo da variação do TAP proporcional a variação de tensão na barra
 t = t + deltat; % Valor do TAP do Transformador 
 
�
MODELAGEM NO POWERWORLD:
Figura 6 – Modelagem do Sistema no PowerWorld
Resultados:
Tabela 7 – Resultados POWERWORLD
	Bus
	 
	 
	 
	 
	 
	 
	 
	 
	Number
	Nom kV
	PU Volt
	Volt (kV)
	Angle (Deg)
	Load MW
	Load Mvar
	Gen MW
	Gen Mvar
	1
	230
	1
	230
	0
	50
	30,99
	186,98
	122,95
	2
	230
	0,98242
	225,957
	-0,98
	170
	105,35
	 
	 
	3
	230
	0,96583
	222,141
	-1,85
	 
	 
	 
	 
	4
	230
	1,02
	234,6
	1,52
	80
	49,58
	318
	186,47
	5
	230
	0,97
	223,1
	-4,37
	200
	123,94
	 
	 
�
RESULTADOS:
Tabela 8 – Resultados Matlab vs Resultado Power Wolrd
	SOFTWARE
	TENSÕES FINAIS NAS BARRAS (PU)
	
	V2
	V3
	V4
	V5
	TAP
	MATLAB
	0.982421
/_-0.975774°
	0.966053
/_-1.846793°
	1.020000
/_1.523655°
	0.970000
/_-4.218671°
		1.031422
	POWERWOLRD
	0,98242/_-0,98°
	0,96583/_--1,85°
	1,02/_1,52°
	0,97/_-4,37°
	1,03359
Comentários: No Powerworld conforme pode ser visto na figura 6 foi inserido um transformador entre as barras 3 e 5 com reatância de 0,02 pu. Neste transformador foi habilitado o Controlador Automático a fim de manter a tensão na barra 5 em 0,97 pu. Os resultados da modelagem do powerwolrd estão sendo mostrados na tabela 7 e comparados com o do matlab na tabela 8. Conforme pode ser observado os resultados estão bem próximos o que evidencia a correta implementação do método no matlab e modelagem do circuito no Powerworld.
Figura 7- Resultado Final do TAP no PowerWolrd 
�
QUESTÃO 6:
IMPLEMENTAÇÃO MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Programa de Pós Graduação em Engenharia Elétrica
% Disciplina Análise de Sistemas de Energia 2016.1
% Metódo Gauss-Seidel
% Questão 6
% Profº: Niraldo
% Alunos: 1 - André Duarte
% 2 - Daniel Vasconcelos
% 3 - Humberto Rocha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
clc;
clear;
 
tic;
nb = 5; %quantidade de barras
 
e = 1; %tolerância
 
% Bases adotadas
Sb = 100; % 100 MVA
Vb = 230; % 230 kV
t = 1.05; % Valor pre-especificado do TAP
Yt = 1/(0.02*i);
 
% Montando a matriz admitância de barra
Yb = zeros(nb,nb);
Yb(1,1) = 8.985190 - 44.835953*i;
Yb(2,1) = -3.815629 + 19.078144*i;
Yb(3,1) = -5.169561 + 25.847809*i;
Yb(1,2) = Yb(2,1);
Yb(2,2) = Yb(1,1);
Yb(4,2) = Yb(3,1);
Yb(1,3) = Yb(4,2);
Yb(3,3) = (8.193267 - 40.863838*i)+(t^2)*Yt;
Yb(4,3) = -3.023705 + 15.118528*i;
Yb(2,4) = Yb(3,1);
Yb(3,4) = Yb(4,3);
Yb(4,4) = (8.193267 - 40.863838*i);
Yb(3,5) = -t*Yt;
Yb(5,3) = -t*Yt;
Yb(5,5) = Yt;
 
% Vetor de tensões iniciais nas barras
Vo = [1; 1; 1; 1.02; 1];
 
% Vetor de Potência Ativa e Reativa
Sc = zeros(nb,1);
Sg = zeros(nb,1);
 
% Barra 1 é a de referência
Sc(1) = (-50-30.99*i)/Sb;
% Barra 2 e 3 de carga
Sc(2) = (-170-105.35*i)/Sb;
Sc(5) = (-200-123.94*i)/Sb;
% Barra 4 é a de tensão controlada
Sc(4) = (-80-49.58*i)/Sb;
Sg(4) = 318/Sb;
 
P = real(Sc + Sg);
Q = imag(Sc + Sg);
 
% Metódo Gauss-Seidel
V = Vo;
k= 0;
Sg1o = zeros(nb,1);
while (e >= 10e-5) 
 k = k+1;
 Icalc = zeros(nb,1);
 for m=2:nb
 soma = 0;
 for j=1:nb;
 if (m ~= j) 
 soma = soma + Yb(m,j)*V(j);
 end
 Icalc(m) = Yb(m,j)*V(j)+Icalc(m);
 end 
 if (m == 4) % Calculo de Q e V para a barra PV
 Qcalc(k) = -imag(conj(V(m))*Icalc(m));
 V(m) = ((P(m)-i*Qcalc(k))/conj(V(m)) - soma)/ Yb(m,m);
 V(m) = abs(Vo(m))*V(m)/(abs(V(m))); % Corrigindo o módulo da tensão da barra PV
 else
 V(m) = ((P(m)-i*Q(m))/conj(V(m)) - soma)/ Yb(m,m); 
 end 
 end
 
 % Calculando o valor da corrente e potência gerada na Barra de Referência - Barra 1
 Icalc(1) = 0;
 for j=1:nb;
 Icalc(1) = Yb(1,j)*V(j)+Icalc(1);
 end
 Sg(1) = V(1)*conj(Icalc(1))-Sc(1);
 
 %Potência gerada na Barra 4
 Sg(4) = V(4)*conj(Icalc(4))-Sc(4); 
 
 Sg = [Sg(1); 0; 0; Sg(4); 0]; %Vetor com potência gerada em cada barra
 er = max(abs(real((Sg - Sg1o)))); % Cálculo do erro entre as Potências Calculadas na iteração
 ei = max(abs(imag((Sg - Sg1o))));
 e = max(er,ei);
 Sg1o = Sg;
end
 
% Cálculo dos Fluxos de Potência nas Linhas
Yshunt = zeros(nb);
Yshunt(1,2) = 0.05125*i;
Yshunt(2,1) = Yshunt(1,2);
Yshunt(1,3) = 0.03875*i;
Yshunt(3,1) = Yshunt(1,3);
Yshunt(2,4) = Yshunt(1,3);
Yshunt(4,2) = Yshunt(2,4);
Yshunt(3,4) = 0.06375*i;
Yshunt(4,3) = Yshunt(3,4);
for m=1:nb
 for j=1:nb
 if (m ~= j)
 SL(j,m) = V(j)*(conj((V(j)-V(m))*(-Yb(j,m)))); % Potencia na linha
 S(j,m) = SL(j,m)- abs(V(j))*(abs(V(j))*Yshunt(j,m));% Potencia na saida e entrada da barra
 I(j,m) = (S(j,m)/V(j))*Sb*10^3/(sqrt(3)*Vb); % Corrente na saída e entrada da barra
 end
 end
end
 
Iload5to3 = (-t*Yt*V(3,1)+Yt*V(5,1))*(Sb*10^3/(sqrt(3)*Vb)); % Corrente no secundário do transformador.
Iload3to5 = (t*t*Yt*V(3,1)-t*Yt*V(5,1))*(Sb*10^3/(sqrt(3)*Vb)); % Corrente no primário do transformador.
 
Perdas = sum((Sg+Sc))*Sb;
SL = SL*Sb; % Potência nas Linhas em MVA
S = S*Sb;
Sg = Sg*Sb; % Potência gerada por barra em MVA
Stg = sum(Sg)*Sb; % Potência total gerada no sistema em MVA
Stc = sum(Sc)*Sb; % Potência total das cargas do sistema em MVA
t1=toc;
 
fprintf(1,' \n'); 
fprintf(1,'**********************************************\n');
fprintf(1,'*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***\n');
fprintf(1,'**********************************************\n');
fprintf(1,' \n'); 
fprintf(1,'Quantidade de iterações: %d \n\n',k); 
fprintf(1,'Tensões Finais nas barras em PU:\n V1 = %f/_%f°\n V2 = %f/_%f°\n V3 = %f/_%f°\n V4 = %f/_%f°\n V5 = %f/_%f°\n\n',abs(V(1)),rad2deg(angle(V(1))),abs(V(2)),rad2deg(angle(V(2))),abs(V(3)),rad2deg(angle(V(3))),abs(V(4)),rad2deg(angle(V(4))),abs(V(5)),rad2deg(angle(V(5)))); 
fprintf(1,'Correntes entrando e saindo em cada Barra em A:\n\n I12 = %f/_%f° I21 = %f/_%f°\n I13 = %f/_%f° I31 = %f/_%f°\n I24 = %f/_%f° I42 = %f/_%f°\n I34 = %f/_%f° I43 = %f/_%f°\n I35 = %f/_%f° I53 = %f/_%f°\n\n',abs(I(1,2)),rad2deg(angle(I(1,2))),abs(I(2,1)),rad2deg(angle(I(2,1))),abs(I(1,3)),rad2deg(angle(I(1,3))),abs(I(3,1)),rad2deg(angle(I(3,1))),abs(I(2,4)),rad2deg(angle(I(2,4))),abs(I(4,2)),rad2deg(angle(I(4,2))),abs(I(3,4)),rad2deg(angle(I(3,4))),abs(I(4,3)),rad2deg(angle(I(4,3))), abs(Iload3to5),rad2deg(angle(Iload3to5)),abs(Iload5to3),rad2deg(angle(Iload5to3))); 
fprintf(1,'Potência Ativa Gerada nas barras em MW:\n P1 = %f\n P2 = %f\n P3 = %f\n P4 = %f\n P5 = %f\n\n',real(Sg(1)),real(Sg(2)),real(Sg(3)),real(Sg(4)),real(Sg(5)));
fprintf(1,'Potência Reativa Gerada nas barras em MVAr:\n Q1 = %f\n Q2 = %f\n Q3 = %f\n Q4 = %f\n Q5 = %f\n\n',imag(Sg(1)),imag(Sg(2)),imag(Sg(3)),imag(Sg(4)),imag(Sg(5)));
fprintf(1,'TAP FINAL: %f\n\n', t);
fprintf(1,'Tempo total de execução: %f s\n', t1);
RESULTADO:
**********************************************
*** RESULTADO FINAL - MÉTODO GAUSS SEIDEL ***
**********************************************
 
Quantidade de iterações: 24 
Tensões Finais nas barras em PU:
 V1 = 1.000000/_0.000000°
 V2 = 0.982421/_-0.975666°
 V3 = 0.966148/_-1.847507°
 V4 = 1.020000/_1.523758°
 V5 = 0.988574/_-4.133399°
Correntes entrando e saindo em cada Barra em A:
 I12 = 112.070072/_29.968162° I21 = 126.584089/_-139.931326°
 I13 = 301.156805/_35.004139° I31 = 312.278310/_-142.089221°
 I24 = 385.790474/_-149.628104° I42 = 376.230573/_27.820318°
 I34 = 315.392019/_-146.050813° I43 = 298.870752/_28.883412°
 I35 = 627.326473/_-35.919792° I53 = 597.453784/_144.080208°
Potência Ativa Gerada nas barras em MW:
 P1 = 186.947193
 P2 = 0.000000
 P3 = 0.000000
 P4 = 318.006976
 P5 = 0.000000
Potência Reativa Gerada nas barras em MVAr:
 Q1 = 122.111692
 Q2 = 0.000000
 Q3 = 0.000000
 Q4 = 185.964060
 Q5 = 0.000000
TAP FINAL: 1.050000
Tempo total de execução: 0.001927 s
COMENTÁRIOS: A diferença entre esta questão e a questão 5 é que nessa o que é pre-especificado agora é o valor do TAP, ou seja, de t que é de 1,05 e o valor de |V5| é estimado incialmente em 1 pu. Nesse caso não é necessário em cada iteração calcular o valor do TAP e sim somente o de V5, também não será necessário em cada iteração calcular o valor de Ybarra, o que torna consequentemente este problema mais simples.
�
MODELAGEM NO POWERWORLD:
Figura 8 – Modelagem do Sistema no PowerWorld
Resultados:
Tabela 9 – Resultados POWERWORLD
	Bus
	 
	 
	 
	 
	 
	 
	 
	 
	Number
	Nom kV
	PU Volt
	Volt (kV)
	Angle (Deg)
	Load MW
	Load Mvar
	Gen MW
	Gen Mvar
	1
	230
	1
	230
	0
	50
	30,99
	186,98
	122,95
	2
	230
	0,98242
	225,957
	-0,98
	170
	105,35
	 
	 
	3
	230
	0,96583
	222,141
	-1,85
	 
	 
	 
	 
	4
	230
	1,02
	234,6
	1,52
	80
	49,58
	318
	186,47
	5
	230
	0,9854
	226,642
	-4,37
	200
	123,94
	 
	 
�
RESULTADOS:
Tabela 10 – Resultados Matlab vs Resultado Power Wolrd
	SOFTWARE
	TENSÕES FINAIS NAS BARRAS (PU)
	
	V2
	V3
	V4
	V5
	TAP
	MATLAB
	0.982421
/_-0.975666°
	0.966148
/_-1.847507°
	1.020000
/_1.523758°
	0.988574
/_-4.133399°
		1.05
	POWERWOLRD
	0,98242/_-0,98°
	0,96583/_--1,85°
	1,02/_1,52°
	0,9854/_-4,37°
	1,05
Comentários: No Powerworld conforme pode ser visto na figura 6 foi inserido um transformador entre as barras 3 e 5 com reatância de 0,02 pu. Neste transformador não foi habilitado o Controlador Automático a fim de manter o valor do TAP fixo em 1,05. Os resultados da modelagem do powerwolrd estão sendo mostrados na tabela 9 e comparados com o do matlab na tabela 10. Conforme pode ser observado os resultados estão bem próximos o que evidencia a correta implementação do método no matlab e modelagem do circuito no Powerworld.
�UNIVERSIDADE FEDERAL DA BAHIA – UFBA
PROGRAMA DE PÓS GRADUAÇÃO EM ENGENHARIA ELÉTRICA
ENGA80-ANÁLISE DE SISTEMAS DE ENERGIA ELÉTRICA 2016.1
ANDRÉ JUAREZ
DANIEL COELHO VASCONCELOS
HUMBERTO ROCHA
_1533316918.unknown
_1533316922.unknown
_1533316924.unknown
_1533316925.unknown
_1533316926.unknown
_1533316923.unknown
_1533316920.unknown
_1533316921.unknown
_1533316919.unknown
_1533316916.unknown
_1533316917.unknown
_1533316915.unknown