Matrices dispersas y Ybus con transformaciones lineales

En esta ocasión traigo un breve tutorial realizado en iPython de como calcular la matriz Ybus sin la necesidad de ciclos que sólo entorpecen el programa. Se utilizan matrices de conectividad y transformaciones lineales para llegar al objetivo. La teoría la he documentado en un pdf.

Formación de la matriz Ybus

In [1]:

# Importamos las librerías necesarias
from cvxopt import matrix, spmatrix, sparse, mul, div, exp
import numpy as np
from cmath import pi

Datos del sistema. En este caso se utilizarán los datos del sistema de 5 nodos documentado en el libro de Stagg. Todos los valores que se utilizan están en por unidad (pu)

In [2]:

# Número de nodos del sistema
m = 5
# Número de líneas
nl = 7
# Resistencia de las líneas
R = matrix([0.02, 0.08, 0.06, 0.06, 0.04, 0.01, 0.08])
# Reactancia inductiva de las líneas
X = matrix([0.06, 0.24, 0.18, 0.18, 0.12, 0.93, 0.24])
# Un medio de la suceptancia de carga en las líneas
B = matrix([0.03, 0.025, 0.02, 0.02, 0.015, 0.01, 0.025])
# En este caso no se tienen inyecciones de potencia activa ni reactiva. Las matrices serán de ceros
Bs = matrix(0,(m,1),'d')
Gs = matrix(0,(m,1),'d')
# Vector complejo (se requiere que sea complejo más adelante) 
# con el tap en el que se encuentra el transformador. Para el caso de líneas es 1
Tap = matrix(1,(nl,1),'z')
# Angulo de defasamiento para los transformadores defasadores. Para el caso de líneas es 0
Angulo = matrix(0,(nl,1),'d')
# Lista que indica si está activo el elemento. 1=Activo, 0= Inactivo o abierto
Activo = matrix(1,(nl,1),'d')
# Índices de los nodos de envío. Se utiliza n-1, siendo n el nodo real, debido a la indexación 
# de Pyton que comienza desde el 0 y no desde el 1
NP = [0, 0, 1, 1, 1, 2, 3]
# Índices de los nodos de llegada
NQ = [1, 2, 2, 3, 4, 3, 4]

Se calcula la admitancia serie de cada elemento:

In [3]:

Ys = div(Activo, (R + 1j * X)) 
print Ys

[ 5.00-j15.00]
[  1.25-j3.75]
[  1.67-j5.00]
[  1.67-j5.00]
[  2.50-j7.50]
[  0.01-j1.08]
[  1.25-j3.75]

La susceptancia de carga:

In [4]:

Bc = mul(Activo, B)
print Bc

[ 0.03]
[ 0.03]
[ 0.02]
[ 0.02]
[ 0.01]
[ 0.01]
[ 0.03]

Se aplica el defase debido a los transformadores defasadores

In [5]:

Tap = mul(Tap, exp(Angulo * pi * 1j / 180))
print Tap

[ 1.00-j0.00]
[ 1.00-j0.00]
[ 1.00-j0.00]
[ 1.00-j0.00]
[ 1.00-j0.00]
[ 1.00-j0.00]
[ 1.00-j0.00]

Calulando la admitancia propia de los nodos de envío

In [6]:

Ypp = Ys + 1j * Bc
print Ypp

[ 5.00-j14.97]
[  1.25-j3.73]
[  1.67-j4.98]
[  1.67-j4.98]
[  2.50-j7.49]
[  0.01-j1.07]
[  1.25-j3.73]

Admitancia propia de los nodos de llegada. La función A.H.T realiza una transpuesta conjugada y después vuelve a transponer resultando el conjugado de la matriz A.

In [7]:

Yqq = div(Ypp, mul(Tap, Tap.H.T))
print Yqq

[ 5.00-j14.97]
[  1.25-j3.73]
[  1.67-j4.98]
[  1.67-j4.98]
[  2.50-j7.49]
[  0.01-j1.07]
[  1.25-j3.73]

Admitancia mutua del nodo de envio al nodo de llegada

In [8]:

Ypq = -div(Ys, Tap.H.T)
print Ypq

[-5.00+j15.00]
[ -1.25+j3.75]
[ -1.67+j5.00]
[ -1.67+j5.00]
[ -2.50+j7.50]
[ -0.01+j1.08]
[ -1.25+j3.75]

Admitancia mutua del nodo de llegada al nodo de envio:

In [9]:

Yqp = -div(Ys, Tap)
print Yqp

[-5.00+j15.00]
[ -1.25+j3.75]
[ -1.67+j5.00]
[ -1.67+j5.00]
[ -2.50+j7.50]
[ -0.01+j1.08]
[ -1.25+j3.75]

Considerando que las inyecciones de potencia nodales están sobre la base de 1Vpu:

In [10]:

Ysh = (Gs + 1j * Bs)
print Ysh

[ 0.00-j0.00]
[ 0.00-j0.00]
[ 0.00-j0.00]
[ 0.00-j0.00]
[ 0.00-j0.00]

Aquí empezamos a manejar matrices dispersas. Las cuales se crean de la forma spmatrix(valores, ind_renglones, ind_columnas, (n,m)). Donde el elemento i-ésimo de ind_columnas es la columna donde estará el valor del elemento i-ésimo. De igual forma para los renglones. "n" es el número de renglones de la matriz y "m" el número de columnas. Creamos una matriz de conectividad para los nodos de envio:

In [11]:

Cp = spmatrix(1, range(nl), NP, (nl, m))
print Cp

[ 1.00   0     0     0     0  ]
[ 1.00   0     0     0     0  ]
[  0    1.00   0     0     0  ]
[  0    1.00   0     0     0  ]
[  0    1.00   0     0     0  ]
[  0     0    1.00   0     0  ]
[  0     0     0    1.00   0  ]

Ahora la matriz de conectividad para los nodos de llegada:

In [12]:

Cq = spmatrix(1, range(nl), NQ, (nl, m))
print Cq

[  0    1.00   0     0     0  ]
[  0     0    1.00   0     0  ]
[  0     0    1.00   0     0  ]
[  0     0     0    1.00   0  ]
[  0     0     0     0    1.00]
[  0     0     0    1.00   0  ]
[  0     0     0     0    1.00]

Se crea un doble arreglo con los índices de las lineas:

In [13]:

i = 2 * range(nl)
print i

[0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6]

Se construyen las mastrices dispersas Yp y Yq tal que Yf*V sea el vector de corrientes inyectadas a los nodos de envio y Yq*V se el vector de corrientes inyectadas a los nodos de llegada. La función np.r_ concatena las dos matrices dadas.

In [14]:

Yp = spmatrix(np.r_[Ypp, Ypq], i, np.r_[NP, NQ], (nl, m), 'z')
Yq = spmatrix(np.r_[Yqp, Yqq], i, np.r_[NP, NQ], (nl, m), 'z')
print Yp
print Yq

[ 5.00-j14.97 -5.00+j15.00      0            0            0      ]
[  1.25-j3.73      0        -1.25+j3.75      0            0      ]
[     0         1.67-j4.98  -1.67+j5.00      0            0      ]
[     0         1.67-j4.98      0        -1.67+j5.00      0      ]
[     0         2.50-j7.49      0            0        -2.50+j7.50]
[     0            0         0.01-j1.07  -0.01+j1.08      0      ]
[     0            0            0         1.25-j3.73  -1.25+j3.75]

[-5.00+j15.00  5.00-j14.97      0            0            0      ]
[ -1.25+j3.75      0         1.25-j3.73      0            0      ]
[     0        -1.67+j5.00   1.67-j4.98      0            0      ]
[     0        -1.67+j5.00      0         1.67-j4.98      0      ]
[     0        -2.50+j7.50      0            0         2.50-j7.49]
[     0            0        -0.01+j1.08   0.01-j1.07      0      ]
[     0            0            0        -1.25+j3.75   1.25-j3.73]

Se agrupa la matriz Ybus multiplicando cada matriz por su matriz de conectividad y añadiendo el efecto de las susceptacias y conductancias nodales:

In [15]:

Ybus = Cp.T * Yp + Cq.T * Yq + spmatrix(Ysh, range(m), range(m), (m, m), 'z')
print Ybus

[  6.25-j18.70  -5.00+j15.00   -1.25+j3.75       0             0      ]
[ -5.00+j15.00  10.83-j32.41   -1.67+j5.00   -1.67+j5.00   -2.50+j7.50]
[  -1.25+j3.75   -1.67+j5.00    2.93-j9.77   -0.01+j1.08       0      ]
[      0         -1.67+j5.00   -0.01+j1.08    2.93-j9.77   -1.25+j3.75]
[      0         -2.50+j7.50       0         -1.25+j3.75   3.75-j11.21]

Notar que para este caso la matriz no es tan dispersa (sólo 6 elementos nulos). Sin embargo en sistemas más grandes las matrices son altamente dispersas. Con esto queda finalizado un breve tutorial de como formar la matriz Ybus sin el uso de ciclos, utilizando transformaciones lineales y sacando ventaja a las matrices dispersas. Autor: Uriel Fernando Sandoval