%-------Universidade Federal de Sergipe ------------------%-------Disciplina: %-------Disc iplina: Sistemas Elétricos de Potência-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.66614.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 fprintf( '\n\niteração V1 th1 V2 th2 V3 th3 V4 th4 \n' \n'); ); fprintf('%6i fprintf( '%6i %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',... %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)= Jacob(2,4)= Jacob(3,4)= Jacob(1,5)= Jacob(2,5)= Jacob(3,5)= Jacob(1,6)= Jacob(2,6)= Jacob(3,6)=
LH(2,2); LH(3,2); LH(4,2); LH(2,3); LH(3,3); LH(4,3); LH(2,4); LH(3,4); 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)= Jacob(2,4)= Jacob(3,4)= Jacob(1,5)= Jacob(2,5)= Jacob(3,5)= Jacob(1,6)= Jacob(2,6)= Jacob(3,6)=
LH(2,2); LH(3,2); LH(4,2); LH(2,3); LH(3,3); LH(4,3); LH(2,4); LH(3,4); 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); TV=Jacob2*DPQ'; % TV; W(2)=W(2)+TV(1); W(3)=W(3)+TV(2); W(4)=W(4)+TV(3);
% Até aqui está OK Vetor de 6 posições 3 fases e 3 tensões fase2,3 e 4 e tensões 2,3 e 4 %atualização dos valores das fases %atualização dos valores das fases %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)= Jacob(2,4)= Jacob(3,4)= Jacob(1,5)= Jacob(2,5)= Jacob(3,5)= Jacob(1,6)= Jacob(2,6)= Jacob(3,6)=
LH(2,2); LH(3,2); LH(4,2); LH(2,3); LH(3,3); LH(4,3); LH(2,4); LH(3,4); 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) end DPQ(1)= DPQ(2)= DPQ(3)=
= (V(L)^2)*G(L,L) + Acumulator; Pprog(1)-Pcal(2); Pprog(2)-Pcal(3); Pprog(3)-Pcal(4);
% valores de potência para serem adicionados % valores de potência para serem adicionados % 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)= Jacob(2,4)= Jacob(3,4)= Jacob(1,5)= Jacob(2,5)= Jacob(3,5)= Jacob(1,6)= Jacob(2,6)= Jacob(3,6)=
LH(2,2); LH(3,2); LH(4,2); LH(2,3); LH(3,3); LH(4,3); LH(2,4); LH(3,4); 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); TV=Jacob2*DPQ'; % TV; W(2)=W(2)+TV(1); W(3)=W(3)+TV(2); W(4)=W(4)+TV(3);
% Até aqui está OK Vetor de 6 posições 3 fases e 3 tensões fase2,3 e 4 e tensões 2,3 e 4 %atualização dos valores das fases %atualização dos valores das fases %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 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
%10.6f
%10.6f\n',...
for a2=1:3 for d2=1:3 Jacob(d2+3,a2)= N(a2+1,d2+1); end end Jacob(1,4)= Jacob(2,4)= Jacob(3,4)= Jacob(1,5)= Jacob(2,5)= Jacob(3,5)= Jacob(1,6)= Jacob(2,6)= Jacob(3,6)=
LH(2,2); LH(3,2); LH(4,2); LH(2,3); LH(3,3); LH(4,3); LH(2,4); LH(3,4); 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); TV=Jacob2*DPQ'; % TV; W(2)=W(2)+TV(1); W(3)=W(3)+TV(2); W(4)=W(4)+TV(3);
% Até aqui está OK Vetor de 6 posições 3 fases e 3 tensões fase2,3 e 4 e tensões 2,3 e 4 %atualização dos valores das fases %atualização dos valores das fases %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)*s in(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)*co s(W(1))) + V(2)*( G(1,2)*sin(W(1)- W(2))-B(1,2)*co s(W(1)-W(2)) )+ V(3)*(G(1,3)*s in(W(1)-W(3))-B (1,3)*cos(W(1)-W (3)) ) );
P1 = 1.2776 Q1 = 0.7586
