fluxo de carga 14 barras

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

FLUXO DE CARGA 14 BARRAS/Banco_dados_14_barras.m%% Dados do sistema de 14 barras
%% Potencia base do sistema (MVA)
Sbase = 100;
%% Dados de barras
% 1 2 3 4 5 6 *7 8 *9 10 11 *12 *13 *14 15 16 *17 *18
% nb tipo Pc Qc Pg Qg Gsh Bsh area tensao angulo Vnom Vmax Vmin Qmax Qmin Pmax Pmin
% MW MVAr MW MVAr MVAr p.u. graus kV p.u. p.u. MVAr MVAr MW MW
barras = [
 1 1 0.0 0.0 0.0 0.0 0.0 0.0 1 1.0600 0.0000 345.0 1.1000 0.9000 9999 -9999 250.0 10.0
 2 2 0.0 12.7 18.3 0.0 0.0 0.0 1 1.0450 0.0000 345.0 1.1000 0.9000 50.0 -40.0 250.0 10.0
 3 2 94.2 19.0 0.0 0.0 0.0 0.0 1 1.0100 0.0000 345.0 1.1000 0.9000 40.0 0.0 270.0 10.0
 4 3 47.8 -3.9 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 5 3 7.6 1.6 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 6 2 11.2 7.5 0.0 0.0 0.0 0.0 1 1.0700 0.0000 345.0 1.1000 0.9000 24.0 -6.0 270.0 10.0
 7 3 0.0 0.0 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 8 2 0.0 0.0 0.0 0.0 0.0 0.0 1 1.0900 0.0000 345.0 1.1000 0.9000 24.0 -6.0 270.0 10.0
 9 3 29.5 16.6 0.0 0.0 0.0 19.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 10 3 9.0 5.8 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 11 3 3.5 1.8 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 12 3 6.1 1.6 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 13 3 13.5 5.8 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0
 14 3 14.9 5.0 0.0 0.0 0.0 0.0 1 1.0000 0.0000 345.0 1.1000 0.9000 00.0 0.0 000.0 0.0];
% tipo: 1 - slack, 2 - PV, 3 - PQ
%% Dados das ramos
% 1 2 3 4 5 6
% de para r x b tap
% p.u. p.u. p.u. 
linhas = [
 1 2 0.03876 0.11834 0.02640 1.000
 1 2 0.03876 0.11834 0.02640 1.000
 1 5 0.05403 0.22304 0.04920 1.000
 2 3 0.04699 0.19797 0.04380 1.000
 2 4 0.05811 0.17632 0.03740 1.000
 2 5 0.05695 0.17388 0.03400 1.000
 3 4 0.06701 0.17103 0.03460 1.000
 4 5 0.01335 0.04211 0.01280 1.000
 4 7 0.00000 0.20912 0.00000 0.950
 4 9 0.00000 0.55618 0.00000 0.961
 5 6 0.00000 0.25202 0.00000 1.000
 6 11 0.09498 0.19890 0.00000 1.000
 6 12 0.12291 0.25581 0.00000 1.000
 6 13 0.06615 0.13027 0.00000 1.000
 7 8 0.00000 0.17615 0.00000 1.000
 7 9 0.00000 0.11001 0.00000 1.000
 9 10 0.03181 0.08450 0.00000 1.000
 9 14 0.12711 0.27038 0.00000 1.000
 10 11 0.08205 0.19207 0.00000 1.000
 12 13 0.22092 0.19988 0.00000 1.000
 13 14 0.17093 0.34802 0.00000 1.000];
FLUXO DE CARGA 14 BARRAS/Programa_Fluxo_Potencia_14_Rev1.m% FLUXO DE POTÊNCIA - MÉTODO DE NEWTON-RAPHSON
%%
%Elaboração de um programa capaz de calcular o fluxo e potência,foi 
%utilizado como base teorica as notas de aulas da disciplina ASEE
clc
clear all
%%
%legenda de impressão
fprintf('\n________________________________________________________________\n');
fprintf('UNIVERSIDADE ESTADUAL PAULISTA');
fprintf('\nCAMPUS ILHA SOLTEIRA');
fprintf('\nPROGRAMA DE PÓS-GRADUÇÃO EM ENGENHARIA ELÉTRICA');
fprintf('\nDISCIPLINA: ANALISE DE SISTEMAS DE ENERGIA ELÉTRICA');
fprintf('\nDOCENTE: DILSON AMANCIO ALVES');
fprintf('\nDISCENTES: BASAKUAU NKOMI NKOSI JÚNIOR e VAGNER ANTONIO DE MORAES DA CRUZ');
fprintf('\nFLUXO DE POTÊNCIA - MÉTODO DE NEWTON-RAPHSON');
fprintf('\n________________________________________________________________');
%% Escolha do sistema
eval('Banco_dados_14_barras'); %busca banco de dados das barras
erro=0.0001; %especifica a tolerância do sistema
% Dados de barras
Barra = barras(:,1); % Busca o número de cada barra
Tipo = barras(:,2); % Busca o tipo de cada barra
Pc = barras(:,3)/Sbase; % Busca o valor da potência ativa consumida em cada barra
Qc = barras(:,4)/Sbase; % Busca o valor da potência reativa consumida em cada barra
Pg = barras(:,5)/Sbase; % Busca o valor da potência ativa gerada nas barras
Qg = barras(:,6)/Sbase; % Busca o valor da potência reativa gerada nas barras
bsh = barras(:,8)/Sbase; % Busca o valor do elemento shunt das barras
V = barras(:,10); % Busca o valor de tensões das barras
theta = barras(:,11)*(pi/180); % Busca o valor de ângulos em radianos 
Vmax = barras(:,13); % Busca o valor de tensão máxima
Vmin = barras(:,14); % Busca o valor de tensão mínima
Qmax = barras(:,15)/Sbase; % Busca o valor da potência reativa máxima que pode ser gerado nas barras PV
Qmin = barras(:,16)/Sbase; % Busca o valor da potência reativa mínima que pode ser gerado nas barras PV
Vesp=V; % Salva o valor específico das barras
%% Dados das linhas do sistema
sai = linhas(:,1); % busca a linha de partida
chega = linhas(:,2); % busca a linha de chegada
r = linhas(:,3); % Resistência da linha
x = linhas(:,4); % Reatância da linha
bshl = linhas(:,5)/2; % Sucespetância shunt k e m
tap = linhas(:,6); % Tap trafos
nb = length(Barra); % Quantidade de barras do sistema 
nl = length(sai); % Quantidade de linhas 
 
%% Para o calculo do fluxo deve-se montar a matriz admitância (Y) do sistema
% como Y=G+jB, será montado as matrizes de condutancia(G) e susceptancia (B) nodais separadamente, isso também 
% facilitará os cálculos.
 G = zeros(nb,nb); % Matriz condutância 
 B = zeros(nb,nb); % Matriz suceptância 
 
for k = 1:nb;
 B(k,k) = bsh(k);% Elementos das diagonais principal
end 
% Matriz condutância e susceptância
for l = 1:nl;
 k = sai(l);
 m = chega(l);
 %Matriz condutância
 gkm = r(l)/(r(l)^2+x(l)^2); %calculo da condutância da linha
 	G(k,k) = G(k,k)+tap(l)^2*gkm; % Diagonal principal de G
 G(m,m) = G(m,m)+gkm;
 G(k,m) = G(k,m)-tap(l)*gkm; % Elementos fora da diagonal principal
 G(m,k) = G(m,k)-tap(l)*gkm;
 % Matriz suceptância 
 bkm = -x(l)/(r(l)^2+x(l)^2); %calculo da susceptância da linha
	B(k,k) = B(k,k)+tap(l)^2*bkm + bshl(l); % Elementos da diagonal principal de B
 	B(m,m) = B(m,m)+bkm+bshl(l);
 B(k,m) = B(k,m)-tap(l)*bkm; %Elementos fora da diagonal principal
 B(m,k) = B(m,k)-tap(l)*bkm; 
end
%% MÉTODO DE NEWTON
%Para iniciarmos o calculo do fluxo vamos escolher V igual a 1 p.u. e theta
% igual a 0 para todas as barras PQ e theta igual a 0 para barras PV.
% Os tipos de barra estão definidos assim: 1-slack, 2-PV, 3-PQ, temos:
% iniciar as iterações
iteracao = 0;
for k = 1:nb;
 if Tipo(k)>= 2 ; % Barra PV ou PQ
 theta(k) = 0;
 if Tipo(k) == 3; % Barra PQ 
 V(k) = 1.0;
 end
 end
end
%Temos que encontrar os deltas de P e Q menor ou igual a tolerância. 
%Calculo das potências especificadas, para barras PQ e PV.
for i = 1:nb;
 Pesp(i) = Pg(i) - Pc(i); % Potência ativa especificada barras PV e PQ.
Qesp(i) = Qg(i) - Qc(i); % Potência reativa especificada barra PQ.
end
%%
% Calculo das potencias P e Q nas barras 
Pks = zeros(nb,1);
Qks = zeros(nb,1);
%calculo do elemento relativo a barra que esta sendo calculada onde m igual a k
for k = 1:nb
 Pks(k) = G(k,k)*V(k)^2;% representa potencia ativa calculada
 Qks(k) = -B(k,k)*V(k)^2;% representa potencia reativa calculada
end
%calculo dos elementos vizinhos onde m é diferente de k
for i = 1:nl;
 k = sai(i);
 m = chega(i);
 t_km = (theta(k) - theta(m));
 Pks(k) = Pks(k) + V(k)*V(m)*(G(k,m)*cos(t_km)+B(k,m)*sin(t_km));
 Qks(k) = Qks(k) + V(k)*V(m)*(G(k,m)*sin(t_km)-B(k,m)*cos(t_km));
 Pks(m) = Pks(m) + V(k)*V(m)*(G(k,m)*cos(-t_km)+B(k,m)*sin(-t_km));
 Qks(m) = Qks(m) + V(k)*V(m)*(G(k,m)*sin(-t_km)-B(k,m)*cos(-t_km));
end
% Assim calculados P e Q especificados e Pks podemos realizar o teste de
% convergencia 
deltaP = zeros(nb,1);
deltaQ = zeros(nb,1);
for k = 1:nb
 if Tipo(k) >= 2
 deltaP(k) = Pesp(k) - Pks(k);
 end
 if Tipo(k) == 3
 deltaQ(k) = Qesp(k) - Qks(k);
 end
end
%encontra o maior error de mismatch
maxdeltaP = max(abs(deltaP)); % Potência ativa
maxdeltaQ = max(abs(deltaQ)); % Potencia reativa
%imprime o erro máximo
fprintf('_____________________________\n');
fprintf('máximo erro primeira iteração\n');
fprintf(' deltaP deltaQ\n');
disp([maxdeltaP maxdeltaQ]);
%iniciar os valores de controle e limite de PQ e PV
ultrap=zeros(nb,1);
%entra no processo de convergencia
while((maxdeltaP > erro)||(maxdeltaQ > erro));
 
% Equacionamento para a construção da matriz Jacobiana 
 H = zeros(nb,nb); M=H; N=H; L=H;
 for k = 1:nb;
 H(k,k) = -Qks(k)-V(k)*V(k)*B(k,k);
 N(k,k) = (Pks(k)+V(k)*V(k)*G(k,k))/V(k);
 M(k,k) = Pks(k)-V(k)*V(k)*G(k,k);
 L(k,k) = (Qks(k)-V(k)*V(k)*B(k,k))/V(k);
 if Tipo(k) == 1;
 H(k,k) = 10e10;
 end
 if Tipo(k) ~= 3;
 L(k,k) = 10e10;
 end
 end
 
 for k = 1:nb;
 for m = 1:nb;
 if k ~=m; 
 t_km = theta(k) - theta(m);
 H(k,m) = V(k)*V(m)*(G(k,m)*sin(t_km)-B(k,m)*cos(t_km));
 H(m,k) = V(k)*V(m)*(-G(k,m)*sin(t_km)-B(k,m)*cos(t_km));
 N(k,m) = V(k)*(G(k,m)*cos(t_km)+B(k,m)*sin(t_km));
 N(m,k) = V(m)*(G(k,m)*cos(t_km)-B(k,m)*sin(t_km));
 M(k,m) = -V(k)*V(m)*(G(k,m)*cos(t_km)+B(k,m)*sin(t_km));
 M(m,k) = -V(k)*V(m)*(G(k,m)*cos(t_km)-B(k,m)*sin(t_km));
 L(k,m) = V(k)*(G(k,m)*sin(t_km)-B(k,m)*cos(t_km));
 L(m,k) = V(m)*(-G(k,m)*sin(t_km)-B(k,m)*cos(t_km));
 end
 end
 end
J = [H N; M L];% Montagem Matriz Jacobiana
deltaPQ = [deltaP; deltaQ];% Matriz delta de P e Q
deltaVT = inv(J)*deltaPQ; %calculo da matriz delta de V e theta 
%atualizando os valores de theta e V
 for k = 1:nb;
 if Tipo(k) >1;
 theta(k) = theta(k) + deltaVT(k);
 end
 if Tipo(k) == 3;
 V(k) = V(k) + deltaVT(k+nb);
 end
 end
% Atualizar os valores de estado P, Theta e Q
 for k = 1:nb;
 Pks(k) = G(k,k)*V(k)*V(k);
 Qks(k) = -B(k,k)*V(k)*V(k);
 end
 for k = 1:nb;
 for m = 1:nb;
 if k ~=m;
 t_km = theta(k) - theta(m);
 Pks(k) = Pks(k) + V(k)*V(m)*(G(k,m)*cos(t_km)+B(k,m)*sin(t_km));
 Qks(k) = Qks(k) + V(k)*V(m)*(G(k,m)*sin(t_km)-B(k,m)*cos(t_km));
 end
 end
 end
 for k = 1:nb;
 if Tipo(k) ~=1;
 deltaP(k) = Pesp(k) - Pks(k);
 end
 if Tipo(k) == 3;
 deltaQ(k) = Qesp(k) - Qks(k);
 end
 end
 maxdeltaP = max(abs(deltaP)); %error máximo de potência ativa
 maxdeltaQ = max(abs(deltaQ)); % error máximo potência reativa
 
iteracao=iteracao+1; %atualiza o número de iterações
%Limite de tensão da barra PQ 
 for k=1:nb ; 
 if (ultrap(k)==1) && (V(k)>Vesp(k)) ;% na primeira iteração o valor da ultrapassagem é O.
 Tipo(k)=2; % troca para barra PV se ultrapassar Vmáx
 V(k)=Vesp(k); % fixa V
 ultrap(k)=0; % zera a ultrapassagem
 
 elseif (ultrap(k)==2) && (V(k)<Vesp(k)); 
 Tipo(k)=2; % troca para barra PV se ultrapassar Vmin
 V(k)=Vesp(k); % fixa V
 ultrap(k)=0; % zera a ultrapassagem
 
 end
 end
 % Controle de tensão em barra PV 
 for k=1:nb;
 if (Tipo(k)==2) ; 
 if Qc(k)+Qks(k)>Qmax(k) && iteracao > 1 && ultrap(k)==0;
 Tipo(k)=3;
 Qesp(k)=Qmax(k)-Qc(k);
 ultrap(k)=1; 
 end
 if Qc(k)+Qks(k)<Qmin(k) && iteracao > 1 && ultrap(k)==0;
 Tipo(k)=3;
 Qesp(k)=Qmin(k)-Qc(k);
 ultrap(k)=2; 
 end
 end
 end
end
fprintf('_____________________________\n');
fprintf('máximo erro após a iteração\n');
fprintf(' deltaP deltaQ\n');
disp([maxdeltaP maxdeltaQ]);
%SUBSISTEMA 2
%Calculo das potencia ativa, reativas e perdas nas linhas
Pkm = zeros(nl,1);
Qkm = zeros(nl,1);
Pmk = zeros(nl,1);
Qmk = zeros(nl,1);
Pperdas = zeros(nl,1);
Qperdas = zeros(nl,1);
for i = 1:nl;
 k = sai(i);
 m = chega(i);
 t_km = (theta(k) - theta(m)); 
 gkm = r(i)/(r(i)^2+x(i)^2);
 bkm = -x(i)/(r(i)^2+x(i)^2);
 Vkm = V(k) * V(m);
 t = tap(i);
 Pkm(i) = (t^2*V(k)*V(k)*gkm - t*Vkm*(gkm*cos(t_km)+bkm*sin(t_km)))*Sbase;
 Qkm(i) = (-t^2*V(k)*V(k)*(bkm + bshl(i))-t*Vkm*(gkm*sin(t_km)-bkm*cos(t_km)))*Sbase;
 Pmk(i) = (V(m)*V(m)*gkm - t*Vkm*(gkm*cos(t_km)-bkm*sin(t_km)))*Sbase;
 Qmk(i) = (-V(m)*V(m)*(bkm + bshl(i))+t*Vkm*(gkm*sin(t_km)+bkm*cos(t_km)))*Sbase;
 Pperdas(i) = Pkm(i) + Pmk(i);
 Qperdas(i) = Qkm(i) + Qmk(i);
end
 
%Relatorio
fprintf('\n________________________________________________________________');
fprintf('\nO Processo iterativo convergiu na %2d iteração ', iteracao);
fprintf('\n \n\n');
%Resultado nas barras do sistema
fprintf('RESULTADO DAS BARRAS\n\n');
fprintf('N° V Theta(º) P Q \n');
fprintf('_____________________________________________________\n');
for i = 1:nb;
 fprintf('%3d \t% 5.4f \t% 5.4f \t% 5.4f \t% 5.4f \n',Barra(i), V(i), theta(i)*180/pi, Pks(i), Qks(i) );
end
%Resultado nas linhas do sistema
fprintf('\n\nRESULTADOS NAS LINHAS\n\n');
fprintf('Linhas Pkm Qkm Pmk Qmk Pperdas Qperdas\n');
fprintf('__________________________________________________________________________________\n');
for i = 1:nl
 fprintf('%1d - %1d \t% 5.4f \t% 5.4f \t% 5.4f \t% 5.4f \t% 5.4f \t% 5.4f\n', sai(i), chega(i), Pkm(i), Qkm(i), Pmk(i), Qmk(i), Pperdas(i), Qperdas(i));
 
end
TotalPp=sum(Pperdas(1:nl));
TotalQp=sum(Qperdas(1:nl));
fprintf('_________________________');
fprintf('\nTOTAL DE PERDAS \n');
fprintf('_________________________\n');
fprintf(' Pperdas Qperdas\n');
fprintf(' MVA MVAr\n');
fprintf('_________________________\n');
fprintf('\t% 5.4f \t% 5.4f\n', TotalPp, TotalQp);
fprintf('\nNa iteracao %3d a barra %3d violou o limite maximo de tensao \n',iteracao,k);
fprintf('\nNa iteracao %3d a barra %3d ultrapassou o limite maximo de reativo \n',iteracao,k);
fprintf('\nNa iteracao: %3d a barra %3d retorno para os limites \n', iteracao,k);