Baixe o app para aproveitar ainda mais
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
Compartilhar