Post on 29-Jan-2016
1
Análisis de Flujo de Carga
Barra Tensión Angulo ------Carga------ ---Generación--- Shunt Mag. grados MW MVAr MW MVAr MVAr Carga_1 1.001 -2.938 200.0 30.0 0.0 0.0 0.0 Carga_2 1.029 -3.427 200.0 20.0 0.0 0.0 0.0 Carga_3 1.009 -13.732 100.0 30.0 0.0 0.0 0.0 Carga_4 0.893 -23.205 400.0 100.0 0.0 0.0 0.0 Gen_1 1.050 -0.709 0 0 500.0 161.3 0.0 Gen_2 1.050 -11.968 0 0 200.0 174.8 0.0 Slack 1.000 0.000 0.0 0.0 340.1 -22.6 0.0
Flujo en las líneas y pérdidas
--Línea-- -Flujo en la línea- --Pérdidas-- desde hasta MW Mvar MVA MW Mvar Carga_1 Carga_3 134.416 -28.964 137.501 4.205 -2.128 Carga_2 4.336 -41.077 41.305 0.156 -17.693 Carga_2 Carga_4 242.202 86.285 257.113 14.930 64.411 Carga_1 -4.180 23.384 23.754 0.156 -17.693 . . . . . . . . . . . . . .
SlackCarga_1 Carga_2
Carga_4Carga_3
Gen_2
Gen_1
PG
G
G
V0°|V|
|V|
Q
P
Q
P
QP
Q
Dada una red
Mediante resolución de las ecuaciones de flujo de carga determino las siguientes incognitas:
Una vez resueltas las barras, mediantes las ecuaciones fundamentales de circuitos, determino:
Presentación del problema
Luego aplicando en forma directa las ecuaciones de la red determino:
2
Expresiones fundamentales de la red
Vi
V1
V2
Vn
yi1
yi2
yin
yi0
Ii ...
ij
ij
n
i
nnninn
iniiii
ni
ni
n
i
niniiiiniiii
niiniiiiiii
y
y
V
V
V
V
YYYY
YYYY
YYYY
YYYY
I
I
I
I
VyVyVyVyyyyI
VVyVVyVVyVyI
ij
ii
2
1
21
21
222221
111211
2
1
2211210
22110
Y diagonal la de fuera elementos
Y diagonal la de elementos
:por dados estan elementos sus y
nodal admitancia matriz denomina le se arriba matriz La
.
.
.
.
....
....
....
....
....
....
....
....
.
.
.
.
.
:red una de barras n las para matricial forma en loexpresando
...)...(
:terminos los oreordenand
)(...)()(
n
jjiijijjii
n
jjiijijjii
n
jjijjijiiii
iiii
n
jjijjiji
n
jjiji
YVVQ
YVVP
VYVjQP
IVjQP
i
VYI
VYI
I
1
1
1
*
1
1
i
)(sin||||||
)(cos||||||
:imaginaria e real partes en Separando
)(||||)(||
:potencia la de expresión la en corriente la doSustituyen
:es barra la en compleja potencia La
)(||||
:tenemospolar forma en ecuación esta Expresando
:escribir Podemos
3
Clasificación de las barras de la red
Las barras son clasificadas generalmente en tres tipos:
• Barra Slack - Es tomada como referencia donde |V| y son especificados, no aporta ecuaciones al algoritmo, si no que una vez calculados los |V| y en el resto de las barras, se calcula Pslack y Qslack :
n
jjiijijjislack
n
jjiijijjislack
YVVQ
YVVP
1
1
)(sin||||||
)(cos||||||
• Barra de carga - o barra PQ, se especifica la potencia activa y reactiva, el módulo y la fase de las tensiones son desconocidas, y se calculan resolviendo el siguiente set de ecuaciones no lineares:
n
jjiijijjii
n
jjiijijjii
YVVQ
YVVP
1
1
)(sin||||||
)(cos||||||
• Barra de generación- o barra PV o barras de tensión controlada, se especifican el módulo de la tensión y la potencia activa, debiendose determinar la fase de la tensión y la potencia reactiva.Los límites de la potencia reactiva son también especificados. Se aplica entonces una única ecuación por barra para el cálculo de la fase de la tensión:
n
jjiijijjii YVVP
1
)(cos||||||
una vez calculadas todas los módulos y fases de las tensiones de todas las barras (o sea convergió algoritmo Newton-Raphson), se calcula Q en todas las barras PV:
n
jjiijijjii YVVQ
1
)(sin||||||
si se viola el límite inferior o superior en alguna/s barras se puede tomar alguna de las siguientes acciones correctivas: 1 - fijar Q=Qlim y liberar la tensión (transformar en una barra PQ) y vuelvo a entrar en el algoritmo N-R. 2 - Aumentar (o disminuir) un escalón porcentual el módulo de la tensión y vuelvo a entrar en el algoritmo N-R).
4
Datos de entrada para resolver el flujo de carga
% Datos de archivo de entrada tomados del Gross, pag. 244%% DATOS DE BARRA% CARGA GENERACION min max Shunt Shunt% BARRA TENSION MW MVAR MW MVAR MVAR MVAR MVAr SUCEPTANCIASL Slack 1 0 0 0 0 0 0 0 0PQ Carga_1 1 200 30 0 0 0 0 0 0PV Gen_1 1.05 60 8 500 0 0 0 0 0PQ Carga_2 1 200 20 0 0 0 0 0 0PV Gen_2 1.05 50 5 200 0 0 0 0 0PQ Carga_3 1 100 30 0 0 0 0 0 0PQ Carga_4 1 400 100 0 0 0 0 0 0%%% DATOS DE LINEAS% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIALinea Carga_1 Carga_3 0.023 0.138 0.271Linea Carga_2 Carga_4 0.023 0.138 0.271Linea Carga_1 Carga_2 0.015 0.092 0.181Linea Carga_3 Carga_4 0.015 0.092 0.181%%% DATOS DE TRANSFORMADORES% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA TAP Trafo Slack Carga_1 0.0012 0.015 1 Trafo Gen_1 Carga_2 0.001 0.012 1 Trafo Gen_2 Carga_3 0.002 0.024 1
SlackCarga_1 Carga_2
Carga_4
Carga_3Gen_2
Gen_1
PG
G
G
V0°
|V|
|V|
Q
P
Q
P
QP
Q
Dada una red
P
Q
P
Q
5
Solución de Ecuaciones Algebraicas No-Lineares - Método de Newton-Raphson
(k)
(k)(k)1)(k
)(
(k)(k)
(k)(k)
)0(
(0))0((1)
(0)
(0)(0)(0))0(
(0)
(0)
2(0)
)0(
2
2(0)
)0((0)
(0)(0)
(0)(0)
x Si
xxx
c
x
)x(c
: Definiendo
Raphson-Newton de algoritmo al identifica ntoprocedimie este de uso Sucesivo
c
x
:ónaproximaci segunda la en resultará inicial estimación la a x Sumando
residuo. el )x(c siendo ,xc
:resultando osdespreciadser pueden
alto orden de terminos los pequeño, muy es xerror el que Asumiendo
...)x(!2
1x)x(
:Taylor de serie en izquierda la de término el oExpandiend
)xx(
:tenemos correcta, solución la de desviación pequeña la es x y inicial, estimación la es x Si
)(
:por dada es ldiemnsiona-uni ecuación una de solución La
kJ
fc
dx
dfJ
dxdf
x
fcdx
df
cdx
fd
dx
dff
cf
cxf
Interpretación gráfica:
C=0
c(0)
x(0)
c(1)
x(1)
J(0)
J(1)
6
Ejemplo 6.1:
a) Búsqueda de la raíz de f(x)=x3-6 x2+9x-4.
cleardx=1; % Se inicializa el error con un valor elevadofun=input('Nombre de la función: '); % Nombre de la función.m donde están las expr. % de f y J. vx=input('Entre la estimación inicial y rango de ploteo [xe xi xf] -> '); x=vx(1);iter = 0; k=1;disp('iter Dc J dx x')% Encabezamiento de resultadoswhile abs(dx) >= 0.001 & iter < 100 % Test de convergencia iter = iter + 1; % No. de iteraciones [f,J]=feval(fun,x); % feval ejecuta la función especificada % en el string fun con el argumento x. yp(k)=f; % Puntos para graficar las xp(k)=x; % pendientes. Dc=0 - f; % Residuo dx= Dc/J; % Se actualiza el error x=x+dx; % Soluciones sucesivas yp(k+1)=0; % Puntos para graficar las xp(k+1)=x; % pendientes. k=k+2; fprintf('%g', iter) % Se muestra iter sin ceros % no significativos disp([Dc, J, dx, x]) % Se completa con el resto de las % variables.end
x=(vx(2):.1:vx(3)); % Rango de x para ploteo.f=feval(fun,x); % Se evalúa f en ese rangoplot(x,f,x,0*x,xp,yp) axis([vx(2) vx(3) min(f) max(f)]) % Se fijan los ejes para x y f.
function[f,J]=pol3(x)f=x.^3-6*x.^2+9*x-4;J=3*x.^2-12*x+9;
7
» te6ej1Nombre de la función: 'pol3'Entre la estimación inicial y rango de ploteo [xe xi xf] -> [6 0 6]iter Dc J dx x1 -50.0000 45.0000 -1.1111 4.8889
2 -13.4431 22.0370 -0.6100 4.2789
3 -2.9981 12.5797 -0.2383 4.0405
4 -0.3748 9.4914 -0.0395 4.0011
5 -0.0095 9.0126 -0.0011 4.0000
6 -0.0000 9.0000 -0.0000 4.0000
0 1 2 3 4 5 6
0
5
10
15
20
25
30
35
40
45
50
Búsqueda de la raíz de f(x)=x3-6 x2+9x-4.
8
Ejemplo 6.1:
b) Estudio de convergencia de f(x)=atg(x).
function[f,J]=atx(x)f=atan(x);J=1./(1+x.*x);
» te6ej1Nombre de la función: 'atx'Entre la estimación inicial y rango de ploteo [xe xi xf] -> [1.4 -20 20]
-20 -15 -10 -5 0 5 10 15 20-1.5
-1
-0.5
0
0.5
1
1.5
9
. : a denomina se , C
:compacta forma en
.
.
..
....
....
..
..
)(.
.
)(
)(
:matricial forma la en o
...)(
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
...)(
...)(
:ordenmayor de terminos los dodesprecian
Taylor, de series en izquierda la de termino el oExpandiend
),...,,(
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
),...,,(
),...,,(
:variables con ecuaciones ahora doConsideran
)0()0((0)
)0(
)0(2
)0(1
21
2
2
2
1
2
1
2
1
1
1
)0(
)0(22
)0(11
)0()0(2
2
)0(1
1
)0(
2)0(2)0(
22
2)0(1
1
2)0(2
1)0(1)0(
22
1)0(1
1
1)0(1
)0()0()0(2
)0(2
)0(1
)0(1
2)0()0()0(
2)0(
2)0(
1)0(
12
1)0()0()0(
2)0(
2)0(
1)0(
11
obianamatriz JacJXJ
x
x
x
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
x
f
fc
fc
fc
cxx
fx
x
fx
x
ff
cxx
fx
x
fx
x
ff
cxx
fx
x
fx
x
ff
cxxxxxxf
cxxxxxxf
cxxxxxxf
nn
n
n
nnn
n
n
nn
nnn
nnnn
nn
nn
nnnn
nn
nn
10
Quedando entonces el algoritmo de Newton-Raphson:
)0()0(2
)0(1 ,...,, iniciales valores los estiman Se nxxx
quemayor es X de elemento algún Si
X
X
Jacobiana la calcula Se
)(.
.
)(
)(
C
(k)
(k))()1(
)(1)((k)
)(
)(
)(22
)(11
(k)
kk
kk
k
knn
k
k
XX
CJ
J
fc
fc
fc
*
*El problema se reduce entonces a resolver sucesivos sistemas de ecuaciones lineares.
En Matlab, la solución del sistema de ecuaciones es obtenida
usando el operador de división de matrices \, o sea \ el cual es
basado en factorización triangular y eliminación Gaussiana, mucho más eficiente
que invertir
XJ C
CJ X
J .
11
Ejemplo 6.2:
Se usa el método de Newton-Raphson para encontrar la intersección de las curvas
1
422
ye
yxx
La siguiente rutina (te6ej2a) genera las gráficas
tita=0:.02:2*pi; % Rango del ángulo de la cfa.r = 2*ones(1, length(tita)); % Vector radio de la cfa.x=-3:.02:1.5; % Rango de x para la segunda ec.y=1- exp(x); % Segunda ec.plot(x,y),gridaxis([-3 3 -3 3]);axis('square'); % Relación de ejes tal que no deformen la cfa.xlabel('x')text(1.1,1.8,' x^2+y^2=4')text(1.2,-2.3,' e^x+y=1')hold on; % Se "fija" la gráfica tal que las sucesivas % se hagan en la misma figura con los mismos ejes.polar(tita, r) % Ploteo polar en un sistema cartesiano.hold off; % Se "libera" la figura
-3 -2 -1 0 1 2 3-3
-2
-1
0
1
2
3
x
y
x2+y2=4
ex+y=1
12
Tomando las derivadas parciales, la matriz Jacobiana resulta:
1
22xe
yxJ
La siguiente rutina (te6ej2b) aplica Newton-Raphson para el sistema arriba
iter = 0; x=input('Entre el vector estimación inicial [x1; x2] -> ');Dx=[1; 1];C=[4; 1];disp('Iter DC Matriz Jacobiana Dx x');while max(abs(Dx)) >= .0001 & iter < 100 iter=iter+1; f = [x(1)^2+x(2)^2; exp(x(1))+x(2)]; DC = C - f; J = [2*x(1) 2*x(2) exp(x(1)) 1]; Dx=J\DC; % Resolución del sistema de ecuaciones x=x+Dx; fprintf('%g', iter) disp([DC, J, Dx, x])end
» te6ej2bEntre el vector estimación inicial [x1; x2] -> [0.5 -1]'Iter DC Matriz Jacobiana Dx x1 2.7500 1.0000 -2.0000 0.8034 1.3034 0.3513 1.6487 1.0000 -0.9733 -1.9733
2 -1.5928 2.6068 -3.9466 -0.2561 1.0473 -0.7085 3.6818 1.0000 0.2344 -1.7389
3 -0.1205 2.0946 -3.4778 -0.0422 1.0051 -0.1111 2.8499 1.0000 0.0092 -1.7296
4 -0.0019 2.0102 -3.4593 -0.0009 1.0042 -0.0025 2.7321 1.0000 0.0000 -1.7296
5 -0.0000 2.0083 -3.4593 -0.0000 1.0042 -0.0000 2.7296 1.0000 -0.0000 -1.7296
13
Tenemos entonces dos ecuaciones por cada barra PQ y una por cada barra PV, suponiendo que:Barra 1 - barra SlackBarra 2 a m - barras PQBarras m+1 a n - barras PVExpandiendo en series de Taylor haciendo estimaciones iniciales para |V| y y despreciando los términos de orden elevado, se llega al siguiente sistema de ecuaciones lineares:
||.
.
||
.
.
||..
||..
......
......||
..||
..
||..
||..
......
......||
..||
..
.
.
.
.
)(
)(2
)(
)(2
)(
2
)()(
2
)(
)(2
2
)(2
)(2
2
)(2
)(
2
)()(
2
)(
)(2
2
)(2
)(2
2
)(2
)(
)(2
)(
)(2
km
k
kn
k
m
km
km
n
km
km
m
kk
n
kkm
kn
kn
n
kn
kn
m
kk
n
kk
km
k
kn
k
V
V
V
Q
V
QQQ
V
Q
V
QQQV
P
V
PPP
V
P
V
PPP
Q
Q
P
P
En forma abreviada:
||43
21
VJJ
JJ
Q
P
El procedimiento para solucionar un flujo de carga con el método de Newton-Raphson es el que sigue:
Para las barras PQEspecifica Pi y Qi
Estima |Vi(0)| y (0) (igual a la slack)
Para las barras PVEspecifica Pi , |Vi| y los limites max y min de Qi
Estima (0) (igual a la slack)
Usando los valoresespecificados y estimados Calculo el vector:
)(
)(2
)(
)(2
.
.
.
.
km
k
kn
k
Q
Q
P
P
n
j
jk
ik
ijijk
jk
iiespk
i
n
j
jk
ik
ijijk
jk
iiespk
i
n
j
jk
ik
ijijk
jk
iiespk
i
YVVPP
YVVQQ
YVVPP
1
)()()()(.
)(
1
)()()()(.
)(
1
)()()()(.
)(
)(cos||||||
:PV barras las para y
)(sin||||||
)(cos||||||
:PQ barras las Para
Para la barra Slack
Se especifica V y
14
n
j
jk
ik
ijijk
jk
iik
i
n
j
jk
ik
ijijk
jk
iik
i
n
j
jk
ik
ijijk
jk
iik
i
YVVPP
YVVQQ
YVVPP
QP
1
)()()()()(
1
)()()()()(
1
)()()()()(
)(cos||||||
:PV barras las para y
)(sin||||||
)(cos||||||
:PQ barras las Para
: y los devetor el Actualizo
Se calculan los elementos de la matriz jacobiana J1, J2, J3 y J4.
Q
P
JJ
JJ
V\
|| 43
21Se resuelve:
Se actualizan los |Vi| y i :|||||| )()()1(
)()()1(
ki
ki
ki
ki
ki
ki
VVV
Mientras halla algún:|Pi
(k)|> o algún|Qi
(k)|>
convergió
n
jjslackijijjslackslack
n
jjslackijijjslackslack
YVVQ
YVVP
1
1
)(sin||||||
)(cos||||||
PV) (barras a 1 de
)(sin||||||1
nmipara
YVVQn
jjiijijjii
Calculo la potencia reactiva en todas las barras PV: Si se violó al menos un límite tomo accióncorrectiva y vuelvo al algoritmo
15
Solución Flujo de Carga Desacoplado Rápido
Para relación X/R alta
P está fuertemente acoplado a y debilmente acoplado a |V|
||||
||
o
||0
0
4
1
4
1
VV
QVJQ
PJP
VJ
J
Q
P
Además considerables simplificaciones a J1 y J4 pueden ser hechas:
)(sin||||)(sin||||||)(sin||||||
:diagonal la de elementos los Para
)(cos||||||
Siendo
2
1,1
1
1
iiiii
n
jjiijijji
n
ijjjiijijji
i
i
n
jjiijijjii
YVYVVYVVP
PJ
YVVP
-QiBii
Siendo Bii la parte imaginaria de los elementos de la diagonal de Y, o sea, la suma de todas lassuceptancias incidentes a la barra i.
iii
i
BV
Q
||P
entonces
|V||V| además ,B típicos potencia de sistemas En
i
i
i2
iii
Q está fuertemente acoplado a |V| y debilmente acoplado a
16
ijij
jijjij
jiji
jiijji
BV
VYVV
YVV
||P
1|| asumiendo obtenida es ciónsimplifica otra ),sin(||||||P
)( entonces ,0 normal operación de scondicione en
ij )sin(||||||P
:diagonal la de fuera elementos los Para
i
iji
ijij
ijj
i
Bii
ijij
i
iiii
i
BVV
Q
BVV
Q
J
||||
:diagonal fuera
||||
:diagonal
: parasimilar forma En 4
Llegamos entonces a que los sistemas de ecuaciones||
||||4
1
VV
QVJQ
PJP
Se pueden plantear como:
PQ. barras las a solo ientescorrespond los y
PQ yPV barras las a ientescorrespond Y de elementos los de asuceptanci la Siendo
||||
||
''
'
''
'
B
B
VBV
Q
BV
P
Siendo B’ y B’’ constantes, estas pueden ser invertidas una única vez antes de iniciar las iteracionesy luego durante el proceso de cálculo los cambios de |V| y son dados en forma directa por:
||
||
||
1''
1'
V
QBV
V
PB
17
El procedimiento para solucionar un flujo de carga con el método de Newton-Raphson desacoplado rápido es el que sigue:
Para las barras PQEspecifica Pi y Qi
Estima |Vi(0)| y (0) (1.00)
Para las barras PVEspecifica Pi , |Vi| y los limites max y min de Qi
Estima (0) (1.00)
||
.
.||
)(
2
)(2
n
kn
k
V
P
V
P
Determinar B’ y B’’ y en consecuencia [B’]-1 y [B’’]-1
Para la barra Slack
Se especifica V y
Usando los valoresespecificados y estimados Calculo los vectores:
n
j
jk
ik
ijijk
jk
iiespk
i
n
j
jk
ik
ijijk
jk
iiespk
i
n
j
jk
ik
ijijk
jk
iiespk
i
YVVPP
YVVQQ
YVVPP
1
)()()()(.
)(
1
)()()()(.
)(
1
)()()()(.
)(
)(cos||||||
:PV barras las para y
)(sin||||||
)(cos||||||
:PQ barras las Para
||
.
.||
)(
2
)(2
m
km
k
V
Q
V
Q
18
||]''[||
||]'[
1
1
V
QBV
V
PB
n
j
jk
ik
ijijk
jk
iik
i
n
j
jk
ik
ijijk
jk
iik
i
n
j
jk
ik
ijijk
jk
iik
i
YVVPP
YVVQQ
YVVPP
QP
1
)()()()()(
1
)()()()()(
1
)()()()()(
)(cos||||||
:PV barras las para y
)(sin||||||
)(cos||||||
:PQ barras las Para
: y los Actualizo
Se actualizan los |Vi| y i :|||||| )()()1(
)()()1(
ki
ki
ki
ki
ki
ki
VVV
Mientras halla algún:|Pi
(k)|> o algún|Qi
(k)|>
convergió
n
jjslackijijjslackslack
n
jjslackijijjslackslack
YVVQ
YVVP
1
1
)(sin||||||
)(cos||||||
PV) (barras a 1 de
)(sin||||||1
nmipara
YVVQn
jjiijijjii
Calculo la potencia reactiva en todas las barras PV: Si se violó al menos un límite tomo accióncorrectiva y vuelvo al algoritmo
|| ||
los ActualizoV
Q
V
P :productos los mediante
|| y calculan Se V
19
Implementación Flujo de Carga
Archivo texto conla configuración de la red
Matrices “usables”por el Matlab
Flujo de Carga
Matriz Admitancia
Resultados
Proceso
Entrada
Salida
[Y][V,P,Q]
20
SlackCarga_1 Carga_2
Carga_4Carga_3Gen_1
Gen_1
PG
G
G
V0°|V|
|V|
% Datos de archivo de entrada tomados del Gross, pag. 244%% DATOS DE BARRA% CARGA GENERACION min max Shunt Shunt% BARRA TENSION MW MVAR MW MVAR MVAR MVAR MVAr SUCEPTANCIASL Slack 1 0 0 0 0 0 0 0 0PQ Carga_1 1 200 30 0 0 0 0 0 0PV Gen_1 1.05 0 0 500 0 0 0 0 0PQ Carga_2 1 200 20 0 0 0 0 0 0PV Gen_2 1.05 0 0 200 0 0 0 0 0PQ Carga_3 1 100 30 0 0 0 0 0 0PQ Carga_4 1 400 100 0 0 0 0 0 0%%% DATOS DE LINEAS% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIALinea Carga_1 Carga_3 0.023 0.138 0.271Linea Carga_2 Carga_4 0.023 0.138 0.271Linea Carga_1 Carga_2 0.015 0.092 0.181Linea Carga_3 Carga_4 0.015 0.092 0.181%%% DATOS DE TRANSFORMADORES% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA TAP Trafo Slack Carga_1 0.0012 0.015 1 Trafo Gen_1 Carga_2 0.001 0.012 1 Trafo Gen_2 Carga_3 0.002 0.024 1
Q
P
Q
P
QP
Q
21
1.0000 1.0000 1.0000 200.0000 30.0000 0 0 2.0000 2.0000 1.0000 200.0000 20.0000 0 0 3.0000 3.0000 1.0000 100.0000 30.0000 0 0 4.0000 4.0000 1.0000 400.0000 100.0000 0 0
5.0000 5.0000 1.0500 0 0 500.0000 0 6.0000 6.0000 1.0500 0 0 200.0000 0 7.0000 7.0000 1.0000 0 0 0 0 1.0000 3.0000 0.0230 0.1380 0.2710 0 0 2.0000 4.0000 0.0230 0.1380 0.2710 0 0 1.0000 2.0000 0.0150 0.0920 0.1810 0 0 3.0000 4.0000 0.0150 0.0920 0.1810 0 0 7.0000 1.0000 0.0012 0.0150 1.0000 0 0 5.0000 2.0000 0.0010 0.0120 1.0000 0 0 6.0000 3.0000 0.0020 0.0240 1.0000 0 0
0 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0
N=
1 4 5 6 7 7 8 11 12 14 0 0
pN=
'Carga_1' 'Carga_2' 'Carga_3' 'Carga_4' 'Gen_1' 'Gen_2' 'Slack'
Barras=
% Datos de archivo de entrada tomados del Gross, pag. 244%% DATOS DE BARRA% CARGA GENERACION min max Shunt Shunt% BARRA TENSION MW MVAR MW MVAR MVAR MVAR MVAr SUCEPTANCIASL Slack 1 0 0 0 0 0 0 0 0PQ Carga_1 1 200 30 0 0 0 0 0 0PV Gen_1 1.05 0 0 500 0 0 0 0 0PQ Carga_2 1 200 20 0 0 0 0 0 0PV Gen_2 1.05 0 0 200 0 0 0 0 0PQ Carga_3 1 100 30 0 0 0 0 0 0PQ Carga_4 1 400 100 0 0 0 0 0 0%%% DATOS DE LINEAS% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIALinea Carga_1 Carga_3 0.023 0.138 0.271Linea Carga_2 Carga_4 0.023 0.138 0.271Linea Carga_1 Carga_2 0.015 0.092 0.181Linea Carga_3 Carga_4 0.015 0.092 0.181%%% DATOS DE TRANSFORMADORES% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA TAP Trafo Slack Carga_1 0.0012 0.015 1 Trafo Gen_1 Carga_2 0.001 0.012 1 Trafo Gen_2 Carga_3 0.002 0.024 1
Se usa cuando se especifica datos generadores (estabilidad)
22
Función red2mat :
Esta función, convierte un archivo ASCII con los datos de la red, en matrices usablespor el Matlab
Argumentos de entrada:• Archivo ASCII con los datos de la red.• Nombre del archivo para almacenar los datos de salida (opcional).
Argumentos de salida:• Matrices N, pN, Barras.
23
El procedimiento para solucionar un flujo de carga con el método de Newton-Raphson desacoplado rápido es el que sigue:
Para las barras PQEspecifica Pi y Qi
Estima |Vi(0)| y (0) (1.00)
Para las barras PVEspecifica Pi , |Vi| y los limites max y min de Qi
Estima (0) (1.00)
||
.
.||
)(
2
)(2
n
kn
k
V
P
V
P
Determinar B’ y B’’ y en consecuencia [B’]-1 y [B’’]-1
Para la barra Slack
Se rspecifica V y
Usando los valoresespecificados y estimados Calculo los vectores:
n
j
jk
ik
ijijk
jk
iiespk
i
n
j
jk
ik
ijijk
jk
iiespk
i
n
j
jk
ik
ijijk
jk
iiespk
i
YVVPP
YVVQQ
YVVPP
1
)()()()(.
)(
1
)()()()(.
)(
1
)()()()(.
)(
)(cos||||||
:PV barras las para y
)(sin||||||
)(cos||||||
:PQ barras las Para
||
.
.||
)(
2
)(2
m
km
k
V
Q
V
Q
24
||]''[||
||]'[
1
1
V
QBV
V
PB
n
j
jk
ik
ijijk
jk
iik
i
n
j
jk
ik
ijijk
jk
iik
i
n
j
jk
ik
ijijk
jk
iik
i
YVVPP
YVVQQ
YVVPP
QP
1
)()()()()(
1
)()()()()(
1
)()()()()(
)(cos||||||
:PV barras las para y
)(sin||||||
)(cos||||||
:PQ barras las Para
: y los Actualizo
Se actualizan los |Vi| y i :|||||| )()()1(
)()()1(
ki
ki
ki
ki
ki
ki
VVV
Mientras halla algún:|Pi
(k)|> o algún|Qi
(k)|>
convergió
n
jjslackijijjslackslack
n
jjslackijijjslackslack
YVVQ
YVVP
1
1
)(sin||||||
)(cos||||||
PV) (barras a 1 de
)(sin||||||1
nmipara
YVVQn
jjiijijjii
Calculo la potencia reactiva en todas las barras PV: Si se violó al menos un límite tomo accióncorrectiva y vuelvo al algortimo
|| ||
los ActualizoV
Q
V
P :productos los mediante
|| y calculan Se V
25
Función ybb :
Esta función, dados los datos de la red, obtiene las matrices Y, B’ y B’’ .
Argumentos de entrada:• Matrices N y pN con los datos de la red.
Argumentos de salida:• Matriz admitancia Y y las matrices suceptancias B’ y B’’.
Red ejemplo:
0.02+j0.04
0.01+j0.03 0.0125+j0.025
400 MW
250 MVar
slack
barra_PV
barra_PQ1
200 MW
1.050° Vpu
|1.04| Vpu
barra_PQ2
% DATOS DE BARRA% CARGA GENERACION min max Shunt% BARRA TENSION MW MVAR MW MVAR MVAR MVAR MVArSL slack 1.05 0 0 0 0 0 0 0PQ barra_PQ1 1 0 0 0 0 0 0 0PQ barra_PQ2 1 400 250 0 0 0 0 0PV barra_PV 1.04 0 0 200 0 0 0 0%% DATOS DE LINEAS% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIALinea slack barra_PQ1 0.02 0.04 0.001Linea slack barra_PV 0.01 0.03 0Linea barra_PQ1 barra_PV 0.0125 0.025 0%% DATOS DE TRANSFORMADORES% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA TAP Trafo barra_PQ2 barra_PQ1 0.012 0.015 0.98
26
Ejecutando [N,pN]=fcm2dat(‘ejemplo4b.m’)
[Y,Bp,Bs]=ybb(N,pN)
N =
1.0000 1.0000 1.0000 0 0 0 0 0 0 0 0 2.0000 2.0000 1.0000 400.0000 250.0000 0 0 0 0 0 0 3.0000 3.0000 1.0400 0 0 200.0000 0 0 0 0 0 4.0000 4.0000 1.0500 0 0 0 0 0 0 0 0 4.0000 1.0000 0.0200 0.0400 0.0010 0 0 0 0 0 0 4.0000 3.0000 0.0100 0.0300 0 0 0 0 0 0 0 1.0000 3.0000 0.0125 0.0250 0 0 0 0 0 0 0 2.0000 1.0000 0.012 0.0150 0.9800 0 0 0 0 0 0
pN =
1 2 3 3 4 4 5 7 8 8
Y =
59.9617 -94.3265i -33.1840 +41.4800i -16.0000 +32.0000i -10.0000 +20.0000i -33.1840 +41.4800i 32.5203 -40.6504i 0 0 -16.0000 +32.0000i 0 26.0000 -62.0000i -10.0000 +30.0000i -10.0000 +20.0000i 0 -10.0000 +30.0000i 20.0005 -50.0000i
» Bp
Bp =
-94.3265 41.4800 32.0000 41.4800 -40.6504 0 32.0000 0 -62.0000
» Bs
Bs =
-94.3265 41.4800 41.4800 -40.6504
27
yt
Vi VjVx
1:aIi Ij
j
i
tt
tt
j
i
ji
jt
it
j
iij
jt
iti
x
xiti
i
ji
jx
V
V
a
y
a
ya
yy
I
I
II
Va
yV
a
yI
IIa
I
Va
yVyI
V
VVyI
I
aII
Va
V
2
2
:matricial forma en y de sexpresione las oescribiendRe
:a llegamos ecuación, esta en ndosubstituye ,1
:que tenemos También
:tenemos , ndoSubstituye
)(
:por dada es corriente La
1
:tenemos asumidas, corriente de sdireccione las Para
Se considera entonces un transformador con admitancia yt en serie con un transformador ideal de relación 1:a siendo a el valor por unidad de la posición del tap, el que permitirá pequeños ajustes de tensión en escalones del orden del 1 al 1.5% de la tensión nominal, el modelo del transformador ahora es :
ayt
tya
a
1
tya
a
2
1
28
Función flunrdr:
Esta función, dados los datos de la red, ejecuta flujo de carga por el métodoNewton-Raphson desacoplado rápido.
Argumentos de entrada:• Matrices N y pN con los datos de la red o archivo.dat generado por fcm2dat.• Nombre del archivo de salida para almacenar los resultados (opcional).
Argumentos de salida:• Para cada una de las barras: Módulo y ángulo de la tensión.
Potencia activa y reactiva de carga. Potencia activa y reactiva de generación.
Compensación.• General: Máximo error de potencia una vez que convergió. Número de iteraciones.
29
Mientras iter<iter máx:
Desde Barra 1 hasta última Barra PV
p f(Y, |V|, )
También convergió
Q?
Todos los Pson < que ?
si
Hay BarrasPV?
Busco en cual de las Barras PV hay violación de los límites
si
si
Hubo violaciónen algún límite?
Aumento la tensión en la barra donde se violó el límite inferior y bajo la tensión donde se violó el límite superior.Aviso que es como si no hubiera convergido
Fuerzo una salida del while.Avisando que es “normal”
si
no
no
Calculo resolviendo el sistema.Actualizo
P Pesp-p calcPV P/|V|
A
B
no
Ano
30
Desde Barra 1 hasta última Barra PV
q f(Y, |V|, )
También convergió
P?
Todos los Qson < que ?
si
Hay BarrasPV?
Busco en cual de las Barras PV hay violación de los límites
si
si
Hubo violaciónen algún límite?
Aumento la tensión en la barra donde se violó el límite inferior y bajo la tensión donde se violó el límite superior.Aviso que es como si no hubiera convergido
Fuerzo una salida del while.Avisando que es “normal”
si
si
no
no
Calculo |V| resolviendo el sistemaActualizo |V|.
A
Q Qesp-q calcQV Q/|V| (solo barras PQ)
Iter=iter+1B
Warning: No convergióen iter max iteraciones
FIN
no
31
» [N,pN,Barras]=red2mat('ejemplo4b.m');» [mv,an,Pd,Qd,Pg,Qg,Qsh,maxerror,iter]=flunrdr(N,pN);» diary ejemplo.tab» tabflu
Flujo en las líneas y pédidas
--Línea-- -Flujo en la línea- --Pérdidas-- desde hasta MW Mvar MVA MW Mvar barra_PQ1 slack -190.357 -114.738 222.262 10.691 21.282 barra_PV -243.851 -178.033 301.926 12.333 24.667 barra_PV slack -56.184 -15.582 58.304 0.314 0.943 barra_PQ1 256.184 202.700 326.677 12.333 24.667 slack barra_PQ1 201.048 136.020 242.738 10.691 21.282 barra_PV 56.498 16.524 58.865 0.314 0.943 Total de pérdidas 23.339 46.891
Flujo en las transformadores
-Transformador- -Flujo transformador- --Tap-- desde hasta MW Mvar MVA barra_PQ1 barra_PQ2 434.196 292.673 471.653 0.980 barra_PQ2 barra_PQ1 -399.994 -249.922 523.626 0.980
Utilizando el archivo del ejemplo anterior:
32
Función flujo:
Ejecuta todos los comandos para realizar un flujo de carga y desplegar los resultados.Argumentos de entrada:• Nombre del archivo ASCII con los datos de la red.
function[]=flujo(archivo)global SbSb=100;[N,pN,Barras]=red2mat(archivo);[mv,an,Pd,Qd,Pg,Qg,Qsh,maxerror,iter,Y]=flunrdr(N,pN);
archivo=[strtok(archivo,'.') '.flu']; % crea el nombre del archivo de salída % igual al de entrada pero con extensión .tabdiary(archivo)tabbartabfludiary
» flujo('ejemplo4b.m') Máximo error en la potencia = 0.0981324 No. de Iteraciones = 11
Barra Tensión Angulo ------Carga------ ---Generación--- Shunt Mag. grados MW MVAr MW MVAr MVAr barra_PQ1 0.961 -3.022 0.0 0.0 0.0 0.0 0.0 barra_PQ2 0.883 -5.006 400.0 250.0 0.0 0.0 0.0 barra_PV 1.040 -0.803 0.0 0.0 200.0 187.1 0.0 slack 1.050 0.000 0.0 0.0 257.5 152.5 0.0 Total 400.0 250.0 457.5 339.7 0.0
Flujo en las líneas y pédidas
--Línea-- -Flujo en la línea- --Pérdidas-- desde hasta MW Mvar MVA MW Mvar barra_PQ1 slack -190.357 -114.738 222.262 10.691 21.282 barra_PV -243.851 -178.033 301.926 12.333 24.667 barra_PV slack -56.184 -15.582 58.304 0.314 0.943 barra_PQ1 256.184 202.700 326.677 12.333 24.667 slack barra_PQ1 201.048 136.020 242.738 10.691 21.282 barra_PV 56.498 16.524 58.865 0.314 0.943 Total de pérdidas 23.339 46.891
Flujo en las transformadores
-Transformador- -Flujo transformador- --Tap-- desde hasta MW Mvar MVA barra_PQ1 barra_PQ2 434.196 292.673 471.653 0.980 barra_PQ2 barra_PQ1 -399.994 -249.922 523.626 0.980
33
Es interesante aprovechar el potencial del Matlab para estudios paramétricos, por ejemplo seconsiderando el siguiente archivo (pag84.m) que corresponde a la red simplificada de Nueva Zelandia:
% Datos de archivo de entrada tomados del Arrillaga, cap. IV%% DATOS DE BARRA% CARGA GENERACION min max Shunt% BARRA TENSION MW MVAr MW MVAR MVAr MVAr MVArPQ INV220 1 200 51 0 0 0 0 0PQ ROX220 1 150 60 0 0 0 0 0 PQ MAN220 1 0 0 0 0 0 0 0PV MAN014 1.045 0 0 690 0 0 0 0PQ TIW220 1 420 185 0 0 0 0 0SL ROX011 1.05 0 0 0 0 0 0 0PQ BEN220 1 500 200 0 0 0 0 0 PV BEN016 1.06 0 0 0 0 0 0 0 PQ AVI220 1 0 0 0 0 0 0 0PV AVI011 1.045 0 0 200 0 0 0 0 PV OHAU 1.05 0 0 350 0 0 0 0 PQ LIV220 1 150 60 0 0 0 0 0PV ISL220 1.0 500 300 0 0 0 0 0PQ BRM220 1 100 60 0 0 0 0 0PQ TEK011 1.05 0 0 150 0 0 0 0PQ TEK220 1 0 0 0 0 0 0 0 PQ TWI220 1 0 0 0 0 0 0 0%%% DATOS DE LINEAS% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA SUCEPTANCIALinea INV220 MAN220 0.013 0.09 0.25Linea INV220 MAN220 0.013 0.09 0.25Linea MAN220 TIW220 0.01 0.1 0.29Linea MAN220 TIW220 0.01 0.1 0.29Linea INV220 TIW220 0.002 0.01 0.04Linea INV220 TIW220 0.002 0.01 0.04Linea INV220 ROX220 0.01 0.11 0.17 Linea ROX220 TWI220 0.016 0.14 0.24Linea ROX220 TWI220 0.016 0.14 0.24Linea ROX220 LIV220 0.03 0.12 0.18Linea BEN220 TWI220 0.004 0.03 0.07Linea LIV220 AVI220 0.007 0.03 0.05Linea AVI220 BEN220 0.004 0.05 0.02Linea AVI220 BEN220 0.004 0.05 0.02Linea LIV220 ISL220 0.03 0.18 0.35Linea TWI220 TEK220 0.002 0.01 0.02Linea TEK220 ISL220 0.02 0.13 0.35Linea TWI220 BRM220 0.02 0.14 0.45Linea BRM220 ISL220 0.002 0.01 0.05Linea TWI220 ISL220 0.02 0.14 0.45%%% DATOS DE TRANSFORMADORES% BARRA_1 BARRA_2 RESISTENCIA REACTANCIA TAP Trafo MAN220 MAN014 0.0006 0.016 1 Trafo ROX220 ROX011 0.002 0.04 1 Trafo TWI220 OHAU 0.004 0.032 1 Trafo AVI220 AVI011 0.0015 0.045 1 Trafo BEN220 BEN016 0.0012 0.032 1 Trafo TEK220 TEK011 0.003 0.056 1
34
Para variar un parámetro de una barra o línea o transformador, necesitamos saber dicho elementoen que fila está de la matriz N, una vez ubicada la fila podemos hacer el estudio paramétrico variandoel o los valores de las columnas que nos interesan.
Función filaN:
Esta función encuentra la fila de N donde está la barra o línea o transformador requerido.
Argumentos de entrada:• Nombre de una o dos barras.• Vector celda con los nombres de todas las barras.• Matriz N (No es necesaria si se ingresó solo una barra)
Argumentos de salida:• Fila o filas de la matriz N donde se encuentra la barra o línea o transformador de interés.
function[fN]=filaN(b1,b2,Barras,N)
fl=0;if iscell(b2), % Si el segundo argumento de entrada es una celda Barras=b2; % es porque se ingreso una única barra, por lo tanto fl=1; % b2 es en realidad el vector celda con los nombres deend % todas las barras.
nb1=find(strcmpi(Barras,b1)); % Se busca el número de barra que corresponde b1if fl==1, % Si se ingresó una única barras, el número de esta se corresponde fN=nb1; % con la fila, pronto, me voy de la función. returnendnb2=find(strcmpi(Barras,b2)); % Se busca el número que le corresponde a la barra b2
fb=(N(:,1:2)==nb1); % fb es un matriz de dos columnas x filas de N, donde hay unos % cuando coincide el nombre de la barra 1
fb=fb+2*(N(:,1:2)==nb2); % A la matriz anterior le sumo dos cuando hay coincidencia % con el nombre de la barra 2.
fN=find(sum(fb')==3); % Las filas de fN que sumen tres son las buscadas.
2
1
3
4
5
35
Por ejemplo: cual es la ubicación dentro de la matriz N de la barra AVI220:
» [N,pN,Barras]=fcm2dat('pag84.m');» f=filan('AVI220',Barras)
f =
6
Para el caso de líneas: cual es la ubicación dentro de la matriz N de la línea TEK220-ISL220:
» f=filan('TEK220','ISL220',Barras,N)
f =
34
Siguiendo por pasos la ejecución de filaN para este último ejemplo:
nb1 = 10
1 -
2 - nb2 = 16
3 -
fb = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
fb = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 1 2 0 0 0 2 0 2 0 0 0 0 0 0 0 0 0 0 1 0
fN = 34
4 - 5 -
36
Ejemplo 6.3a:
Dada la red definida en ‘pag84.m’, el siguiente ejemplo estudia el módulo de la tensióny el ángulo en la barra ‘INV220’ al variar la carga en esta de 0 a 500 MW y de 0 a 250 MVAr
pinv=linspace(0,500,5); % Rango potencia activa (5 puntos)qinv=linspace(0,250,5); % Rango potencia reactiva (5 puntos)
clear v af=0;c=0;
global Sb;Sb=100;
[N,pN,Barras]=red2mat('pag84.m');INV220=filaN('INV220',Barras); % Se 'ubica' a INV220 en N
for var1=pinv, % Variación de la potencia activa c=c+1; N(INV220,4)=var1; % Se modifica el elemento de N correspondiente % a la potencia activa de INV220 for var2=qinv, N(INV220,5)=var2; % Se modifica el elemento de N correspondiente % a la potencia reactiva de INV220 f=f+1; [mv an]=flunrdr(N,pN); % Flujo de cargas v(f,c)=mv(INV220); % Se alacena módulo de la tensión a(f,c)=an(INV220); % y el ángulo end f=0;end
[x,y]=meshgrid(pinv,qinv); % Se crea el grid para el plote 3Dsurf(x,y,v,a) % Ploteo de la superficie: v contra x,yshading interp % Interpolación de coloresview([40 35]) % Variación del ángulo de la visióncolorbar('v') % Escala de colores de a (ángulo) en pos. verticalxlabel('Potencia Activa')ylabel('Potencia Reactiva')zlabel(Tensión)title('Estudio de Flujo de Carga sobre la Barra INV220')
37
38
Ejemplo 6.3b:
Idéntico al anterior pero usando flunrdrI.
pinv=linspace(0,500,5); % Rango potencia activa (5 puntos)qinv=linspace(0,250,5); % Rango potencia reactiva (5 puntos)
clear v af=0;c=0;
load pag84.dat -mat % Se cargan los datos previamente calculados con fcm2datINV220=filaN('INV220',Barras); % Se 'ubica' a INV220 en N[Y,Bp,Bs]=ybb(N,pN);
for var1=pinv, % Variación de la potencia activa c=c+1; N(INV220,4)=var1; % Se modifica el elemento de N correspondiente % a la potencia activa de INV220 for var2=qinv, N(INV220,5)=var2; % Se modifica el elemento de N correspondiente % a la potencia reactiva de INV220 f=f+1; [mv an]=flunrdrI(N,pN,Y,Bp,Bs); % Flujo de cargas v(f,c)=mv(INV220); % Se alacena módulo de la tensión a(f,c)=an(INV220); % y el ángulo end f=0;end
[x,y]=meshgrid(pinv,qinv); % Se crea el grid para el plote 3Dsurf(x,y,v,a) % Ploteo de la superficie: v contra x,yshading interp % Iterolación de coloresview([40 35]) % Variación del ángulo de la visióncolorbar('v') % Escala de colores de a (angulo) en pos. verticalxlabel('Potencia Activa')ylabel('Potencia Reactiva')zlabel('Tension')title('Estudio de Flujo de Carga sobre la Barra INV220')