Metodo Newton Rapshon no Matlab para um sistema de Quatro barras

Prévia do material em texto

%-------Universidade Federal de Sergipe ------------------- 
 
%-------Disciplina: Sistemas Elétricos de Potência--------- 
 
 
Elabore uma rotina em Matlab para resolução do fluxo de carga pelo método de Newton-Raphson para o 
sistema abaixo. 
 
 
Aluno: 
%-------Dados de entrada do Sistema elétrico com Quatro barras----------------------- 
 
YBarra = [3-12i, (-2+8i), (-1+4i), 0;(-2+8i), (3.666-14.64i),-0.666+2.64i,-1+4i;-1+4i,- 0.666+2.64i,3.666-
14.64i,-2+8i;0 ,-1+4i , -2+8i,3-12i]; 
G=[3,-2,-1,0;-2,3.666,-0.666,-1;-1,-0.666,3.666,-2;0,-1,-2,3]; %Matriz condutância 
B=[-12,8,4,0;8,-14.64,2.64,4;4,2.64,-14.64,8;0,4,8,-12]; %Matriz Suspectância 
nBarras=4; % número de barras 
V=[1.06;1;1;1]; % Tensão nas barras 
W=[0;0;0;0]; % fase das Tensões nas barras 
Pprog=[-0.5;-0.4;-0.3]; % A ordem é 2, 3, 4 em relação as barras 
Qprog=[-0.2;-0.3;-0.1]; % A ordem é 2, 3, 4 em relação as barras 
MagYbar=[12.37, 8.25, 4.123, 0; 8.25, 15.09, 2.72, 4.123; 4.123, 2.72, 15.09, 8.25; 0, 4.123, 8.25, 12.37]; 
FasYbar=(pi/180)*[-75.9, 104, 104, 0; 104, -75.9, 104, 104;104, 104, -75.9, 104; 0, 104, 104, -75.9]; 
 
% Teste de Convergência 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*cos(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Pcal(L) = (V(L)^2)*G(L,L) + Acumulator; 
end 
DPQ(1)= Pprog(1)-Pcal(2); % valores de potência que serão adicionados 
DPQ(2)= Pprog(2)-Pcal(3); % valores de potência que serão adicionados 
DPQ(3)= Pprog(3)-Pcal(4); % valores de potência que serão adicionados 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*sin(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Qcal(L) = (-1)*(V(L)^2)*B(L,L) - Acumulator; 
end 
DPQ(4)= Qprog(1)-Qcal(2); % valores de potência que serão adicionados 
DPQ(5)= Qprog(2)-Qcal(3); % valores de potência que serão adicionados 
DPQ(6)= Qprog(3)-Qcal(4); % valores de potência que serão adicionados 
fprintf('\n\niteração V1 th1 V2 th2 V3 th3 
V4 th4 \n'); 
fprintf('%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... 
 0,V(1),W(1),V(2),W(2),V(3),W(3),V(4),W(4)); 
 
 
 
 
%Primeira Iteração 
%|-----------| 
%| M | LH | 
%|-----------| 
%| N | UH | 
%|-----------| 
%Calculando M a matriz que contém as derivadas parciais da Potência ativa em relação as fases 
for i=2:4 
 for j=1:4 
 if i~=j 
 M(i,j)= -(abs(V(j)*V(i))*MagYbar(i,j))*sin(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
M(1,1)=0;M(2,2)=0; M(3,3)=0; M(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 M(k,k)= -M(k,l)+M(k,k); 
 end 
 end 
end 
 
 
%Calculando N, a matriz que contém as derivadas parciais da Potência reativa em relação as Fases 
for i=2:4 
 for j=1:4 
 if i~=j 
 N(i,j)= -abs((V(i)*V(j))*MagYbar(i,j))*cos(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
N(1,1)=0; N(2,2)=0; N(3,3)=0; N(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 N(k,k)= -N(k,l)+N(k,k); 
 end 
 end 
end 
 
%Calculando LH matriz que contém as derivadas parciais da Potência ativa em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 LH(i,j)=-N(i,j); 
 else 
 LH(i,i) = N(i,i)+2*(V(i)^2)*G(i,i); 
 end 
 end 
end 
 
%Calculando UH matriz que contém as derivadas parciais da Potência reativa em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 UH(i,j)= M(i,j); 
 else 
 UH(i,i)= -M(i,i)-2*(V(i)^2)*B(i,i); 
 end 
 end 
end 
%Obtenção da Matriz do Jacobiano 
 
for a=1:6 
 for d=1:6 
 if (d<=3 && a<=3) 
 Jacob(d,a)= M(a+1,d+1); 
 end 
 end 
end 
 
for a2=1:3 
 for d2=1:3 
 Jacob(d2+3,a2)= N(a2+1,d2+1); 
 end 
end 
 
Jacob(1,4)= LH(2,2); 
Jacob(2,4)= LH(3,2); 
Jacob(3,4)= LH(4,2); 
Jacob(1,5)= LH(2,3); 
Jacob(2,5)= LH(3,3); 
Jacob(3,5)= LH(4,3); 
Jacob(1,6)= LH(2,4); 
Jacob(2,6)= LH(3,4); 
Jacob(3,6)= LH(4,4); 
 
 
for a4=1:3 
 for d4=1:3 
 Jacob(d4+3,a4+3)= UH(a4+1,d4+1); 
 end 
end 
Jacob2=inv(Jacob); % Matriz inversa do Jacobiano 
TV=Jacob2*DPQ'; % Vetor de 6 posições 3 fases e 3 tensões fase 2, 3 e 4 e tensões 2, 3 e 4 
W(2)=W(2)+TV(1); %atualização dos valores das fases 
W(3)=W(3)+TV(2); %atualização dos valores das fases 
W(4)=W(4)+TV(3); %atualização dos valores das fases 
 
V(2)=V(2)*( 1 + TV(4)); 
V(3)=V(4)*( 1 + TV(5)); 
V(4)=V(4)*( 1 + TV(6)); 
 
 
% Teste de Convergência para Primeira Iteração 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*cos(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Pcal(L) = (V(L)^2)*G(L,L) + Acumulator; 
end 
DPQ(1)= Pprog(1)-Pcal(2); % valores de potência para serem adicionados 
DPQ(2)= Pprog(2)-Pcal(3); % valores de potência para serem adicionados 
DPQ(3)= Pprog(3)-Pcal(4); % valores de potência para serem adicionados 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*sin(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Qcal(L) = (-1)*(V(L)^2)*B(L,L) - Acumulator; 
end 
DPQ(4)= Qprog(1)-Qcal(2); % valores de potência para serem adicionados 
DPQ(5)= Qprog(2)-Qcal(3); % valores de potência para serem adicionados 
DPQ(6)= Qprog(3)-Qcal(4); % valores de potência para serem adicionados 
fprintf('%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... 
 1,V(1),W(1),V(2),W(2),V(3),W(3),V(4),W(4)); 
 
%Segunda Iteração 
%|-----------| 
%| M | LH | 
%|-----------| 
%| N | UH | 
%|-----------| 
 
%Calculando M a matriz que contém as derivadas parciais em relação as fases 
for i=2:4 
 for j=1:4 
 if i~=j 
 M(i,j)= -(abs(V(j)*V(i))*MagYbar(i,j))*sin(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
M(1,1)=0;M(2,2)=0; M(3,3)=0; M(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 M(k,k)= -M(k,l)+M(k,k); 
 end 
 end 
end 
 
%Calculando N matriz que contém as derivadas parciais em relação as Tensões 
for i=2:4 
 for j=1:4 
 if i~=j 
 N(i,j)= -abs((V(i)*V(j))*MagYbar(i,j))*cos(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
N(1,1)=0; N(2,2)=0; N(3,3)=0; N(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 N(k,k)= -N(k,l)+N(k,k); 
 end 
 end 
end 
 
%Calculando LH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 LH(i,j)=-N(i,j); 
 else 
 LH(i,i) = N(i,i)+2*(V(i)^2)*G(i,i); 
 end 
 end 
end 
 
%Calculando UH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 UH(i,j)= M(i,j); 
 else 
 UH(i,i)= -M(i,i)-2*(V(i)^2)*B(i,i); 
 end 
 end 
end 
 
%Obtenção da Matriz do Jacobiano 
 
for a=1:6 
 for d=1:6 
 if (d<=3 && a<=3) 
 Jacob(d,a)= M(a+1,d+1); 
 end 
 end 
end 
 
for a2=1:3 
 for d2=1:3 
 Jacob(d2+3,a2)= N(a2+1,d2+1); 
 end 
end 
 
Jacob(1,4)= LH(2,2); 
Jacob(2,4)= LH(3,2); 
Jacob(3,4)= LH(4,2); 
Jacob(1,5)= LH(2,3); 
Jacob(2,5)= LH(3,3); 
Jacob(3,5)= LH(4,3); 
Jacob(1,6)= LH(2,4); 
Jacob(2,6)= LH(3,4); 
Jacob(3,6)= LH(4,4); 
 
 
for a4=1:3 
 for d4=1:3 
 Jacob(d4+3,a4+3)= UH(a4+1,d4+1); 
 end 
end 
 
Jacob; 
Jacob2=inv(Jacob); % Até aqui está OK 
TV=Jacob2*DPQ'; % Vetor de 6 posições 3 fases e 3 tensões fase2,3 e 4 e tensões 2,3 e 4 
TV; 
W(2)=W(2)+TV(1); %atualização dos valores das fases 
W(3)=W(3)+TV(2); %atualização dos valores das fases 
W(4)=W(4)+TV(3); %atualização dos valores das fases 
 
V(2)=V(2)*( 1 + TV(4)); 
V(3)=V(4)*(1 + TV(5)); 
V(4)=V(4)*( 1 + TV(6)); 
 
% Teste de Convergência p/ Segunda Iteração 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*cos(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Pcal(L) = (V(L)^2)*G(L,L) + Acumulator; 
end 
DPQ(1)= Pprog(1)-Pcal(2); % valores de potência para serem adicionados 
DPQ(2)= Pprog(2)-Pcal(3); % valores de potência para serem adicionados 
DPQ(3)= Pprog(3)-Pcal(4); % valores de potência para serem adicionados 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*sin(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Qcal(L) = (-1)*(V(L)^2)*B(L,L) - Acumulator; 
end 
DPQ(4)= Qprog(1)-Qcal(2); % valores de potência para serem adicionados 
DPQ(5)= Qprog(2)-Qcal(3); % valores de potência para serem adicionados 
DPQ(6)= Qprog(3)-Qcal(4); % valores de potência para serem adicionados 
fprintf('%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... 
 2,V(1),W(1),V(2),W(2),V(3),W(3),V(4),W(4)); 
 
%Terceira Iteração 
%|-----------| 
%| M | LH | 
%|-----------| 
%| N | UH | 
%|-----------| 
%Calculando M a matriz que contém as derivadas parciais em relação as fases 
for i=2:4 
 for j=1:4 
 if i~=j 
 M(i,j)= -(abs(V(j)*V(i))*MagYbar(i,j))*sin(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
M(1,1)=0;M(2,2)=0; M(3,3)=0; M(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 M(k,k)= -M(k,l)+M(k,k); 
 end 
 end 
end 
 
%Calculando N matriz que contém as derivadas parciais em relação as Tensões 
for i=2:4 
 for j=1:4 
 if i~=j 
 N(i,j)= -abs((V(i)*V(j))*MagYbar(i,j))*cos(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
N(1,1)=0; N(2,2)=0; N(3,3)=0; N(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 N(k,k)= -N(k,l)+N(k,k); 
 end 
 end 
end 
 
%Calculando LH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 LH(i,j)=-N(i,j); 
 else 
 LH(i,i) = N(i,i)+2*(V(i)^2)*G(i,i); 
 end 
 end 
end 
 
%Calculando UH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 UH(i,j)= M(i,j); 
 else 
 UH(i,i)= -M(i,i)-2*(V(i)^2)*B(i,i); 
 end 
 end 
end 
%Obtenção da Matriz do Jacobiano 
 
for a=1:6 
 for d=1:6 
 if (d<=3 && a<=3) 
 Jacob(d,a)= M(a+1,d+1); 
 end 
 end 
end 
 
for a2=1:3 
 for d2=1:3 
 Jacob(d2+3,a2)= N(a2+1,d2+1); 
 end 
end 
 
Jacob(1,4)= LH(2,2); 
Jacob(2,4)= LH(3,2); 
Jacob(3,4)= LH(4,2); 
Jacob(1,5)= LH(2,3); 
Jacob(2,5)= LH(3,3); 
Jacob(3,5)= LH(4,3); 
Jacob(1,6)= LH(2,4); 
Jacob(2,6)= LH(3,4); 
Jacob(3,6)= LH(4,4); 
 
 
for a4=1:3 
 for d4=1:3 
 Jacob(d4+3,a4+3)= UH(a4+1,d4+1); 
 end 
end 
Jacob2=inv(Jacob); % Inversa do Jacobiano 
TV=Jacob2*DPQ'; % Vetor de 6 posições 3 fases e 3 tensões fase 2, 3 e 4 e tensões 2, 3 e 4 
TV; 
W(2)=W(2)+TV(1); %atualização dos valores das fases 
W(3)=W(3)+TV(2); %atualização dos valores das fases 
W(4)=W(4)+TV(3); %atualização dos valores das fases 
 
V(2)=V(2)*( 1 + TV(4)); 
V(3)=V(4)*( 1 + TV(5)); 
V(4)=V(4)*( 1 + TV(6)); 
 
 
% Teste de Convergência para Terceira Iteração 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*cos(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Pcal(L) = (V(L)^2)*G(L,L) + Acumulator; 
end 
DPQ(1)= Pprog(1)-Pcal(2); % valores de potência para serem adicionados 
DPQ(2)= Pprog(2)-Pcal(3); % valores de potência para serem adicionados 
DPQ(3)= Pprog(3)-Pcal(4); % valores de potência para serem adicionados 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*sin(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Qcal(L) = (-1)*(V(L)^2)*B(L,L) - Acumulator; 
end 
DPQ(4)= Qprog(1)-Qcal(2); % valores de potência para serem adicionados 
DPQ(5)= Qprog(2)-Qcal(3); % valores de potência para serem adicionados 
DPQ(6)= Qprog(3)-Qcal(4); % valores de potência para serem adicionados 
fprintf('%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... 
 3,V(1),W(1),V(2),W(2),V(3),W(3),V(4),W(4)); 
 
%Quarta Iteração 
%|-----------| 
%| M | LH | 
%|-----------| 
%| N | UH | 
%|-----------| 
 
%Calculando M a matriz que contém as derivadas parciais em relação as fases 
for i=2:4 
 for j=1:4 
 if i~=j 
 M(i,j)= -(abs(V(j)*V(i))*MagYbar(i,j))*sin(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
M(1,1)=0;M(2,2)=0; M(3,3)=0; M(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 M(k,k)= -M(k,l)+M(k,k); 
 end 
 end 
end 
 
%Calculando N matriz que contém as derivadas parciais em relação as Tensões 
for i=2:4 
 for j=1:4 
 if i~=j 
 N(i,j)= -abs((V(i)*V(j))*MagYbar(i,j))*cos(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
N(1,1)=0; N(2,2)=0; N(3,3)=0; N(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 N(k,k)= -N(k,l)+N(k,k); 
 end 
 end 
end 
 
%Calculando LH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 LH(i,j)=-N(i,j); 
 else 
 LH(i,i) = N(i,i)+2*(V(i)^2)*G(i,i); 
 end 
 end 
end 
 
%Calculando UH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 UH(i,j)= M(i,j); 
 else 
 UH(i,i)= -M(i,i)-2*(V(i)^2)*B(i,i); 
 end 
 end 
end 
%Obtenção da Matriz do Jacobiano 
 
for a=1:6 
 for d=1:6 
 if (d<=3 && a<=3) 
 Jacob(d,a)= M(a+1,d+1); 
 end 
 end 
end 
 
for a2=1:3 
 for d2=1:3 
 Jacob(d2+3,a2)= N(a2+1,d2+1); 
 end 
end 
 
Jacob(1,4)= LH(2,2); 
Jacob(2,4)= LH(3,2); 
Jacob(3,4)= LH(4,2); 
Jacob(1,5)= LH(2,3); 
Jacob(2,5)= LH(3,3); 
Jacob(3,5)= LH(4,3); 
Jacob(1,6)= LH(2,4); 
Jacob(2,6)= LH(3,4); 
Jacob(3,6)= LH(4,4); 
 
 
for a4=1:3 
 for d4=1:3 
 Jacob(d4+3,a4+3)= UH(a4+1,d4+1); 
 end 
end 
 
Jacob; 
Jacob2=inv(Jacob); % Até aqui está OK 
TV=Jacob2*DPQ'; % Vetor de 6 posições 3 fases e 3 tensões fase2,3 e 4 e tensões 2,3 e 4 
TV; 
W(2)=W(2)+TV(1); %atualização dos valores das fases 
W(3)=W(3)+TV(2); %atualização dos valores das fases 
W(4)=W(4)+TV(3); %atualização dos valores das fases 
 
V(2)=V(2)*( 1 + TV(4)); 
V(3)=V(4)*( 1 + TV(5)); 
V(4)=V(4)*( 1 + TV(6)); 
 
% Teste de Convergência p/ Quarta Iteração 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*cos(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Pcal(L) = (V(L)^2)*G(L,L) + Acumulator; 
end 
DPQ(1)= Pprog(1)-Pcal(2); % valores de potência para serem adicionados 
DPQ(2)= Pprog(2)-Pcal(3); % valores de potência para serem adicionados 
DPQ(3)= Pprog(3)-Pcal(4); % valores de potência para serem adicionados 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*sin(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Qcal(L) = (-1)*(V(L)^2)*B(L,L) - Acumulator; 
end 
DPQ(4)= Qprog(1)-Qcal(2); % valores de potência para serem adicionados 
DPQ(5)= Qprog(2)-Qcal(3); % valores de potência para serem adicionados 
DPQ(6)= Qprog(3)-Qcal(4); % valores de potência para serem adicionados 
fprintf('%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... 
 4,V(1),W(1),V(2),W(2),V(3),W(3),V(4),W(4)); 
 
%Quinta Iteração 
%|-----------| 
%| M | LH | 
%|-----------| 
%| N | UH | 
%|-----------| 
 
%Calculando M a matriz que contém as derivadas parciais em relação as fases 
for i=2:4 
 for j=1:4 
 if i~=j 
 M(i,j)=-(abs(V(j)*V(i))*MagYbar(i,j))*sin(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
M(1,1)=0;M(2,2)=0; M(3,3)=0; M(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 M(k,k)= -M(k,l)+M(k,k); 
 end 
 end 
end 
 
%Calculando N matriz que contém as derivadas parciais em relação as Tensões 
for i=2:4 
 for j=1:4 
 if i~=j 
 N(i,j)= -abs((V(i)*V(j))*MagYbar(i,j))*cos(FasYbar(i,j)+W(j)-W(i)); 
 end 
 end 
end 
N(1,1)=0; N(2,2)=0; N(3,3)=0; N(4,4)=0; 
for k=2:4 
 for l=1:4 
 if l~=k 
 N(k,k)= -N(k,l)+N(k,k); 
 end 
 end 
end 
 
%Calculando LH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 LH(i,j)=-N(i,j); 
 else 
 LH(i,i) = N(i,i)+2*(V(i)^2)*G(i,i); 
 end 
 end 
end 
 
%Calculando UH matriz que contém as derivadas parciais em relação as Tensões 
for i= 2:4 
 for j=2:4 
 if i~=j 
 UH(i,j)= M(i,j); 
 else 
 UH(i,i)= -M(i,i)-2*(V(i)^2)*B(i,i); 
 end 
 end 
end 
%Obtenção da Matriz do Jacobiano 
 
for a=1:6 
 for d=1:6 
 if (d<=3 && a<=3) 
 Jacob(d,a)= M(a+1,d+1); 
 end 
 end 
end 
 
for a2=1:3 
 for d2=1:3 
 Jacob(d2+3,a2)= N(a2+1,d2+1); 
 end 
end 
 
Jacob(1,4)= LH(2,2); 
Jacob(2,4)= LH(3,2); 
Jacob(3,4)= LH(4,2); 
Jacob(1,5)= LH(2,3); 
Jacob(2,5)= LH(3,3); 
Jacob(3,5)= LH(4,3); 
Jacob(1,6)= LH(2,4); 
Jacob(2,6)= LH(3,4); 
Jacob(3,6)= LH(4,4); 
 
 
for a4=1:3 
 for d4=1:3 
 Jacob(d4+3,a4+3)= UH(a4+1,d4+1); 
 end 
end 
 
Jacob; 
Jacob2=inv(Jacob); % Até aqui está OK 
TV=Jacob2*DPQ'; % Vetor de 6 posições 3 fases e 3 tensões fase2,3 e 4 e tensões 2,3 e 4 
TV; 
W(2)=W(2)+TV(1); %atualização dos valores das fases 
W(3)=W(3)+TV(2); %atualização dos valores das fases 
W(4)=W(4)+TV(3); %atualização dos valores das fases 
 
V(2)=V(2)*( 1 + TV(4)); 
V(3)=V(4)*( 1 + TV(5)); 
V(4)=V(4)*( 1 + TV(6)); 
 
% Teste de Convergência p/ Quinta Iteração 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*cos(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Pcal(L) = (V(L)^2)*G(L,L) + Acumulator; 
end 
DPQ(1)= Pprog(1)-Pcal(2); % valores de potência para serem adicionados 
DPQ(2)= Pprog(2)-Pcal(3); % valores de potência para serem adicionados 
DPQ(3)= Pprog(3)-Pcal(4); % valores de potência para serem adicionados 
 
for L= 2:4 
Acumulator=0; 
for C = 1:4 
 if C~=L 
 Acumulator = abs(V(L)*V(C))* MagYbar(L,C)*sin(FasYbar(L,C)+W(C)-W(L))+ Acumulator; 
 end 
end 
Qcal(L) = (-1)*(V(L)^2)*B(L,L) - Acumulator; 
end 
DPQ(4)= Qprog(1)-Qcal(2); % valores de potência para serem adicionados 
DPQ(5)= Qprog(2)-Qcal(3); % valores de potência para serem adicionados 
DPQ(6)= Qprog(3)-Qcal(4); % valores de potência para serem adicionados 
 
%Finalmente os resultados das iterações 
 V; % Este Vetor Apresenta as tensões nas barras corrigidas 
 W; % Este Vetor Apresenta as fases das tensões nas barras corrigidas em graus 
fprintf('%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... 
 5,V(1),W(1),V(2),W(2),V(3),W(3),V(4),W(4)); 
 
%Potência Ativa e Reativa na Barra Oscilante 
 
P1=V(1)*( V(1)*(G(1,1)*cos(W(1))) + V(2)*( G(1,2)*cos(W(1)-W(2))+B(1,2)*sin(W(1)-W(2)) )+ 
V(3)*(G(1,3)*cos(W(1)-W(3))+B(1,3)*sin(W(1)-W(3)) ) ); 
Q1=V(1)*( V(1)*(-B(1,1)*cos(W(1))) + V(2)*( G(1,2)*sin(W(1)-W(2))-B(1,2)*cos(W(1)-W(2)) )+ 
V(3)*(G(1,3)*sin(W(1)-W(3))-B(1,3)*cos(W(1)-W(3)) ) ); 
 
 
 
 
P1 = 1.2776 
 
Q1 = 0.7586