Código Flujo de carga Newton Raphson

UNIVERSIDAD INDUSTRIAL DE SANTANDER
Escuela de Ingenierías Eléctrica, Electrónica y de Telecomunicaciones
Perfecta Combinación entre Energía e Intelecto
ESTUDIO FLUJO DE CARGAS

Cristian Leonardo Manrique Pérez
PROFESOR: GERARDO LATORRE BAYONA
SISTEMAS DE POTENCIA

Universidad Industrial de Santander
Facultad de Ingeniería Físico Mecánicas
Escuela De Ingeniería Eléctrica, Electrónica Y Telecomunicaciones
Bucaramanga
2014

PLANTEAMIENTO

Flujo de cargas y cortocircuito para el sistema de la figura 1.

Figura 1. Sistema de potencia

Procedimiento paso a paso para resolver el problema por el método de Newton Raphson
El primer paso es calcular la Ybarra del sistema sin incluir cargas ni generadores, el segundo paso es determinar los tipos de barra del sistema cuáles son las variables de entrada (grados de libertad) y cuáles son las incógnitas (resultados) ver tabla 1, el tercer paso determinar los datos de potencias especificadas en cada barra y su respectivo calculo en por unidad, el cuarto paso es hacer una descripción del sistema ver tabla 2, el quinto paso por medio de proceso iterativo determinar el Jacobiano, el vector de correcciones, hacer las correcciones en los ángulos y en las tensione, para un factor de tolerancia de 1,0*10^-3, el sexto paso con los resultados de tensiones y ángulos determinar corrientes inyectadas en las barras, flujo de corrientes por las líneas, flujo de potencias en las líneas, perdidas en las líneas, potencias inyectadas en las barras y potencias generadas en las barras.

Sistema flujo de cargas
TIPO DE BARRA
BARRA
GENERACIÓN
CARGA
TENSIÓN DE BARRA
Tap t
INCÓGNITAS

Pg [MW]
Qg [MVAr]
Pc [MW]
Qc [MVAr]
V pu
Delta

slack
1
_
_
0
0
1
0
1
Pg1; Qg1
PQ
2
0
0
0
0
1
0
1
V2; ɗ2
PQ
3
0
0
40,00
19,37
1
0
1
V3; ɗ3
PQ
4
0
0
0
0
1
0
1
V4; ɗ4
PQ
5
0
0
0
0
1
0
1
V5; ɗ5
PQ
6
18,02
_
0
0
1
0
1
V6; ɗ6
PV
7
18,02
_
0
0
1
0
1
Qg7; ɗ7
PV
8
0
0
7,00
3,92
1
0
1
Qg8; ɗ8

Tabla 1. Datos sistema flujo de cargas

DESCRIPCIÓN DEL SISTEMA
Número de Barras
8
Número de Transformadores
2
Número Líneas
2
Número de cargas
2
Niveles de tensión [KV]
Vb1 [KV]
220

Vb2 [KV]
110

Vb3 [KV]
44
SB [MVA]
100
Zb1 [Ὼ]
121
Zb2 [Ὼ]
19,36
IDENTIFICACIÓN BARRAS
CANTIDAD
Slack
1
1
PV
7-8
2
PQ
2-3-4-5-6
5
PQV
0
-
Ecuaciones de P
7

Ecuaciones de Q
5

Sistema de ecuaciones
12
12X12
Número de incógnitas
16

Tabla 2. Descripción del sistema

Cálculos que se deben realizar antes de iniciar el proceso iterativo
Para el cálculo de la Ybarra:
Para la línea de 110 KV, encontrar la capacitancia que se presenta en la línea, tomando como referencia el tipo de conductor y la disposición de las torres que de la línea de transmisión. Se supone que no existen perdidas por efecto corona por ello se asume la perditancia igual a cero.

DISTANCIA MEDIA GEOMÉTRICA MUTUA
Entre fases
Dm=Deq=3D12*D23*D31 [m]
D12,D23,D31=distancia entre fases
DISTANCIA MEDIA GEOMÉTRICA INDUCTIVA PROPIA
Para un conductor por fase.
Dsh=RMG [cm]
RMG: Radio Medio Geométrico
DISTANCIA MEDIA GEOMÉTRICA CAPACITIVA PROPIA
Para un conductor por fase.
DsC=r [cm]
r: Radio del conductor

INDUCTANCIA POR FASE
L=2*10-4LnDm Dsh [H/km/fase]
REACTANCIA INDUCTIVA
XL=2*π*f*L [ /km/fase]
f: Frecuencia
CAPACITANCIA POR FASE
C=18*109*LnDm DsC-1*1000 [F/km/fase]
REACTANCIA CAPACITIVA
Xc=12*π*f*C [ /km/fase]
SUCEPTANCIA CAPACITIVA
B=1Xc [S/km]
CONDUCTANCIA O PERDITANCIA
Gk=PPV2*10-3 [S/km]
Pp: Pérdidas de potencia en kilo-Vatios por kilometro
V: Tensión por fase en kilo-voltios

RESISTENCIA EN [Ω/KM] DEL CONDUCTOR
Rc =RL=10ρS [ /Km]
Este valor es dado en tablas del conductor penguin
IMPEDANCIA SERIE POR UNIDAD DE LONGITUD POR FASE

Z=Rc+jXL [

Rc: Resistencia en [Ω/Km] del Conductor

ADMITANCIA EN PARALELO POR UNIDAD DE LONGITUD POR FASE A NEUTRO
Y=G+jB [S/km/fase]
Cable Penguin
AWG
Diametro [mm]
RMG [mm]
Peso [kg/Km]
Resistencia 75º C [Ohm/km]
Capacidad [A]
Capacidad CC [KA]
4/0
14,31
4,61
433
0,396
355
16,1
R1(+) [Ohm/km]
0,396
Dm [m]
7,3296

X1(+) [Ohm/km]
0,5558
L [H/Km]
0,00147429

B1(+) [Ohm/km]
3,0214E-06
XL [Ω/km]
0,555794183

R(o) [Ohm/km]
1,188
C [F/Km]
8,01452E-09

X(o) [Ohm/km]
1,66738255
B [S/Km]
3,0214E-06

XBo) [Ohm/km]
9,06421E-06

Tabla 3. Cálculos línea 110 KV

Se supuso para la línea de 44 KV sin efecto capacitivo
Las impedancias en por unidad de los transformadores y de las líneas en bases del sistema

De la Barra
A la Barra
Tensión [KV]
Longitud [Km]
R1
X1
B1
CAP Normal [A]
Capacidad de la líneas Limite Térmico [MVA]
1
2
_
_
0
0,0597
0
472,377
180
2
3
_
_
0
-0,0086
0
944,755
180
2
4
_
_
0
0,1536
0
787,296
60
3
5
110
85
0,2782
0,3904
0,0155
355
39,05
4
8
44
3
0,0236
0,0465
0,0000
300
13,2
5
6
_
_
0
0,302
0
131,216
25
5
7
_
_
0
0,302
0
131,216
25
Tabla 4. Impedancia en las bases del sistema
Para el flujo de cargas:
Determinar la potencia en las barras en por unidad
Determinar la dimensión de las sub-matrices del Jacobiano

DIMENSIÓN DE SUB-MATRICES DEL JACOBIANO (FILAS X COLUMNA)
H
7
7
N
7
5
J
5
7
L
5
5
Jacobiano

H
N

J
L

Tabla 5. Jacobiano
Se inició con un flat clásico la iteración (tensiones y ángulos de las barras en 1 pu y ángulo en 0 grados respectivamente)
Para cada Iteraciones determinar:
Valores de potencia calculada
Los errores de los elementos del Jacobiano
Las correcciones y los valores corregidos

Datos arrojados en las iteraciones por matlab
____________________________________
Numero de iteraciones
1

Vector de Errores en pu
-0.0000
-0.4000
-0.0000
-0.0000
0.1802
0.1802
-0.0700
-0.0000
-0.1315
0
0.0622
-0.0392

Tolerancia
0.4000

Vector de Correccion en angulos en grados
-0.3749
-0.3554
-0.9909
9.8888
13.0069
13.0069
-1.1244

Valores Corregidos de los angulos en grados
0
-0.3749
-0.3554
-0.9909
9.8888
13.0069
13.0069
-1.1244

Vector de Correccion en tension
-0.0185
-0.0162
-0.0245
0.0306
-0.0280

Valores Corregidos de tensiones en pu
1.0000
0.9815
0.9838
0.9755
1.0306
1.0000
1.0000
0.9720

Potencias Activa calculada en las barras
1.0e-11 *

0.0205
0.1424
0.1488
0.0533
0.0063
0.0041
0.0041
0.0178

Potencias Reactiva calculada en las barras
0
0.0000
-6.2200
0
-6.2200
0
0
0

2

0.0023
-0.0192
-0.0002
-0.0131
-0.0054
-0.0054
-0.0029
-0.0007
-0.0350
-0.0010
-0.0443
-0.0007

Tolerancia
0.0443

-0.1390
-0.1191
-0.1705
-0.5740
-0.6460
-0.6460
-0.1785

0
-0.5139
-0.4745
-1.1615
9.3149
12.3608
12.3608
-1.3029

-0.0020
-0.0017
-0.0023
-0.0062
-0.0024

1.0000
0.9795
0.9821
0.9732
1.0245
1.0000
1.0000
0.9696

10.7573
-0.2312
-38.0767
0.0198
1.3054
18.5631
18.5631
-6.7070

30.9995
0.0706
-15.8747
0.0956
4.4302
-9.6420
-9.6420
-3.8541

3

0.0001
0.0002
-0.0001
-0.0004
-0.0001
-0.0001
0.0001
-0.0006
0.0012
-0.0001
-0.0020
0.0002

Tolerancia
0.0020

-0.0006
-0.0005
-0.0008
-0.0042
-0.0045
-0.0045
-0.0008

0
-0.5145
-0.4750
-1.1622
9.3107
12.3564
12.3564
-1.3037

1.0e-03 *

0.0263
0.0183
0.0453
-0.2303
0.0569

1.0000
0.9795
0.9821
0.9732
1.0243
1.0000
1.0000
0.9697

14.7163
-0.0086
-40.0206
0.0098
0.0395
18.0259
18.0259
-7.0076

34.3767
0.0582
-19.4861
0.0075
0.1975
-7.6287
-7.6287
-3.9388

4

1.0e-04 *

-0.0012
0.0801
0.0433
-0.0678
-0.0100
-0.0100
-0.0417
0.1959
-0.1653
0.0434
-0.4753
-0.0852

Tolerancia
4.7531e-05

1.0e-06 *

-0.2631
-0.2283
-0.2709
-0.6497
-0.5518
-0.5518
-0.2687

0
-0.5145
-0.4750
-1.1622
9.3107
12.3564
12.3564
-1.3037

1.0e-05 *

-0.0556
-0.0338
-0.1255
-0.5560
-0.1792

1.0000
0.9795
0.9821
0.9732
1.0243
1.0000
1.0000
0.9697

14.7346
0.0000
-40.0008
-0.0004
0.0007
18.0201
18.0201
-6.9996

34.3329
-0.0020
-19.3683
-0.0004
0.0048
-7.5527
-7.5527
-3.9191
___________________________________

Cálculos de las potencias generadas en las barras Slack, PV y PQ.

POTENCIA INYECTADAS EN LAS BARRAS Y GENERADAS EN LAS BARRAS
" # Barra " Pinyectada " Qinyectada " Pgenerada " Qgenerada "
___________________________________________________________________
1 14.7346 34.3338 14.7346 34.3338
___________________________________________________________________
2 -0.0000 0.0001 -0.0000 0.0001
___________________________________________________________________
3 -40.0000 -19.3701 -0.0000 -0.0001
___________________________________________________________________
4 0.0000 0.0000 0.0000 0.0000
___________________________________________________________________
5 0.0000 0.0001 0.0000 0.0001
___________________________________________________________________
6 18.0200 -7.5508 18.0200 -7.5508
___________________________________________________________________
7 18.0200 -7.5508 18.0200 -7.5508
___________________________________________________________________
8 -7.0000 -3.9200 -0.0000 -0.0000
______________________________________________________________________

Cálculo de flujos de potencia en las líneas, flujo de corriente por las líneas y por los transformadores
-----Flujo de potencias por las Líneas--------

Flujo de Potencia Activa por las líneas MW
" # Barras " Potencia Activa "
(2,1) -14.7346
(1,2) 14.7346
(3,2) -7.7184
(4,2) -7.0162
(2,3) 7.7184
(5,3) 36.0400
(2,4) 7.0162
(8,4) -7.0000
(3,5) -32.2816
(6,5) 18.0200
(7,5) 18.0200
(5,6) -18.0200
(5,7) -18.0200
(4,8) 7.0162

Flujo de Potencia Reactiva líneas MVAr
" # Barras " Potencia Reactiva "
(2,1) -33.5005
(1,2) 34.3338
(3,2) -29.5266
(4,2) -3.9518
(2,3) 29.4435
(5,3) -10.8819
(2,4) 4.0570
(8,4) -3.9200
(3,5) 16.1561
(6,5) -7.5508
(7,5) -7.5508
(5,6) 8.7037
(5,7) 8.7037
(4,8) 3.9519

FLUJO DE CORRIENTES POR LAS LINEAS
# Barras " IL pu " Angulo en Grados
_______________________________
1 2 0.3736 -66.7731
_______________________________
2 3 0.3107 -75.8253
_______________________________
2 4 0.0827 -30.5526
_______________________________
3 5 0.3676 26.1118
_______________________________
4 8 0.0827 -30.5527
_______________________________
5 6 0.1954 35.0913
_______________________________
5 7 0.1954 35.0913
___________________________

Cálculo de las pérdidas en el sistema
Perdidas de Potencias en las lineas Avtiva – Reactiva
" # Barras " PpL MW " PqL MVAr "
_______________________________
1 2 0.0000 0.8334
_______________________________
2 3 0.0000 -0.0830
_______________________________
2 4 -0.0000 0.1052
_______________________________
3 5 3.7584 5.2742
_______________________________
4 8 0.0162 0.0318
_______________________________
5 6 0.0000 1.1528
_______________________________
5 7 0.0000 1.1528
_______________________________

Cálculo de tensiones en las barras y tensiones inducidas en los generadores
Newton Raphson Analisis de flujos en sistemas de potenccia
# Barra " V pu " Angulo en Grados
_________________________
1 1.0000 0.0000
_________________________
2 0.9795 -0.5145
_________________________
3 0.9821 -0.4750
_________________________
4 0.9732 -1.1622
_________________________
5 1.0243 9.3107
_________________________
6 1.0000 12.3564
_________________________
7 1.0000 12.3564
_________________________
8 0.9697 -1.3037

Tensiones y angulo de tensiones inducida en los generadores en las barras
# Barra " Eg pu " Angulo en Grados
_________________________
2 0.0000 0.0000
_________________________
3 0.0000 0.0000
_________________________
4 0.0000 0.0000
_________________________
5 0.0000 0.0000
_________________________
6 1.0670 21.4866
_________________________
7 1.0670 21.4866
_________________________
8 0.0000 0.0000

Código en Matlab:
clear all
clc
%% Newton Raphson Analisis de flujo de cargas sistema Diseñado de 20 barras
% Cristian Leonardo Manrique Perez
% Sistemas de Potencia
% PROFESOR: GERARDO LATORRE BAYONA
% UIS 2014
%% Y Barra
% "De la Barra "A la Barra" R pu " X pu " B/2 pu "X' tap pu "
Lineas =[1 2 0 0.0597 0 1;
2 3 0 -0.0086 0 1;
2 4 0 0.1536 0 1;
3 5 0.2782 0.3904 0.0622 1;
4 8 0.0236 0.0465 0 1;
5 6 0 0.302 0 1;
5 7 0 0.302 0 1];
% Barra de inicio
Ib= Lineas(:,1);
% Barra Fnal
Fb= Lineas(:,2);
% Resistencia
r = Lineas(:,3);
% Reactancia
x = Lineas(:,4);
% Suceptancia
b = Lineas(:,5);
% Tap valor
a = Lineas(:,6);
% Z impedancias
z = r + i*x;
% Inverso de impedancia ADMITANCIA
y = 1./z;
% Suceptancia
b = i*b;
% numero de barras
nb = max(max(Ib),max(Fb));
% Numero de Lineas
nl = length(Ib);
% Matriz Basia Y barra, Valor inicial
Y = zeros(nb,nb);
% MATRIZ YBARRA
for k = 1:nl
Y(Ib(k),Fb(k)) = Y(Ib(k),Fb(k)) - y(k)/a(k);
Y(Fb(k),Ib(k)) = Y(Ib(k),Fb(k));
end
for m = 1:nb
for n = 1:nl
if Ib(n) == m
Y(m,m) = Y(m,m) + y(n)/(a(n)^2) + b(n);
elseif Fb(n) == m
Y(m,m) = Y(m,m) + y(n) + b(n);
end
end
end
% Matrix Ybarra
Y;
% Matrix Zbarra
Z=inv(Y);
%% Datos de las barras
% Tipo de Barra
% 1 - Slack
% 2 - PV
% 3 - PQ
% Tipo " Barra " Pg[MW] " Qg[MVAr] " Pc[MW] " Qc[MVAr] " V pu " Delta " Qmin " Qmax "
Datos=[1 1 0 0 0 0 1.0 0 0 0;
3 2 0 0 0 0 1.0 0 0 0;
3 3 0 0 40 19.37 1.0 0 0 0;
3 4 0 0 0 0 1.0 0 0 0;
3 5 0 0 0 0 1.0 0 0 0;
2 6 18.02 0 0 0 1.0 0 -11.2 13.5;
2 7 18.02 0 0 0 1.0 0 -11.2 13.5;
3 8 0 0 7.00 3.92 1.0 0 0 0];
%% Datos del sistema en pu Sistema en por uniddad
% Base sistema [MVA]
Ba = 100;
% Tipo de Barra
type = Datos(:,1);
% Numero de Barra
bus = Datos(:,2);
% Pg pu
Pg = Datos(:,3)/Ba;
% Qg pu
Qg = Datos(:,4)/Ba;
% Pc pu
Pc = Datos(:,5)/Ba;
% Qc pu
Qc = Datos(:,6)/Ba;
% Voltaje especificado en las barras
V = Datos(:,7);
% Angulo delta de las barras
del = Datos(:,8);
% Limite minimo de reactiva
Qmin = Datos(:,9)/Ba;
% Limite maximo de reactiva
Qmax = Datos(:,10)/Ba;
%% Sistemas de Ecuaciones
% Pi = Pgi - Pci
P = Pg - Pc;
% Qi = QGi - Qci
Q = Qg - Qc;
% P especificada
Psp = P;
% Q especificada
Qsp = Q;
% Matrx de conductacias
G = real(Y);
% Matrix de Suceptancias
B = imag(Y);
%% Identificacion de tipo de barras
% PV ----------------------(Tomando como pv la barra Slack)
pv = find(type == 2 " type == 1);
% PQ ---(Tomando como PQ la barra PQV, para un valor determinado del tap t)
pq = find(type == 3);
% Numero de barras PV
npv = length(pv);
% Numeros de barras PQ
npq = length(pq);
%% Jacobiano Metodo Newton Raphson
nbus =nb;
Tol = 1;
Iter = 1;
while (Tol > 0.1*10e-3) % Valor de toleracia
% disp('Numero de iteraciones')
% disp(Iter)
P = zeros(nbus,1);
Q = zeros(nbus,1);
% calculo de P y Q
for i = 1:nbus
for k = 1:nbus
P(i) = P(i) + V(i)*abs(Y(i,k))*V(k)*cos(del(i)-del(k)-angle(Y(i,k)));
Q(i) = Q(i) + V(i)*abs(Y(i,k))*V(k)*sin(del(i)-del(k)-angle(Y(i,k)));
end
end
% Verificndo los limites de Reactiva
if Iter <= 7 && Iter > 2 % Mirando unicamente siete iteraciones
for n = 2:nbus
if type(n) == 2
QG = Q(n)+Qc(n);
if QG < Qmin(n)
V(n) = V(n) + 0.01;
elseif QG > Qmax(n)
V(n) = V(n) - 0.01;
end
end
end
end
% Calculo de la correcion de las variables
dPa = Psp-P;
dQa = Qsp-Q;
k = 1;
dQ = zeros(npq,1);
for i = 1:nbus
if type(i) == 3
dQ(k,1) = dQa(i);
k = k+1;
end
end
dP = dPa(2:nbus);
M = [dP; dQ]; % Vector de Errores
% Jacobiano
% H - Derivadas de las potecias inyectadas al sistema variables - angulos
H = zeros(nbus-1,nbus-1);
for i = 1:(nbus-1)
m = i+1;
for k = 1:(nbus-1)
n = k+1;
if n == m
for n = 1:nbus
H(i,k) = H(i,k) - (V(m)*abs(Y(m,n))*V(n)*sin(del(m)-del(n)-angle(Y(m,n))));
end
H(i,k) = H(i,k) - V(m)^2*B(m,m);
else
H(i,k) = V(m)*abs(Y(m,n))*V(n)*sin(del(m)-del(n)-angle(Y(m,n)));
end
end
end
% disp('H')
% disp(H)
% N - Derivada de las potecias inyectadas al sistema con tensiones
N = zeros(nbus-1,npq);
for i = 1:(nbus-1)
m = i+1;
for k = 1:npq
n = pq(k);
if n == m
for n = 1:nbus
N(i,k) = N(i,k) + V(m)*abs(Y(m,n))*V(n)*cos(del(m)-del(n)-angle(Y(m,n)));
end
N(i,k) = N(i,k) + V(m)^2*G(m,m);
else
N(i,k) = V(m)*abs(Y(m,n))*V(n)*cos(del(m)-del(n)-angle(Y(m,n)));
end
end
end
% disp('N')
% disp(N)
% J - Derivadas de las potecias reactivas inyectadas al sistema variables - angulos
J = zeros(npq,nbus-1);
for i = 1:npq
m = pq(i);
for k = 1:(nbus-1)
n = k+1;
if n == m
for n = 1:nbus
J(i,k) = J(i,k) + V(m)*abs(Y(m,n))*V(n)*cos(del(m)-del(n)-angle(Y(m,n)));
end
J(i,k) = J(i,k) - V(m)^2*G(m,m);
else
J(i,k) = -(V(m)*abs(Y(m,n))*V(n)*cos(del(m)-del(n)-angle(Y(m,n))));
end
end
end
% disp('J')
% disp(J)
% L - Derivadas de las potecias reactivas inyectadas al sistema variables - tensiones
L = zeros(npq,npq);
for i = 1:npq
m = pq(i);
for k = 1:npq
n = pq(k);
if n == m
for n = 1:nbus
L(i,k) = L(i,k) + V(m)*abs(Y(m,n))*V(n)*sin(del(m)-del(n)-angle(Y(m,n)));
end
L(i,k) = L(i,k) - V(m)^2*B(m,m);
else
L(i,k) = V(m)*abs(Y(m,n))*V(n)*sin(del(m)-del(n)-angle(Y(m,n)));
end
end
end
% disp('L')
% disp(L)
% Matrix JACOBIANO
JAC = [H N; J L];
% Vector de Correccion
X = inv(JAC)*M;
% Correcion de las variables de estado
dTh = X(1:nbus-1);
dV = X(nbus:end);

del(2:nbus) = dTh + del(2:nbus); % Angulo en las tensiones de las barras
k = 1;
for i = 2:nbus
if type(i) == 3
V(i) = dV(k) + V(i); % magnitud de las tensiones de barra
k = k+1;
end
end

% disp('Vector de Errores en pu')
% disp(M)
Iter = Iter + 1;
Tol = max(abs(M)); % evaluacion de la toleracia
% disp('Tolerancia ')
% disp(Tol)
% disp('Vector de Correccion en angulos en grados')
% disp(dTh*180/pi)
% disp('Valores Corregidos de los angulos en grados')
% disp( 180/pi*del)
% disp('Vector de Correccion en tension')
% disp(dV)
% disp('Valores Corregidos de tensiones en pu')
% disp(V)
% disp('Potencias Activa calculada en las barras')
% disp(P*Ba)
% disp('Potencias Reactiva calculada en las barras')
% disp(Q*Ba)
end
disp('Numero de iteraciones')
disp(Iter)
Del = 180/pi*del;
%% Tensiones y angulo de tensiones en las barras
disp(' Newton Raphson Analisis de flujos en sistemas de potenccia ');
disp(' # Barra " V pu " Angulo en Grados ');
for m = 1:nb
disp('_________________________');
fprintf('%3g', m); fprintf(' %8.4f', V(m)); fprintf(' %8.4f', Del(m));
fprintf('\n');
end
disp('___________________________');
%%
% Funcion de conversion de polar a rectangular
Vm=V.*cos(del) + j*V.*sin(del);
% Corrientes inyectadas a las barras
I1 = Y*Vm;
Im = abs(I1);
Ia = angle(I1)*180/pi;
%% Corrientes Inyectadas a las BARRAS
disp(' Corrientes Inyectadas a las BARRAS');
disp(' # Barra " I pu " Angulo en Grados ');
for m = 1:nb
if Ia(m) >= 90
Ia(m) = Ia(m)-180;
elseif Ia(m) <= -90
Ia(m) = Ia(m)+180;
else
Ia(m)=Ia(m);
end
disp('_________________________');
fprintf('%3g', m); fprintf(' %8.4f', Im(m)); fprintf(' %8.4f', Ia(m));
fprintf('\n');
end
disp('___________________________');

%% Flujos en las barras y en las lineas
% Matrices iniciales
Iij = zeros(nb,nb);
Sij = zeros(nb,nb);
Si = zeros(nb,1);
% Flujos de corrientes por las lineas
disp('FLUJO DE CORRIENTES POR LAS LINEAS')
disp(' # Barras " IL pu " Angulo en Grados ');
for m = 1:nl
p = Ib(m); q = Fb(m);
Iij(p,q) =-(Vm(p) - Vm(q))*Y(p,q); % Y(m,n) = -y(m,n)=1/Z %
Iij(q,p) = -Iij(p,q);
Iijm = abs(-(Vm(p) - Vm(q))*Y(p,q));
Iija = angle(-(Vm(p) - Vm(q))*Y(p,q))*180/pi;

if Iija >= 90
Iija = Iija-180;
elseif Iija <= -90
Iija = Iija+180;
else
Iija=Iija;
end
IijL= [Iijm,Iija];
disp('_______________________________');
fprintf('%3g', p,q); fprintf(' %8.4f', IijL);
fprintf('\n');

end
disp('___________________________');
%% Flujo de potencias en las lineas
for m = 1:nb
for n = 1:nb
if m ~= n
Sij(m,n) = Vm(m)*conj(Iij(m,n))*Ba;
end
end
end
Sij = sparse(Sij);
disp('-----Flujo de potencias por las Lineas--------')
Pij = real(Sij);
disp('Flujo de Potencia Activa por las lineas')
disp('" # Barras " Potencia Activa "');
disp(Pij);
disp('Flujo de Potencia Reactiva por las lineas')
Qij = imag(Sij);
disp('" # Barras " Potencia Reactiva "');
disp(Qij)
%Perdidas en las lineas
Lij = zeros(nl,1);
disp('Perdidas de Potencias en las lineas Avtiva - Reactiva')
disp('" # Barras " PpL " PqL "');
for m = 1:nl
p = Ib(m); q = Fb(m);
Lij(m) = Sij(p,q) + Sij(q,p);
Lpij =real(Lij); Lqij =imag(Lij);
disp('_______________________________');
fprintf(' %3g', p , q); fprintf(' %8.4f', Lpij(m)); fprintf(' %8.4f', Lqij(m));
fprintf('\n');
end
disp('_______________________________');
%potencia inyectada en las barras
disp('POTENCIA INYECTADAS EN LAS BARRAS Y GENERADAS EN LAS BARRAS')
disp(' " # Barra " Pinyectada " Qinyectada " Pgenerada " Qgenerada "');
for i = 1:nb
for k = 1:nb
Si(i) = Si(i) + conj(Vm(i))* Vm(k)*Y(i,k);
end
Pi = real(Si)*Ba;
Qi = -imag(Si)*Ba;
Pg = Pi+Pc*Ba;
Qg = Qi+Qc*Ba;

disp('___________________________________________________________________');
fprintf(' %3g', i); fprintf(' %8.4f', Pi(i)); fprintf(' %8.4f', Qi(i));fprintf(' %8.4f', Pg(i)) ;fprintf(' %8.4f', Qg(i));
fprintf('\n');

end
disp('______________________________________________________________________');
%% TENSIONES INDUCIDAS EN LOS GENERADORES
Eg = zeros(nb,1);
Ig = zeros(nb,1);
Id = zeros(nb,1);
xd=j*0.941;
xq=j*0.7;
rs=0.00343;
for k = 1:npv
m = pv(k);
Ig(m)=Ig(m)+(conj(Si(m))/(V(m).*cos(del(m)) - j*V(m).*sin(del(m))));
Id(m)=Id(m)+abs(Ig(m))*cos(pi+angle(Ig(m))-del(m))*(cos(-pi+del(m)) + j*sin(-pi+del(m)));
Eg(m)=Eg(m)+Id(m)*(xd-xq)+Ig(m)*(rs+xq)+(V(m).*cos(del(m)) + j*V(m).*sin(del(m)));

end
%% Tensiones y angulo de tensiones inducida en los generadores en las barras
disp('Tensiones y angulo de tensiones inducida en los generadores en las barras')
disp(' # Barra " Eg pu " Angulo en Grados ');
for m = 2:nb
disp('_________________________');
fprintf('%3g', m); fprintf(' %8.4f', abs(Eg(m))); fprintf(' %8.4f', angle(Eg(m))*180/pi);
fprintf('\n');
end
disp('___________________________');

Código Flujo de carga Newton Raphson

Recommend Documents