Buscar

Exemplo de Fluxo de Potência Newton Raphson (Matlab)

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

clc;
clear all;
%Tipo 1- Carga 2- Tensão 3-Referência
%Barras num_barra Tipo Tensão Fase(graus) PGerada QGerada PCarga QCarga cap_paralelo 
barra = [1 3 0.95 0.0 0.0 0.0 0.0 0.0 0.0
 2 2 0.95 0.0 .40 0.0 0.217 0.127 0.0
 3 2 0.96 0.0 0.0 0.0 0.942 0.190 0.0
 4 1 1.00 0.0 0.0 0.0 0.478 0.039 0.0
 5 1 1.00 0.0 0.0 0.0 0.076 0.016 0.0
 6 2 0.95 0.0 0.0 0.0 0.30 0.300 0.0
 7 3 1.00 0.0 0.0 0.0 0.0 0.0 0.0
 8 2 0.95 0.0 0.0 0.0 0.0 0.0 0.0
 9 1 1.00 0.0 0.0 0.0 0.295 0.166 0.0
 10 1 1.00 0.0 0.0 0.0 0.090 0.058 0.0
 11 1 1.00 0.0 0.0 0.0 0.035 0.018 0.0
 12 1 1.00 0.0 0.0 0.0 0.061 0.016 0.0
 13 1 1.00 0.0 0.0 0.0 0.135 0.058 0.0
 14 1 1.00 0.0 0.0 0.0 0.149 0.050 0.0];
 
%OrigemDestino Resistência Reatância capacitor_serie Fase gama (graus) 
trechos = [
 1 2 0.193 0.0592 0.059 0
 1 5 0.0540 0.02230 0.049 0
 2 3 0.04699 0.01980 0.044 0
 2 4 0.05811 0.1763 0.037 0 
 2 5 0.05695 0.1739 0.034 0 
 3 4 0.06701 0.1710 0.035 0 
 4 5 0.01335 0.0421 0.013 0
 4 7 0.000 0.2091 0.0 0
 4 9 0.000 0.5562 0.0 0
 5 6 0.000 0.2520 0.0 0
 6 11 0.09498 0.1989 0.000 0 
 6 12 0.12291 0.2558 0.000 0
 6 13 0.06615 0.1303 0.000 0 
 7 8 0.000 0.1761 0.000 0
 7 9 0.000 0.1100 0.000 0
 9 10 0.03181 0.0845 0.000 0
 9 14 0.12711 0.2704 0.000 0
10 11 0.08205 0.1921 0.000 0
12 13 0.22092 0.1999 0.000 0
13 14 0.17093 0.3480 0.000 0];
 
%	Determinação da dimensão das matrizes
[num_barras, colun1] = size(barra);
[num_trechos , colun2] = size(trechos);
%	Obtenção dos dados de barra
for k = 1:num_barras %Para resolver a letra c deve-se manualmente colocar a linha adicional
 cod_barra(k) = barra(k,1); %Codificação da linha
 tipo(k) = barra(k,2); %Indica tipo de barra: PQ, PV, ou Vteta
 V(k) = barra(k,3); %Tensão na barra em pu
 teta(k) = barra(k,4) * pi/180; %teta em radianos
 Pesp_g(k) = barra(k,5); %Potência ativa especificada gerada
 Qesp_g(k) = barra(k,6); %Potência reativa especificada gerada
 Pesp_c(k) = barra(k,7); %Potência ativa especificada carga
 Qesp_c(k) = barra(k,8); %Potência reativa especificada carga
 Cap_paralelo(k) = barra(k,9); %Capacitor em paralelo com a barra (basta adicionar valores para resolver a letra b)
 Pesp(k) = Pesp_g(k) - Pesp_c(k); %Potência ativa especificada (gerador - carga)
 Qesp(k) = Qesp_g(k) - Qesp_c(k); %Potência reativa especificada (gerador - carga)
 n_iteracoes(barra(k,1)) = k; %Número de iterações
end
%	Obtenção dos dados do trecho entre as barras
for a = 1:num_trechos
	barra_origem(a) = trechos(a,1); %Determina origem
	barra_destino(a) = trechos(a,2); %Determina destino
	r(a) = trechos(a,3); %Resitência
	x(a) = trechos(a,4); %Reatância
	bs(a) = trechos(a,5) / 2; %Montagem dos capacitores shunt no modelo pi
 gama(a) = trechos(a,6)*pi/180; %Ângulo de defasagem nas barras em radianos (dado em graus)
end
Y = zeros(num_barras,num_barras); %Cria um vetor de zeros para a admitancia com a dimensão do sistema
%Matriz de admitâncias
for k = 1:num_barras
	Y(k,k) = i*Cap_paralelo(k);
end
for a = 1:num_trechos %Montagem da matriz de admitâncias a partir dos dados fornecidos da linha e sistema
 k = barra_origem(a);
 m = barra_destino(a);
 y(a) = 1/(r(a) + i*x(a));
 Y(k,k) = Y(k,k) + y(a) + i*bs(a);
 Y(m,m) = Y(m,m) + y(a) + i*bs(a);
 Y(k,m) = Y(k,m) - y(a);
 Y(m,k) = Y(m,k) - y(a);
end
G = real(Y); %Parte real de Y - matriz de condutâncias
B = imag(Y); %Parte imaginária de Y - matriz de susceptâncias
for k=1:14
 for m=1:14
 modY(k,m)=(sqrt(G(k,m)^2+B(k,m)^2)); %Cálculo do módulo de Y
 end
end
for k=1:14
 for m=1:14
 if G(k,m)>0
 faseY(k,m)=(atan(B(k,m)/G(k,m))); %Cálculo da fase de Y
 else
 faseY(k,m)=0; %Prevenção contra indeterminação no caso de divisão por zero
 end
 end
end
v1p = ones(num_barras,1); %Cria uma matriz de uns com num_barras linhas e 1 coluna)
v1q = ones(num_barras,1); %Cria uma matriz de uns com num_barras linhas e 1 coluna)
v2p = zeros(num_barras,1); %Cria uma matriz de zeros com num_barras linhas e 1 coluna)
v2q = zeros(num_barras,1); %Cria uma matriz de zeros com num_barras linhas e 1 coluna)
v3p = zeros(num_barras,1); %Cria uma matriz de zeros com num_barras linhas e 1 coluna)
v3q = zeros(num_barras,1); %Cria uma matriz de zeros com num_barras linhas e 1 coluna)
deltaP = 10^10; %Limite máximo para a caondição
deltaQ = 10^10; %Limite máximo para a condição
ite=0; %Inicializa a variável que será incrementada a cada iteração
max_ite=2; %Número máximo de iterações
while deltaP > 1e-2 || deltaQ > 1e-2 || ite<2
 
 %Potências da carga calculada
 for k=1:num_barras 
 Pcalc_c(k) = G(k,k)*(V(k)^2);
 Qcalc_c(k) = -B(k,k)*(V(k)^2);
 end
 
 for a=1:num_trechos
 k = barra_origem(a);
 m = barra_destino(a);
 
 g = real(y(a)); %condutância do trecho km
 b = imag(y(a)); %susceptância do trecho km
 
 %Cálculo das potências que fluem entre as barras
 Pcalc_c(k) = Pcalc_c(k) + V(k)*V(m)*(-g*cos(teta(k) - teta(m) + gama(a))-b*sin(teta(k) - teta(m) + gama(a)));
 Pcalc_c(m) = Pcalc_c(m) + V(k)*V(m)*(-g*cos(teta(k) - teta(m) + gama(a))+b*sin(teta(k) - teta(m) + gama(a)));
 Qcalc_c(k) = Qcalc_c(k) + V(k)*V(m)*(-g*sin(teta(k) - teta(m) + gama(a))+b*cos(teta(k) - teta(m) + gama(a)));
 Qcalc_c(m) = Qcalc_c(m) - V(k)*V(m)*(-g*sin(teta(k) - teta(m) + gama(a))-b*cos(teta(k) - teta(m) + gama(a)));
 
 end
 
 DELTA_P = zeros(num_barras,1); %Cria uma matriz de zeros com o número de linhas do sistema e uma coluna
 DELTA_Q = zeros(num_barras,1); %Cria uma matriz de zeros com o número de linhas do sistema e uma coluna
 deltaP=0; %Coloca o a variável condicional igual a zero
 deltaQ=0; %Coloca o a variável condicional igual a zero
 contP=0; %Inicializa contador
 contQ=0; %Inicializa contador
 
 for k=1:num_barras
 Pprog(k)=(v1p(k)+v2p(k)*V(k)+v3p(k)*(V(k)^2))*Pesp(k); %dimensiona a matriz
 Qprog(k)=(v1q(k)+v2q(k)*V(k)+v3q(k)*(V(k)^2))*Qesp(k); %dimensiona a matriz
 if tipo(k)~=3 %desconsidera a barra de oscilação
 DELTA_P(k) = Pprog(k) - Pcalc_c(k);
 if abs(DELTA_P(k)) > deltaP %Estabelece a tolerância
 deltaP=DELTA_P(k);
 contP=cod_barra(k);
 end
 end
 if tipo(k)<=1 %considera somente as barras de carga
 DELTA_Q(k) = Qprog(k) - Qcalc_c(k);
 if abs(DELTA_Q(k)) > deltaQ %Estabelece a tolerância
 deltaQ=DELTA_Q(k);
 contQ=cod_barra(k);
end
 end
 end
 
 DELTA_PQ = [DELTA_P;DELTA_Q]; %Montagem da matriz de variações
 
 if abs(deltaP)>1e-2 || abs(deltaQ)>1e-2 %Condição para cálculo da Jacobiana
 
 %Cálculo da Jacobiana
 
 %Jacobiana=| H | N |
 % | M | L |
 
 H = zeros(num_barras,num_barras); %Cria as matrizes com suas dimensões
 M = zeros(num_barras,num_barras); %Cria as matrizes com suas dimensões
 N = zeros(num_barras,num_barras); %Cria as matrizes com suas dimensões
 L = zeros(num_barras,num_barras); %Cria as matrizes com suas dimensões
 
 %Cálculo das diagonais
 for k=1:num_barras
 H(k,k) = -Qcalc_c(k)-(V(k)^2)*B(k,k);
 if tipo(k)==3
 H(k,k) = 10^15; %Válido para barra de oscilação (fase constante)
 end
 N(k,k) = (Pcalc_c(k)+(V(k)^2)*G(k,k))/V(k) - (v2p(k)+2*v3p(k)*V(k)*Pesp(k));
 M(k,k) = Pcalc_c(k)-(V(k)^2)*G(k,k);
 L(k,k) = (Qcalc_c(k)-(V(k)^2)*B(k,k))/V(k) - (v2p(k)+2*v3p(k)*V(k)*Qesp(k));
 if tipo(k)>=2 %Válido para barra de tensão e oscilação (fase constante)
 L(k,k)=10^15;
 end
 end
 
 for a=1:num_trechos
 k = barra_origem(a);
 m = barra_destino(a);
 
 %Componentes de H (Derivada da potência ativa em relação aos ângulos)
 H(k,m) = H(k,m)+(V(k)^2)*(G(k,m)*sin(teta(k) - teta(m) + gama(a))-B(k,m)*cos(teta(k) - teta(m) + gama(a)));
 H(m,k) = H(m,k)+(V(k)*V(m)*(-G(k,m)*sin(teta(k) - teta(m) + gama(a))-B(k,m)*cos(teta(k) - teta(m) + gama(a))));
 
 %Componentes de N (Derivada da potência ativa em relação às tensões)
 N(k,m) = N(k,m)+(V(k)*(G(k,m)*cos(teta(k) - teta(m) + gama(a))+B(k,m)*sin(teta(k) - teta(m) + gama(a))));
 N(m,k) = N(m,k)+(V(m)*(G(k,m)*cos(teta(k) - teta(m) + gama(a))-B(k,m)*sin(teta(k) - teta(m) + gama(a))));
 
 %Componentes de M (Derivada da potência reativa em relação aos ângulos)
 M(k,m) = M(k,m)-V(k)*V(m)*(G(k,m)*cos(teta(k) - teta(m) + gama(a))+B(k,m)*sin(teta(k) - teta(m) + gama(a)));
 M(m,k) = M(m,k)-V(k)*V(m)*(G(k,m)*cos(teta(k) - teta(m) + gama(a))-B(k,m)*sin(teta(k) - teta(m) + gama(a)));
 
 %Componentes de L (Derivada da potência reativa em relação às tensões)
 L(k,m) = L(k,m)+V(k)*(G(k,m)*sin(teta(k) - teta(m) + gama(a))-B(k,m)*cos(teta(k) - teta(m) + gama(a)));
 L(m,k) = L(m,k)-V(m)*(G(k,m)*sin(teta(k) - teta(m) + gama(a))+B(k,m)*cos(teta(k) - teta(m) + gama(a)));
 
 end
 
 Jacobiana = [ H N ; M L ]; %Estruturação da matriz Jacobiana
 DELTA_V = inv(Jacobiana)*DELTA_PQ; %Cálculo da matriz de variação de V e teta
 
 V = V + DELTA_V(num_barras+1:2*num_barras)'; %atualização das tensões
 teta = teta + DELTA_V(1:num_barras)'; %atualização dos ângulos
 
 ite = ite+1;
 
 if ite>max_ite %Condição de limite de iterações (evitar loop infinito)
 deltaP=0;
 deltaQ=0;
 fprintf('Limite de iterações excedido');
 end
 end 
end %fim while 
%Resolução do segundo subsistema
%Cálculo das potências
for k=1:num_barras
 Pcalc_c(k) = G(k,k)*(V(k)^2);
 Qcalc_c(k) = -B(k,k)*(V(k)^2);
end
for a=1:num_trechos
 
 k = barra_origem(a);
 m = barra_destino(a);
 
 g=real(y(a));
 b=imag(y(a));
 
 Pcalc_c(k) = Pcalc_c(k) + V(k)*V(m)*(-g*cos(teta(k) - teta(m) + gama(a))-b*sin(teta(k) - teta(m) + gama(a)));
 Pcalc_c(m) = Pcalc_c(m) + V(k)*V(m)*(-g*cos(teta(k) - teta(m) + gama(a))+b*sin(teta(k) - teta(m) + gama(a)));
 Qcalc_c(k) = Qcalc_c(k) + V(k)*V(m)*(-g*sin(teta(k) - teta(m) + gama(a))+b*cos(teta(k) - teta(m) + gama(a)));
 Qcalc_c(m) = Qcalc_c(m) - V(k)*V(m)*(-g*sin(teta(k) - teta(m) + gama(a))-b*cos(teta(k) - teta(m) + gama(a)));
 
end
 
%Fluxo de potência
for a=1:num_trechos
 
 k = barra_origem(a);
 m = barra_destino(a);
 
 g=real(y(a));
 b=imag(y(a));
 
 %Potência entre as barras
 Pkm(a) = (V(k)^2)*g-V(k)*V(m)*(g*cos(teta(k) - teta(m) + gama(a))+b*sin(teta(k) - teta(m) + gama(a)));
 Pmk(a) = (V(m)^2)*g-V(k)*V(m)*(g*cos(teta(k) - teta(m) + gama(a))-b*sin(teta(k) - teta(m) + gama(a)));
 Qkm(a) = -(V(k)^2)*(b+bs(a))+V(k)*V(m)*(b*cos(teta(k) - teta(m) + gama(a))-g*sin(teta(k) - teta(m) + gama(a)));
 Qmk(a) = -(V(m)^2)*(b+bs(a))+V(k)*V(m)*(b*cos(teta(k) - teta(m) + gama(a))+g*sin(teta(k) - teta(m) + gama(a)));
 
end
%	Resultados Obtidos
%Tensões, ângulos e potências
fprintf('\nNúmero de iterações\n');
ite
fprintf('\nTensões, ângulos e potências\n');
fprintf('\nResultados Obtidos\n');
fprintf('\nTensão nas barras (pu)\n');
V
fprintf('\nÂngulo teta em graus\n');
teta=teta*180/pi;
teta
fprintf('\nPotência ativa calculada (pu)\n');
Pcalc_c
fprintf('\nPotência reativa calculada (pu)\n');
Qcalc_c
%Fluxo de potência
fprintf('\nFluxo de Potência\n');
fprintf('\n De Para Pkm Qkm\n');
for a = 1:num_trechos
	fprintf(' |%3d %6d %9.2f %8.2f |\n',barra_origem(a),barra_destino(a),Pkm(a)*100,Qkm(a)*100)
end

Teste o Premium para desbloquear

Aproveite todos os benefícios por 3 dias sem pagar! 😉
Já tem cadastro?