Estimaci
´
on de Sincrofasores en Redes El
´
ectricas
Inteligentes: de Modelos a Restricciones de
Dise
˜
no
Francisco Messina
†1
, Leonardo Rey Vega
†∗2
y Cecilia G. Galarza
†∗3
†
Facultad de Ingenier
´
ıa, Universidad de Buenos Aires
1,2,3
{fmessina,lrey,cgalar}@fi.uba.ar
∗
Centro de Simulaci
´
on Computacional, CONICET
Abstract—We first present a detailed review of two popular
classes of model based synchrophasor estimators. They are
the discrete Fourier transform (DFT) and the Taylor-Fourier
transform (TFT) class estimators, which arise from constant
or polynomial phasor models, respectively. Limitations of
these techniques, arising from the models themselves, are
exposed. These limitations make the original DFT and TFT
approaches unfitted for situations where the phasor can not
be properly described in such simple terms. This issue can be
addressed with the recently introduced convex semi-infinite
programming (CSIP) based estimators, which are also
reviewed. In particular, we present the constraints associated
with the tests which are most critical for the model based
methods, showing how to precisely control the estimators
performance in those cases. Finally, we show that the DFT
and TFT estimators are particular instances of the CSIP
estimator, so that there exists a connection between these
two apparently different approaches. This opens the door for
future analyses and developments.
Resumen— Presentamos en primer lugar una rese
˜
na
de dos populares clases de m
´
etodos de estimaci
´
on de
sincrofasores. Ellas son la clase de m
´
etodos basados en la
transformada de Fourier discreta (DFT) y la transformada
de Taylor-Fourier (TFT), que surgen de un modelo de
fasor constante y polin
´
omico, respectivamente. Se exponen
las limitaciones de estos m
´
etodos, que surgen por los
propios modelos asumidos, lo cual los hace poco apropiados
en situaciones donde el fasor no puede representarse en
t
´
erminos tan simples. Este problema puede resolverse con los
estimadores basados en optimizaci
´
on convexa semi-infinita
(CSIP), un enfoque novedoso que tambi
´
en es descripto. En
particular, enfatizamos las restricciones asociadas con las
situaciones m
´
as cr
´
ıticas para los m
´
etodos basados en modelos,
mostrando c
´
omo controlar precisamente el desempe
˜
no de
los estimadores en dichos casos. Finalmente, mostramos que
los estimadores DFT y TFT son instancias particulares del
estimador CSIP, de modo existe una conexi
´
on entre estos dos
enfoques aparentemente diferentes. Esto abre la puerta para
futuros an
´
alisis y desarrollos.
I. INTRODUCCI
´
ON
A. Contexto
En los
´
ultimos a
˜
nos ha cobrado un gran inter
´
es, princi-
palmente por su gran impacto potencial a nivel mundial,
el paradigma de red el
´
ectrica inteligente (en ingl
´
es Smart
Grid o SG) [1]. Una SG busca alcanzar altos niveles de efi-
ciencia y confiabilidad mediante la integraci
´
on de distintos
recursos como sensores inteligentes para monitoreo, fuentes
de energ
´
ıa renovables, generaci
´
on distribuida de energ
´
ıa,
control robusto y autom
´
atico, seguridad cibern
´
etica, etc. La
confluencia de estos elementos involucra diversas
´
areas de
ingenier
´
ıa el
´
ectrica y electr
´
onica, as
´
ı como
´
areas de ciencias
de la computaci
´
on.
Las unidades de medici
´
on fasorial (PMUs) son una com-
ponente fundamental de las SGs [2] pues permiten obtener
pr
´
acticamente en tiempo real los par
´
ametros m
´
as impor-
tantes de las se
˜
nales de tensi
´
on y corriente en nodos cr
´
ıticos
de la red el
´
ectrica. Concretamente, dichos par
´
ametros son:
amplitud y fase (o conjuntamente fasor), frecuencia, y tasa
de variaci
´
on de frecuencia (ROCOF). La gran ventaja de esta
tecnolog
´
ıa, con respecto a los sistemas tradicionales, radica
en que las mediciones de los PMUs est
´
an sincronizadas al
est
´
andar de tiempo UTC, t
´
ıpicamente por medio de un re-
ceptor GPS. Esta caracter
´
ıstica es lo que permite integrar en
forma coherente mediciones de PMUs en grandes regiones
geogr
´
aficas, lo cual aporta informaci
´
on muy valiosa de la
red el
´
ectrica. De hecho, dicha informaci
´
on tiene un impacto
directo en las diversas aplicaciones de estos sistemas de
medici
´
on de gran
´
area (WAMS) [3].
Idealmente, las se
˜
nales de tensi
´
on son sinusoidales puras
de frecuencia 50 o 60 Hz (dependiendo del pa
´
ıs), pero en
la pr
´
actica existir
´
an peque
˜
nas perturbaciones con respecto
a los par
´
ametros nominales que es importante observar
peri
´
odicamente. Estas perturbaciones pueden ser, por ejem-
plo, corrimientos en la frecuencia, modulaciones en amplitud
o en fase, rampas en frecuencia, saltos en amplitud o en
fase, etc. A su vez, la se
˜
nal puede estar contaminada por
arm
´
onicas, interarm
´
onicas y ruido. Otro requisito importante
es que la latencia de las mediciones debe ser suficientemente
peque
˜
na, lo cual es particularmente cr
´
ıtico en aplicaciones
de control y protecci
´
on. Por esta raz
´
on, las mediciones de
los PMUs deben satisfacer un conjunto de especificaciones
estrictas, que est
´
an establecidas en el est
´
andar de sincrofa-
sores IEEE Std. C37.118.1-2011 [4] y su posterior enmienda
IEEE C37.118.1a-2014 [5]. En adelante, por conveniencia,
nos referiremos a estos est
´
andares simplemente como el Std.
B. Definici
´
on de Sincrofasor, Frecuencia y ROCOF
La representaci
´
on fasorial de se
˜
nales senoidales es muy
com
´
un en el an
´
alisis de circuitos electr
´
onicos y sistemas
el
´
ectricos. Sea
x(t) = a cos(2πf
0
t + φ), (1)
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
79
Recibido: 02/08/17; Aceptado: 12/09/17
una onda senoidal, donde f
0
es la frecuencia nominal, a es
la amplitud y φ la fase. Su representaci
´
on fasorial
1
se define
como
X = a e
jφ
. (2)
Cuando el origen temporal t = 0 coincide exactamente con
el paso de un segundo en el est
´
andar de tiempo UTC, X
se denomina representaci
´
on sincrofasorial de x(t). En este
caso, φ representa la diferencia de fase con respecto a un
coseno de frecuencia f
0
sincronizado con el tiempo UTC.
En forma m
´
as general, la amplitud y la fase dependen del
tiempo: a = a(t), φ = φ(t). En este caso, la se
˜
nal de tensi
´
on
o corriente se puede escribir como
x(t) = a(t) cos[2πf
0
t + φ(t)], (3)
y su sincrofasor es
X(t) = a(t) e
jφ(t)
. (4)
Entonces, en general, los sincrofasores var
´
ıan en el tiempo.
La frecuencia y la ROCOF instant
´
aneas (con respecto a
los valores nominales f
0
y 0, respectivamente) se definen
a partir de la primera y segunda derivadas de la fase:
f(t) = f
0
+
1
2π
dφ(t)
dt
, (5)
ROCOF(t) =
1
2π
d
2
φ(t)
dt
2
. (6)
C. Definici
´
on de M
´
etricas de Desempe
˜
no
El Std. define tres m
´
etricas de desempe
˜
no que se utilizan
para evaluar todos los tests
2
: error vectorial total (TVE),
error en frecuencia (FE) y error en tasa de variaci
´
on de
frecuencia (RFE). El TVE se define del siguiente modo:
TVE(t) =
|
b
X(t) − X(t)|
|X(t)|
, (7)
donde X(t) es el sincrofasor en el tiempo t y
b
X(t) es su
estimaci
´
on para dicho instante de tiempo. Es decir, no es
otra cosa que el m
´
odulo del error relativo. Gr
´
aficamente,
se puede interpretar como la distancia entre X(t) y
b
X(t)
normalizada por |X(t)|. Generalmente, para estar en norma
el fasor estimado debe ser tal que el TVE sea menor a alg
´
un
valor, por ejemplo, el 1 %. Notar que cuanto mayor es el
error en amplitud, menor es la tolerancia en el error de fase
y viceversa. El FE se define como el error absoluto entre
la frecuencia instant
´
anea f(t) de la se
˜
nal y su estimaci
´
on
b
f(t):
FE(t) = |f(t) −
b
f(t)|. (8)
El RFE se define como el error absoluto entre la ROCOF
instant
´
anea f
0
(t) de la se
˜
nal y su estimaci
´
on
b
f
0
(t):
RFE(t) = |f
0
(t) −
b
f
0
(t)| (9)
1
Muchas veces se trabaja con el fasor RMS definido como X
RMS
=
X/
√
2. Sin embargo, aqu
´
ı preferimos por simplicidad trabajar con el fasor
X. Obviamente siempre se puede obtener una representaci
´
on a partir de la
otra con un simple cambio de escala.
2
La excepci
´
on es el test de saltos en amplitud y fase, donde se definen
otras m
´
etricas de desempe
˜
no para evaluar las respuestas transitorias. Ellas
son: el sobrepico, los tiempos de respuesta del TVE, FE y RFE, y el tiempo
de retardo.
D. Clases de PMUs
El Std. define dos clases de desempe
˜
no para el PMU: clase
M y clase P. La clase M, o de Monitoreo, es la clase m
´
as
adecuada para aplicaciones que requieren mediciones de alta
exactitud y admiten una latencia considerable en el reporte
de las mediciones. En cambio, la clase P, o de Protecci
´
on,
se utiliza en aplicaciones que exigen respuestas r
´
apidas pero
de menor exactitud. Como es de esperarse, los requisitos
de respuesta r
´
apida y alta exactitud constituyen uno de los
compromisos m
´
as importantes en esta aplicaci
´
on.
E. Reporte de Estimaciones
El PMU reporta los datos a un concentrador de datos para
su posterior uso con una tasa dada. Cabe notar que dicha tasa
no est
´
a relacionada con la tasa de muestreo del conversor
AD dentro del PMU, sino que de acuerdo al Std. debe ser
igual a un subm
´
ultiplo de la frecuencia nominal del sistema.
Para f
0
= 50 Hz, se exige que la tasa de reporte F
s
pueda
tomar los valores 10, 25 y 50 reportes por segundo (fps).
Cada tasa de reporte representa un modo de funcionamiento
diferente del dispositivo pues tiene distintos requisitos.
F. Consideraciones de Latencia
La latencia del PMU se define como la diferencia tem-
poral entre el instante en que los datos se encuentran
disponibles a la salida del PMU y el instante exacto del
tiempo de reporte indicado por la etiqueta temporal. Los
l
´
ımites m
´
aximos dependen principalmente de la tasa de
reporte y la clase del PMU, ya que estos factores determinan
el tipo de filtrado que se realiza. Concretamente, los l
´
ımites
de latencia son 2/F
s
y 7/F
s
para el desempe
˜
no clase P y
clase M respectivamente.
G. Motivaci
´
on y contribuci
´
on del art
´
ıculo
Los m
´
etodos m
´
as importantes de estimaci
´
on de sincrofa-
sores, basados en modelos de fasor constante o polin
´
omico,
tienen claras limitaciones impuestas por los mismos mod-
elos. Luego de disctuir esto, presentamos un enfoque muy
general introducido recientemente basado en t
´
ecnicas de op-
timizaci
´
on convexa semi-infinita
3
y mostramos la conexi
´
on
que tiene con los m
´
etodos previos. Cabe mencionar que
existen otros art
´
ıculos que presentan una comparaci
´
on de
distintos m
´
etodos [6]–[8]. La contribuci
´
on principal de este
art
´
ıculo es ampliar la perspectiva del estado del arte de
la estimaci
´
on de sincrofasores y enriquecerla con nuestras
contribuciones en el
´
area.
H. Organizaci
´
on del art
´
ıculo
El resto del art
´
ıculo est
´
a organizado del siguiente modo.
En la Secci
´
on II presentamos algunas consideraciones gen-
erales, comunes a todos los algoritmos. Luego, en la Secci
´
on
III discutimos en detalle los m
´
etodos basados en la DFT,
y luego mostramos sus limitaciones. El mismo an
´
alisis se
realiza para los m
´
etodos basados en la TFT en la Secci
´
on IV.
Posteriormente, en la Secci
´
on V se presentan los m
´
etodos de
dise
˜
no de filtros basados en CSIP, as
´
ı como su relaci
´
on con
los m
´
etodos previos, mostrando la generalidad del enfoque
de optimizaci
´
on. Finalmente, concluimos el art
´
ıculo en la
Secci
´
on VI.
3
Puede pensarse que el m
´
etodo basado en optimizaci
´
on se construye a
partir de una familia infinita de modelos, como veremos en la Secci
´
on V.
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
80
http://elektron.fi.uba.ar
II. PRELIMINARES
A. Preprocesamiento
Todos los algoritmos que vamos a analizar en este trabajo
parten de mediciones que se tienen de una se
˜
nal trif
´
asica de
tensi
´
on o corriente, que es digitalizada mediante un conver-
sor AD, produciendo una se
˜
nal vectorial que denotaremos
x
abc
[n] = [x
a
[n], x
b
[n], x
c
[n]]
T
, donde a, b, c representan
cada una de las fases del sistema el
´
ectrico, y que se puede
descomponer del siguiente modo
x
abc
[n] = s
abc
[n] + n
abc
[n], (10)
donde s
abc
[n] es la componente fundamental de inter
´
es
y n
abc
[n] representa interferencias (arm
´
onicas e inter-
arm
´
onicas) y ruido. En general, s
abc
[n] puede representarse
como la suma de las componentes sim
´
etricas instant
´
aneas
[9], tal como se muestra en la Ec. (11) (ver parte in-
ferior de la siguiente p
´
agina), donde a
i
[n] ≡ a
i
(nT ) y
φ
i
[n] ≡ φ
i
(nT ), con i = a, b, c, 0, 1, 2, son la amplitud
y fase instant
´
aneas de cada fase (a, b, c) o componente (0,1
o positiva,2 o negativa) correspondiente, ω
0
es la frecuencia
angular nominal del sistema de potencia, T = 2π/(Cω
0
) es
el per
´
ıodo de muestreo, siendo C la cantidad de muestras
que se obtienen por ciclo nominal, y ν
0
= ω
0
T = 2π/C.
Sin embargo, los algoritmos no operan directamente con
la se
˜
nal x
abc
[n], sino que primero la transforman con el
fin de construir mediciones directas del fasor de secuencia
positiva de la componente fundamental, que posteriormente
ser
´
an procesadas para obtener las estimaciones del sincrofa-
sor, la frecuencia y la ROCOF. Como esta etapa es general,
independiente de la elecci
´
on del algoritmo de estimaci
´
on,
se presenta en esta secci
´
on. Dicha conversi
´
on se basa en las
transformaciones
4
de Clarke (o transformaci
´
on αβγ) y de
Park (o transformaci
´
on abc−dqo), que son transformaciones
de coordenadas frecuentemente utilizadas en el an
´
alisis de
sistemas de potencia trif
´
asicos [10]. Ambas son transforma-
ciones lineales. La transformaci
´
on de Clarke queda definida
por la matriz constante
C =
2
3
1 −
1
2
−
1
2
0
√
3
2
−
√
3
2
, (12)
mientras que la transformaci
´
on de Park se puede representar
por la siguiente matriz variante en el tiempo
P [n]
=
2
3
cos(φ
a,0
[n])
cos(φ
b,0
[n]) cos(φ
c,0
[n])
− sin(φ
a,0
[n]) − sin(φ
b,0
[n]) − sin(φ
c,0
[n])
,
(13)
donde φ
a,0
[n] = ν
0
n, φ
b,0
[n] = ν
0
n − 2π/3 y φ
c,0
[n] =
ν
0
n + 2π/3. Operando, obtenemos que
x
dq
[n] =
x
d
[n]
x
q
[n]
= Cx
abc
[n] = a
1
[n]
cos(ν
0
n + φ
1
[n])
sin(ν
0
n + φ
1
[n])
+ a
2
[n]
cos(−ν
0
n − φ
2
[n])
sin(−ν
0
n − φ
2
[n])
+ Cn
abc
[n], (14)
4
Estrictamente, usamos las transformaciones de Clark y Park reducidas
pues tanto la componente γ como la componente cero (o) no contienen
informaci
´
on del sincrofasor y, por lo tanto, no son calculadas.
y
dq
[n] =
y
d
[n]
y
q
[n]
= P [n]x
abc
[n] = a
1
[n]
cos(φ
1
[n])
sin(φ
1
[n])
+ a
2
[n]
cos(−2ν
0
n − φ
2
[n])
sin(−2ν
0
n − φ
2
[n])
+ P [n]n
abc
[n].
(15)
Es f
´
acil ver que hay una simple relaci
´
on entre ambas trans-
formaciones, a saber, P [n] = R(−ν
0
n)C, donde R(−ν
0
n)
es una matriz de rotaci
´
on en un
´
angulo −ν
0
n en sentido
antihorario. En general, los algoritmos se formulan a partir
de la se
˜
nal x
d
[n]+jx
q
[n], pero equivalentemente se pueden
formular a partir de la se
˜
nal y
d
[n] + jy
q
[n], lo cual muchas
veces ofrece una versi
´
on simplificada del m
´
etodo en cuesti
´
on
y revela m
´
as claramente la esencia del mismo.
Observando la ecuaci
´
on (15), es claro que el primer
t
´
ermino se puede identificar como el sincrofasor X[n],
mientras que los otros t
´
erminos representan el efecto de la
componente negativa, interferencias y ruido. Cabe destacar
algunas ventajas de esta transformaci
´
on [11]:
1) A diferencia de una demodulaci
´
on monof
´
asica, no
genera un t
´
ermino de doble frecuencia para se
˜
nales
de entrada balanceadas, es decir, con a
0
[n] = 0 y
a
2
[n] = 0.
2) Filtra por completo la componente cero que, como
mencionamos previamente, no contiene informaci
´
on
´
util para nuestro problema.
3) La frecuencia de la componente positiva es mapeada
a 0 rad./muestra, mientras que la componente neg-
ativa es trasladada a −2ν
0
rad./muestra. Este de-
sacoplamiento en frecuencia tiene el efecto de pro-
ducir una mayor robustez ante desbalanceos de la
se
˜
nal de entrada.
B. Algoritmos de procesamiento por bloques
Los m
´
etodos basados en procesamiento por bloques
parten de un conjunto finito de datos a partir del cual se
deben estimar los par
´
ametros de inter
´
es. Concretamente,
la formulaci
´
on general es la siguiente. A partir de N
muestras de una se
˜
nal de tensi
´
on o corriente se pretende
estimar mediante alguna operaci
´
on matem
´
atica el fasor,
la frecuencia, y la ROCOF correspondientes al tiempo de
reporte. Por simplicidad, supondremos que el n
´
umero de
muestras en un bloque es impar. Entonces, denotaremos por
x = {x
n
}
n∈N
= [x
−R
, . . . , x
0
, . . . , x
R
]
T
∈ C
N
al vector
de muestras en un bloque de la se
˜
nal x
d
[n] + jx
q
[n] y por
y = {y
n
}
n∈N
= [y
−R
, . . . , y
0
, . . . , y
R
]
T
∈ C
N
al vector
de muestras del mismo bloque de la se
˜
nal y
d
[n] + jy
q
[n],
donde R =
N−1
2
y N = {−R, . . . , R}. La relaci
´
on entre
las componentes de ambos vectores es simplemente
y
n
= x
n
e
−jν
0
n
, n ∈ N . (16)
Adem
´
as, asumimos que el tiempo de reporte de las medi-
ciones siempre coincide con el instante de tiempo asociado
a la muestra x
0
o y
0
, es decir, con el centro del bloque.
Este planteo es habitual en la literatura de sincrofasores,
posiblemente porque es el que se propone en el mismo Std.
[4]. Notar que el hecho de que el tiempo de reporte est
´
e en el
centro del bloque hace que la estimaci
´
on sea no causal. Esto
puede llamar la atenci
´
on del lector pero no es problem
´
atico
pues, como hemos mencionado en la Secci
´
on I-F, el Std.
admite cierta latencia en el reporte de las mediciones de las
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
81
http://elektron.fi.uba.ar
PMUs. De hecho,
´
este es un t
´
ıpico problema de suavizado
en procesamiento de se
˜
nales. Sin embargo, es evidente que
el largo del bloque queda delimitado por la latencia m
´
axima
tolerable t
L
. Concretamente, para no violar esta restricci
´
on,
debemos tener necesariamente que
5
R ≤
t
L
T
. (17)
Recordemos que la latencia m
´
axima tolerable depende de
la clase del PMU. Luego, la restricci
´
on sobre el largo del
bloque ser
´
a distinta para aplicaciones de protecci
´
on (clase
P) y monitoreo (clase M).
III. ALGORITMOS BASADOS EN LA DFT
Tradicionalmente, el algoritmo que se ha utilizado con
mayor popularidad para estimar el sincrofasor se basa en
una simple transformada de Fourier discreta (DFT) [12]. En
particular, como
´
unicamente se desea conocer el fasor de la
componente fundamental de la se
˜
nal de entrada, la operaci
´
on
que se realiza para obtener la estimaci
´
on del sincrofasor es
la siguiente:
b
X =
1
N
X
n∈N
x
n
e
−jν
0
n
. (18)
Si reescribimos
b
X en t
´
erminos de la se
˜
nal {y
n
}, tenemos
b
X =
1
N
X
n∈N
y
n
. (19)
Luego, la interpretaci
´
on de este estimador fasorial es muy
sencilla:
b
X no es otra cosa que el promedio de la se
˜
nal
{y
n
}, es decir, de la se
˜
nal {x
n
} demodulada a frecuencia
nominal. Esta observaci
´
on permite hacer algunas conexiones
interesantes que ser
´
an
´
utiles posteriormente:
1) El estimador
b
X se puede interpretar como la salida en
el instante n = 0 de un filtro promediador, es decir,
cuya respuesta impulsiva es
h
n
=
1
N
1{n ∈ N }, (20)
donde 1{·} representa la funci
´
on indicadora. Como
el filtro promediador es simplemente un filtro pasa
bajos, esto sugiere que se pueden utilizar otros filtros
pasa bajos, tal como se muestra en el ap
´
endice C del
Std. [4]. Luego, es posible utilizar t
´
ecnicas de dise
˜
no
de filtros para obtener nuevos estimadores fasoriales.
Esta idea se discute en detalle en la Secci
´
on V.
2) En t
´
erminos estad
´
ısticos, el estimador de la DFT
se puede interpretar como la media muestral de la
secuencia {y
n
}. Dicho estimador es
´
optimo en el
5
En rigor, en la pr
´
actica debemos considerar un margen de latencia
para los retardos propios del c
´
omputo de las estimaciones. Dichos retardos
dependen del hardware utilizado para la implementaci
´
on del algoritmo pero
generalmente son despreciables frente a la latencia propia del esquema de
procesamiento por bloques. Por lo tanto, su contribuci
´
on a la latencia ser
´
a
ignorada en este trabajo.
siguiente sentido. Supongamos que las mediciones se
pueden modelar como
x
n
= s
n
+ w
n
= Xe
jν
0
n
+ u
n
, n ∈ N , (21)
donde X = Ae
jφ
es el sincrofasor, {s
n
} es la se
˜
nal
y {u
n
} es ruido blanco complejo circular Gaussiano
de varianza σ
2
. Las muestras de {y
n
} son
y
n
= X + v
n
, n ∈ N , (22)
donde v
n
= e
−jν
0
n
u
n
. Por definici
´
on, {v
n
} tambi
´
en
es ruido blanco complejo circular Gaussiano. En
forma vectorial, el modelo se puede escribir como
y = X1 + v, (23)
donde X es simplemente un par
´
ametro complejo que
se desea estimar y 1 es un vector de N unos. Luego,
y ∼ CN (X1, σ
2
I). Es bien sabido en teor
´
ıa de
estimaci
´
on que el estimador
b
X = (1
T
1)
−1
1
T
y =
1
N
X
n∈N
y
n
, (24)
es eficiente, es decir, es un estimador insesgado que
satisface la cota de Cram
´
er-Rao y, por lo tanto, es
un estimador insesgado de m
´
ınima varianza [13]. Su
varianza es entonces
Var(
b
X) = F
−1
(X) =
σ
2
N
, (25)
donde F (X) es la informaci
´
on de Fisher de X en y.
3) Es interesante observar que es posible relacionar am-
bas interpretaciones anteriores mediante el concepto
de filtro adaptado, que es el filtro que maximiza
la relaci
´
on se
˜
nal a ruido a su salida y act
´
ua como
correlador de la se
˜
nal y y la se
˜
nal constante 1. Este
concepto es cl
´
asico en procesamiento de se
˜
nales [14].
Una propiedad bien conocida de la DFT es su excelente
rechazo de la componente continua y las arm
´
onicas en
condiciones de frecuencia nominal, que se puede entender
f
´
acilmente a partir de la ortogonalidad de la familia de
se
˜
nales {{e
jhν
0
n
} : h = 0, . . . , N − 1} [15] suponiendo que
K = N/C ∈ N, siendo K la cantidad de ciclos nominales
dentro del bloque. Concretamente, supongamos que la se
˜
nal
es de la forma
s
n
= Xe
jν
0
n
+ X
h
e
jhν
0
n
, n ∈ N , (26)
donde h ∈ {0, 2, . . . , N − 1} es el
´
ındice de arm
´
onica y X
h
es el fasor de la arm
´
onica. Ignorando el ruido
6
, tenemos
entonces que
b
X = X + X
h
1
N
X
n∈N
e
j(h−1)ν
0
n
= X, (27)
6
En lo que sigue, salvo cuando se discuta el comportamiento estad
´
ıstico
de un estimador, siempre ignoraremos el efecto del ruido.
s
abc
[n]
=
a
a
[n] cos(ν
0
n + φ
a
[n])
a
b
[n] cos(ν
0
n + φ
b
[n])
a
c
[n] cos(ν
0
n + φ
c
[n])
= a
0
[n]
cos(ν
0
n + φ
0
[n])
cos(ν
0
n + φ
0
[n])
cos(ν
0
n + φ
0
[n])
+a
1
[n]
cos(ν
0
n + φ
1
[n])
cos(ν
0
n + φ
1
[n] −
2π
3
)
cos(ν
0
n + φ
1
[n]
+
2π
3
)
+a
2
[n]
cos(ν
0
n + φ
2
[n])
cos(ν
0
n + φ
2
[n]
+
2π
3
)
cos(ν
0
n + φ
2
[n] −
2π
3
)
,
(11)
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
82
http://elektron.fi.uba.ar
pues
P
n∈N
e
j(h−1)ν
0
n
= he
jν
0
n
, e
jhν
0
n
i = 0, donde h·, ·i
es el producto interno can
´
onico en C
N
.
A pesar de ser un estimador muy interesante por todas las
propiedades que hemos enunciado previamente, el estimador
de la DFT tiene un comportamiento pobre en condiciones
de frecuencia no nominal [16], lo cual muestra que no es
robusto ante peque
˜
nas discrepancias entre el modelo de se
˜
nal
y la se
˜
nal real. Para comprender este problema, supongamos
ahora que la se
˜
nal s
n
no tiene frecuencia nominal sino que
existe una peque
˜
na desviaci
´
on δ de modo que
s
n
= Xe
j(ν
0
+δ)n
, n ∈ N . (28)
Entonces, tendremos
b
X =
D
R
(δ)
N
X, (29)
donde D
R
(δ) es el kernel de Dirichlet que se puede expresar
como
D
R
(δ) =
X
n∈N
e
jδn
=
sin(
Nδ
2
)
sin(
δ
2
)
. (30)
Luego, vemos que el estimador del sincrofasor para frecuen-
cias no nominales se ve afectado por la ganancia D
R
(δ)/N.
Entonces, podemos observar que el TVE resulta
TVE =
|
b
X − X|
|X|
=
D
R
(δ)
N
− 1
. (31)
Mediante una expansi
´
on en serie de Taylor alrededor de
δ = 0, es f
´
acil ver que
D
R
(δ)
N
− 1 =
1 − N
2
24
δ
2
+ O(N
4
δ
4
), (32)
de modo que el TVE crece aproximadamente en forma
cuadr
´
atica con el corrimiento en frecuencia δ y el largo del
bloque N. Expresiones an
´
alogas para una se
˜
nal monof
´
asica
se pueden encontrar en [16]. Se observa que el compor-
tamiento de la DFT es a
´
un peor en dicho caso por la
presencia de la componente imagen. Del mismo modo,
podemos observar que el rechazo de arm
´
onicas tambi
´
en
sufre en condiciones de frecuencia no nominal pues la
ortogonalidad de las exponenciales complejas deja de ser
v
´
alida.
En la Fig. 1 presentamos un gr
´
afico del TVE en funci
´
on
de ∆f = δ/(2πT ) para bloques de distinto largo, marcando
el l
´
ımite del %1 establecido en el Std. considerando T = 1
ms. Tambi
´
en se muestran los puntos ∆f = ±2 Hz que
se corresponden con los l
´
ımites de frecuencia no nominal
para el PMU clase P. Se puede observar que en todos los
casos el TVE supera la cota del Std. cuando se considera el
rango de frecuencias [45, 55] Hz, que se corresponde con el
PMU clase M. Solamente para el bloque m
´
as corto (de un
ciclo de largo), el TVE satisface la norma para el rango de
frecuencias del PMU clase P.
Una generalizaci
´
on natural de la DFT para reducir estos
efectos de leakage, es utilizar una DFT ventaneada o WDFT.
En dicho caso, el estimador se puede escribir como
b
X =
X
n∈N
w
n
y
n
X
n∈N
w
n
, (33)
Fig. 1: TVE del estimador fasorial basado en la DFT en
funci
´
on del corrimiento en frecuencia ∆f para bloques de
distinto largo.
donde {w
n
} son los coeficientes de la ventana elegida. Es
claro que con esta generalizaci
´
on,
b
X es simplemente un
promedio ponderado de la se
˜
nal {y
n
}. El estudio de la influ-
encia de ventanas en aplicaciones de estimaci
´
on espectral se
ha estudiado ampliamente [17]–[19]. En particular, una clase
rica de ventanas que se ha propuesto utilizar en la aplicaci
´
on
de sincrofasores es la clase cosenoidal [20], que incluye a las
ventanas tradicionales de Hanning, Hamming, Blackman y
Nuttall. Otra ventana atractiva, por ser una aproximaci
´
on de
la ventana
´
optima
7
de Slepian [15], es la ventana de Kaiser,
que se define como
w
n
=
I
0
β
q
1 −
n
R
2
I
0
(β)
, n ∈ N , (34)
donde I
0
(·) es la funci
´
on de Bessel modificada de orden
cero y β es un par
´
ametro no negativo que controla la forma
de la ventana. B
´
asicamente, a medida que β aumenta, el
ancho del l
´
obulo principal de la ventana aumenta y el
´
area
de los l
´
obulos secundarios disminuye. El estimador WDFT
en condiciones de frecuencia no nominal ser
´
a
b
X =
P
n∈N
w
n
Xe
jδn
P
n∈N
w
n
= X
W (−δ)
W (0)
, (35)
donde W (ν) =
P
n∈N
w
n
e
−jνn
es la transformada de
Fourier (DTFT) de {w
n
}. Luego, si la ventana satisface la
condici
´
on
W (−δ)
W (0)
≤
D
R
(δ)
N
, (36)
para un rango apropiado de δ, se tendr
´
a un mejor compor-
tamiento que con el estimador DFT.
Otro algoritmo muy popular para compensar los errores
que aparecen en condiciones de frecuencia no nominal es el
algoritmo IpDFT o DFT interpolada [21]. La motivaci
´
on
del algoritmo es muy sencilla. De la Ec. (35) podemos
apreciar que si conoci
´
eramos el corrimiento en frecuencia δ,
podr
´
ıamos compensar el factor que aparece en el estimador
WDFT perfectamente. En el algoritmo IpDFT, se comienza
por evaluar la WDFT de la secuencia {y
n
} en la frecuencia
nominal y en los bins adyacentes:
Z(0) = XW (−δ), Z(±ε) = XW (±ε − δ), (37)
7
En el sentido de maximizar la relaci
´
on energ
´
ıa del l
´
obulo principal /
energ
´
ıa total de la ventana.
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
83
http://elektron.fi.uba.ar
donde Z(ν) es la DTFT de la secuencia {z
n
}, con z
n
=
w
n
y
n
para todo n ∈ N , y ε = 2π/N . Luego, se comparan
los valores de |Z(ε)| y |Z(−ε)|. Si |Z(ε)| > |Z(−ε)|,
tendremos que δ > 0, y en caso contrario δ < 0. Esto asume
que el corrimiento en frecuencia satisface |δ| < ε/2 [18], lo
cual en general es cierto para la aplicaci
´
on de sincrofasores,
pero puede no valer para bloques muy largos. Supongamos,
por ejemplo, que se satisface la condici
´
on |Z(ε)| > |Z(−ε)|.
Luego, el siguiente paso es evaluar
γ
.
=
|Z(ε)|
|Z(0)|
=
W (−δ + ε)
W (−δ)
= f(δ), (38)
donde f (·) es una funci
´
on que depende de la ventana
elegida. Cabe notar que la igualdad γ = f(δ) es v
´
alida
´
unicamente para el caso ideal en que no hay ruido ni inter-
ferencia. Finalmente, se busca invertir la relaci
´
on anterior
para obtener una estimaci
´
on del corrimiento en frecuencia
b
δ = f
−1
(γ), (39)
que luego se utilizar
´
a para compensar el estimador WDFT
de la Ec. (35). Desde un punto de vista pr
´
actico, a
´
un si f(·)
es una funci
´
on muy compleja para invertir anal
´
ıticamente,
una aproximaci
´
on polin
´
omica de f
−1
(·) puede ser sufi-
cientemente exacta [22]. Cabe notar que por medio del
algoritmo IpDFT, obtenemos adem
´
as de la correcci
´
on del
leakage, un estimador de frecuencia. Alternativamente, para
los estimadores DFT y WDFT, se puede utilizar el estimador
que se propone en el Std. basado en diferencias finitas [5].
Del mismo modo, ninguno de los m
´
etodos basados en la
DFT proporciona un estimador de ROCOF, de modo que
es razonable utilizar nuevamente el m
´
etodo propuesto en el
Std. Cabe notar que dichos m
´
etodos requieren la estimaci
´
on
de dos sincrofasores adicionales
8
por cada tiempo de reporte
que no son reportados pero son necesarios para obtener esti-
maciones de frecuencia y ROCOF suficientemente exactas.
Adem
´
as, como bien se se
˜
nala en el Std., los estimadores
propuestos son muy sensibles al ruido.
IV. ALGORITMOS BASADOS EN EL MODELO DE
TAYLOR-FOURIER
Una limitaci
´
on fundamental de los m
´
etodos basados en
la DFT es que asumen, al menos en forma impl
´
ıcita, un
modelo estacionario. Por lo tanto, su desempe
˜
no en general
es pobre en condiciones din
´
amicas, es decir, cuando el
sincrofasor var
´
ıa en el tiempo significativamente dentro del
bloque de datos [8]. En este sentido, una contribuci
´
on muy
importante en el
´
area de estimaci
´
on de sincrofasores fue
el reconocimiento de la gran ventaja de usar un modelo
de fasor din
´
amico [23]. Concretamente, en dicho trabajo se
plantea un modelo polin
´
omico para el sincrofasor, es decir,
la se
˜
nal {x
n
} toma la siguiente forma
x
n
=
L
X
l=0
X
l
n
l
e
jν
0
n
+ u
n
, n ∈ N , (40)
donde X
0
, . . . , X
L
∈ C son los par
´
ametros del modelo.
El estimador que presentamos a continuaci
´
on se conoce
8
Si el sincrofasor que se desea estimar es X[n], la implementaci
´
on de
los estimadores de frecuencia y ROCOF por diferencias finitas requiere la
estimaci
´
on de X[n − 1] y X[n + 1].
como estimador de Taylor-Fourier (TFT) pues se basa en
una expansi
´
on en Taylor del fasor din
´
amico referido a la
frecuencia fundamental ν
0
. En t
´
erminos de la se
˜
nal {y
n
},
tenemos que
y
n
=
L
X
l=0
X
l
n
l
+ v
n
, n ∈ N . (41)
En forma vectorial podemos escribir entonces
y = N X + v, (42)
donde hemos definido la matriz
N =
1 −R (−R)
2
· · · (−R)
L
1 −R + 1 (−R + 1)
2
· · · (−R + 1)
L
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 R R
2
· · · R
L
,
(43)
y el vector de par
´
ametros X = [X
0
, . . . X
L
]
T
. Notar que el
modelo nuevamente es lineal en los par
´
ametros y, adem
´
as,
y ∼ CN (N X, σ
2
I). Luego, el estimador de m
´
ınimos
cuadrados [23] es eficiente [13]:
b
X = (N
T
N )
−1
N
T
y = Gy, (44)
donde G = (N
T
N )
−1
N
T
. Por lo tanto, su matriz de
covarianza es
Cov(
b
X) = F
−1
(X) = σ
2
GG
T
= σ
2
(N
T
N )
−1
, (45)
donde F (X) es la matriz de informaci
´
on de Fisher de X
en y. Un resultado fundamental de teor
´
ıa de estimaci
´
on
establece que la varianza del estimador de un par
´
ametro
no puede disminuir, y en general crece, cuando se estima
junto con otros par
´
ametros [13]. En nuestro caso, esto
significa que la varianza del estimador TFT del sincrofasor
(el par
´
ametro X
0
) es mayor que la varianza del estimador
DFT (el par
´
ametro X). Del mismo modo, cuanto mayor sea
L mejor ser
´
a la aproximaci
´
on polin
´
omica del sincrofasor en
general pero mayor ser
´
a la varianza por el hecho de que
el modelo tiene m
´
as par
´
ametros. Esto es simplemente una
manifestaci
´
on de la relaci
´
on de compromiso que existe entre
el sesgo y la varianza de un estimador. Concretamente, el
sesgo decrece con L pero la varianza crece con L. El caso
L = 0, que se corresponde con el estimador DFT cl
´
asico da
el estimador de m
´
ınima varianza pero con un sesgo m
´
aximo
para un fasor din
´
amico. Podemos decir que, en dicho caso, el
modelo es muy simple o existe underfitting. Por el contrario,
cuando se utiliza L > 0, se obtiene una mejor aproximaci
´
on
al fasor din
´
amico a expensas de un incremento en la var-
ianza del estimador. Eventualmente, con L suficientemente
grande, la aproximaci
´
on ser
´
a excelente, pero a expensas de
un incremento considerable en la varianza. En este caso, el
modelo es demasiado complejo o existe overfitting. Adem
´
as,
siempre se debe tener L N para evitar que la matriz G
est
´
e mal condicionada [24].
Es interesante observar que a partir de la estimaci
´
on de
los par
´
ametros X
0
, X
1
y X
2
se pueden obtener directamente
estimaciones de la frecuencia angular ω y la ROCOF angular
α, tambi
´
en en el centro del bloque de datos, explotando las
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
84
http://elektron.fi.uba.ar
relaciones entre las derivadas de un fasor din
´
amico y su
frecuencia y ROCOF [23], [25]. Las relaciones son:
bω =
1
T
Im
b
X
1
b
X
0
!
, (46)
bα =
2
T
2
"
Im
b
X
2
b
X
0
!
− Re
b
X
1
b
X
0
!
Im
b
X
1
b
X
0
!#
. (47)
Una generalizaci
´
on directa del estimador TF, que se presenta
en [24], se basa en resolver un problema de cuadrados
m
´
ınimos ponderado (WLS) en lugar de un problema de
LS. La motivaci
´
on es an
´
aloga a la que justifica utilizar
un estimador WDFT. Es decir, el estimador del vector de
par
´
ametros X surge de resolver
b
X = (N
T
W
T
W N )
−1
N
T
W
T
W y = Gy, (48)
con W = diag({w
n
}), donde {w
n
} son los pesos o
coeficientes de la ventana elegida. Por razones obvias, el
estimador de la Ec. (48) se denomina estimador WTF. Por
supuesto, la ecuaci
´
on (48) se reduce a la ecuaci
´
on (44)
cuando se toma una ventana rectangular, es decir, W = I .
En [24], los autores proponen trabajar con la ventana de
Kaiser, que hemos introducido en la Ec. (34).
Es interesante analizar, al igual que para el estimador
basado en la DFT, c
´
omo afecta un corrimiento en frecuencia
a la estimaci
´
on del sincrofasor. Nuevamente, ignoramos el
ruido y partimos del modelo de se
˜
nal presentado en la
(28). Entonces, el problema se reduce a hallar la mejor
aproximaci
´
on polin
´
omica, en el sentido LS o WLS, de una
se
˜
nal exponencial compleja que se puede escribir como
y = Xe(δ), (49)
donde e(δ) = [e
−jδR
, . . . , 1, . . . , e
jδR
]
T
. La soluci
´
on, como
ya vimos, est
´
a dada en general por la ecuaci
´
on (48). Es f
´
acil
verificar que el TVE y el FE entonces est
´
an dados por
TVE = |g
T
0
e(δ) − 1|, (50)
FE =
1
2πT
Im
g
T
1
e(δ)
g
T
0
e(δ)
− δ
, (51)
donde g
T
i
es la fila i-
´
esima de la matriz G con i = 0, . . . , L.
Adem
´
as, se puede demostrar que en este caso el RFE es nulo
pues los estimadores WTFT se pueden interpretar como un
banco de filtros diferenciador de fase lineal [24], [26]. En
la Figs. 2 y 3 mostramos c
´
omo afecta el corrimiento en
frecuencia ∆f = δ/(2πT) a la estimaci
´
on del sincrofasor
y la frecuencia mediante el estimador WTFT utilizando
distintas ventanas de Kaiser, es decir, variando el par
´
ametro
β. En las simulaciones se utiliza L = 3 tal como se propone
en [24]. Se puede observar que el comportamiento del
estimador mejora al incrementar β, es decir, cuando el ancho
de banda del espectro de la la ventana es suficientemente
grande. Adem
´
as, el resultado es mejor para bloques de datos
cortos, lo cual es esperable por la naturaleza local de la
aproximaci
´
on de Taylor. Sin embargo, cabe notar que esto
tambi
´
en significa que el estimador WTFT ser
´
a m
´
as sensible
al ruido y las se
˜
nales de interferencia. En [24], posiblemente
por esta consideraci
´
on, se adopta el valor β = 8 para bloques
de largo K = 4.
Es interesante preguntarse si es posible usar directamente
L = 2, pues para los par
´
ametros que nos interesa estimar
∆f [Hz]
-5 0 5
TVE [%]
0
0.2
0.4
0.6
0.8
1
K = 1
K = 2
K = 3
K = 4
(a) β = 0
∆f [Hz]
-5 0 5
TVE [%]
0
0.2
0.4
0.6
0.8
1
K = 1
K = 2
K = 3
K = 4
(b) β = 2
∆f [Hz]
-5 0 5
TVE [%]
0
0.2
0.4
0.6
0.8
1
K = 1
K = 2
K = 3
K = 4
(c) β = 4
∆f [Hz]
-5 0 5
TVE [%]
0
0.2
0.4
0.6
0.8
1
K = 1
K = 2
K = 3
K = 4
(d) β = 6
∆f [Hz]
-5 0 5
TVE [%]
0
0.2
0.4
0.6
0.8
1
K = 1
K = 2
K = 3
K = 4
(e) β = 8
∆f [Hz]
-5 0 5
TVE [%]
0
0.2
0.4
0.6
0.8
1
K = 1
K = 2
K = 3
K = 4
(f) β = 10
Fig. 2: TVE para el estimador de Taylor-Fourier con L = 3
para distintas ventanas de Kaiser.
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(a) β = 0
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(b) β = 2
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(c) β = 4
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(d) β = 6
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(e) β = 8
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(f) β = 10
Fig. 3: FE para el estimador de Taylor-Fourier con L = 3
para distintas ventanas de Kaiser.
esto es suficiente, tal como muestran las Ecs. (46) y (47).
Esta pregunta tiene sentido desde un punto de estad
´
ıstico
pues, como se mencion
´
o antes, en general el agregado de
par
´
ametros genera un deterioro en la varianza del estimador.
Por lo tanto, tiene sentido intentar usar el valor de L
m
´
ınimo compatible con nuestro problema. Desafortunada-
mente, como puede verse en la Fig. 4, al utilizar L = 2 el
FE crece considerablemente y excede el l
´
ımite impuesto por
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
85
http://elektron.fi.uba.ar
el Std., especialmente si se consideran las restricciones de
un PMU clase M. Adem
´
as, este comportamiento no mejora
demasiado usando ventanas de banda ancha. Esto muestra
que un polinomio cuadr
´
atico es un modelo demasiado
simple como para describir al sincrofasor en condiciones
de frecuencia no nominal, lo cual se debe simplemente al
hecho de que la aproximaci
´
on de m
´
ınimos cuadrados de
la se
˜
nal {e
jδn
} con un polinomio de grado L = 2 no es
suficientemente buena. Por lo tanto, concluimos que L = 3
representa efectivamente una buena relaci
´
on de compromiso
entre el sesgo y la varianza del estimador, o bien, la exactitud
y la complejidad del modelo.
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(a) β = 0
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(b) β = 2
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(c) β = 4
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(d) β = 6
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(e) β = 8
∆f [Hz]
-5 0 5
FE [Hz]
×10
-3
0
1
2
3
4
5
K = 1
K = 2
K = 3
K = 4
(f) β = 10
Fig. 4: FE para el estimador de Taylor-Fourier con L = 2
para distintas ventanas de Kaiser.
Existen otras variaciones y extensiones interesantes del
estimador TFT que presentamos brevemente a continuaci
´
on.
En el art
´
ıculo [27] se ha propuesto referenciar el modelo de
Taylor-Fourier a una frecuencia distinta a la nominal, que
debe ser estimada previamente, para mejorar su compor-
tamiento en los casos en que la frecuencia difiere notable-
mente de la frecuencia nominal. En particular, se propone
realizar una estimaci
´
on en dos etapas. En la primera etapa,
se estima la frecuencia fundamental de la se
˜
nal mediante el
algoritmo IpDFT, descripto en la Secci
´
on III. En la segunda
etapa, se estima el modelo de fasor din
´
amico pero referido
a la frecuencia hallada en el paso anterior en lugar de la
frecuencia nominal.
Otra generalizaci
´
on de los estimadores de Taylor-Fourier
se basa en el concepto de lazos de enganche de fase (PLL)
[28]. En dicho caso, se utiliza un modelo de la forma
y
n
=
L
X
l=0
X
l
n
l
e
jθ
n
+ v
n
, n ∈ N , (52)
donde θ
n
=
P
L
l=0
η
l,k
n
l
es un polinomio cuyos coeficientes
η
l,k
var
´
ıan con el n
´
umero de bloque k en forma recursiva.
Una posible estrategia de adaptaci
´
on que se sugiere en [28]
es utilizar los coeficientes de la fase hallados a partir de la
estimaci
´
on del vector X en el bloque (k − 1)-
´
esimo como
innovaciones para actualizar los coeficientes del polinomio
de fase en el bloque k-
´
esimo. La ventaja de este enfoque
es que permite referenciar el modelo no solamente a una
frecuencia no nominal sino tambi
´
en a derivadas de fase
superiores no nulas. Por ejemplo, para L = 2, el m
´
etodo
permite referenciar el modelo a una ROCOF no nula, lo cual
es
´
util en condiciones din
´
amicas. Luego, la aproximaci
´
on
polin
´
omica se hace solamente sobre el remanente del fasor
din
´
amico que no se puede modelar como una exponencial
compleja con fase cuadr
´
atica.
Finalmente, debemos notar que es posible extender
el modelo de Taylor-Fourier para incluir componentes
arm
´
onicas y, si se estiman previamente sus frecuencias medi-
ante alg
´
un m
´
etodo como la IpDFT, tambi
´
en interarm
´
onicas
[29]. En este caso, la se
˜
nal estar
´
a compuesta por M + 1
componentes de Taylor-Fourier centradas en las frecuencias
ν
m
, es decir
x
n
=
M
X
m=0
L
X
l=0
X
l,m
n
l
e
jν
m
n
+ u
n
, n ∈ N . (53)
Notar que para las componentes arm
´
onicas tenemos la
siguiente relaci
´
on entre las frecuencias: ν
m
= (m + 1)ν
0
.
En forma vectorial, este modelo se puede escribir como
x = AX + u, donde A =
D(ν
0
)N . . . D(ν
M
)N
.
Cuando realizamos la demodulaci
´
on de la frecuencia nom-
inal ν
0
, obtenemos el modelo de la se
˜
nal y como
y = BX + v, (54)
donde
B = D
H
(ν
0
)A =
N . . . D(ν
M
− ν
0
)N
. (55)
Nuevamente podemos definir el estimador de Taylor-Fourier
como la soluci
´
on WLS de este problema, de modo que
obtenemos
b
X = (B
H
B)
−1
B
H
X = H X, (56)
donde ahora el vector de par
´
ametros es X =
[X
0,0
, . . . , X
L,0
, X
0,1
. . . X
L,M
]
T
. Cabe notar que la
matriz H puede precomputarse solamente para el modelo
a frecuencia nominal con una cantidad prefijada de
arm
´
onicas. En dicho caso, la complejidad del algoritmo
crece linealmente con la dimensi
´
on del vector de par
´
ametros,
es decir, (M + 1)(L + 1). Si bien es posible aplicar las
generalizaciones del estimador TFT de una componente
a este caso m
´
as general, hay que notar que el costo
del algoritmo crece muy r
´
apidamente con el n
´
umero de
componentes del modelo. Por ejemplo, la matriz que
hay que invertir en la Ec. (56) es una matriz compleja de
(M+1)(L+1)×(M+1)(L+1), de modo que la complejidad
computacional
9
crece como O([(M +1)(L+1)]
3
). Adem
´
as,
como se discuti
´
o previamente, existen otro problemas a
´
un
m
´
as serios que ocurren cuando la cantidad de par
´
ametros
a estimar es similar a la cantidad de muestras del bloque
9
Una forma habitual de resolver un sistema lineal que involucra una
matriz hermitiana definida positiva es mediante la descomposici
´
on de
Cholesky.
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
86
http://elektron.fi.uba.ar
de datos. Por un lado, cuando el modelo es demasiado
complejo, tenemos un problema de overfitting. Por otro
lado, para que la matriz B
H
B est
´
e bien condicionada
debemos tener (M + 1)(L + 1) N. Esta restricci
´
on
muestra que no es factible utilizar este m
´
etodo con bloques
cortos. Una implementaci
´
on m
´
as eficiente del m
´
etodo, que
explota el hecho de que solamente se deben estimar los
primeros L + 1 par
´
ametros de X, y provee adaptaci
´
on de
la referencia de frecuencia, se presenta en [30]. A
´
un as
´
ı,
cuanto m
´
as detallado es el modelo de la se
˜
nal (determinado
por el grado de los polinomios de Taylor-Fourier y la
cantidad de componentes del modelo), m
´
as costoso es el
algoritmo.
V. M
´
ETODOS BASADOS EN EL DISE
˜
NO DE FILTROS
Una limitaci
´
on evidente de los algoritmos basados en el
modelo de Taylor-Fourier es que tienen un comportamiento
pobre ante cambios abruptos en la amplitud y/o fase del
sincrofasor, pues la aproximaci
´
on polin
´
omica de una funci
´
on
con discontinuidades presenta un sobrepico fundamental por
el cl
´
asico fen
´
omeno de Gibbs [31]. Esto motiva buscar otro
enfoque que no dependa de un modelo tan restrictivo para
resolver el problema de estimaci
´
on de sincrofasores. En este
sentido, podemos decir que los m
´
etodos que se discuten en
esta secci
´
on, basados en criterios de dise
˜
no de filtros, no
est
´
an ligados a un modelo particular. Recientemente, se han
propuesto criterios de dise
˜
no de filtros para estimaci
´
on de
sincrofasores, frecuencia y ROCOF basados en restricciones
en el dominio de la frecuencia [32], [33], que dan por
resultado m
´
ascaras espectrales. Dichas m
´
ascaras espectrales
pueden utilizarse luego como base para el dise
˜
no de fil-
tros, por ejemplo, mediante el cl
´
asico algoritmo de Parks-
McClellan [15]. Sin embargo, la gran limitaci
´
on de estos
criterios es que no incorporan expl
´
ıcitamente los requisitos
que se tienen en el dominio temporal, lo cual significa que
la formulaci
´
on del problema en ambos casos es incompleta.
Una formulaci
´
on completa del problema de estimaci
´
on del
sincrofasor y sus derivadas, lo cual permite obtener la
frecuencia y la ROCOF con ecuaciones similares a las Ecs.
(46) y (47), con restricciones en el dominio del tiempo y
la frecuencia fue presentada en los trabajos [26], [34]. La
idea b
´
asica es que el problema de los dise
˜
nos de los filtros
se puede expresar como un programa de optimizaci
´
on semi-
infinita (SIP) [35]. Los SIPs son problemas de optimizaci
´
on
que involucran una cantidad finita de variables pero una
cantidad infinita de restricciones [36]. Concretamente, un
problema SIP en general se puede escribir como
min
x
f(x)
subject
to g
k
(x, y
k
) ≤ 0, y
k
∈ Y
k
, k = 1, . . . , p,
(57)
donde x ∈ R
n
son las variables de optimizaci
´
on, Y
k
⊂ R
m
k
son conjuntos infinitos (t
´
ıpicamente se asume que adem
´
as
son compactos), f : R
n
→ R es la funci
´
on costo, y g
k
:
R
n
×Y
k
→ R son las funciones de restricci
´
on del problema.
El conjunto de factibilidad del problema es
F = {x ∈ R
n
: g
k
(x, y
k
) ≤ 0 para
todo y
k
∈ Y
k
, k = 1, . . . , p}.
(58)
En particular, si f y g
k
(para cada y
k
∈ Y
k
) son funciones
convexas de x, decimos que tenemos un programa convexo
semi-infinito (CSIP) y el conjunto F es convexo. Esta es la
clase de problema que nos interesa y tiene la siguiente im-
portante propiedad: un m
´
ınimo local es tambi
´
en un m
´
ınimo
global [37]. Obviamente esto impacta fuertemente en las
propiedades de convergencia de los algoritmos que resuelven
problemas CSIP. Existen diversos solvers para resolver este
tipo de problemas de optimizaci
´
on [38], [39].
En nuestro problema, usamos un banco de filtro difer-
enciador que consiste de tres filtros de fase lineal H
i
(ν),
i = 0, 1, 2 [26], [34]. El primer filtro, H
0
(ν), que se utiliza
para obtener el estimador del sincrofasor, se elige como un
filtro FIR tipo I, para producir una respuesta en amplitud par
y un retardo de grupo entero. Del mismo modo, el segundo
filtro, H
1
(ν) , que es un filtro diferenciador de primer
orden, se elige como un filtro tipo III para producir una
respuesta en amplitud impar y un retardo de grupo entero.
Finalmente, el tercer filtro, H
2
(ν), que es un diferenciador
de segundo orden, nuevamente se elige como un filtro tipo
I. Sin p
´
erdida de generalidad, consideramos que los
´
ordenes
de los filtros son todos iguales a N − 1. Las expresiones
para las respuestas en amplitud son entonces
A
i
(ν) = g
T
i
(ν) a
i
, i = 0, 1, 2,
donde g
0
(ν) = [1, cos(ν), . . . , cos(Rν)]
T
, g
1
(ν) =
[sin(ν), . . . , sin(Rν)]
T
, g
2
(ν) = g
0
(ν), y nuevamente R =
(N − 1)/2. Los vectores a
i
representan los coeficientes del
filtro en la representaci
´
on de fase lineal [15].
Claramente, en nuestro caso, las variables de optimizaci
´
on
son los coeficientes a
i
de los filtros. Como se se
˜
nala en
los trabajos [26], [34], si se plantea el problema del dise
˜
no
de los filtros A
i
(ν) en forma secuencial, se obtienen tres
programas CSIP, uno para cada filtro. La elecci
´
on de las
funciones costo para los problemas que se propone en [26]
es
f
i
(a
i
) = λ
i
A
i
− A
id
i
2
2,Ω
1
+ (1 − λ
i
) kA
i
k
2
2,Ω
2
, (59)
donde i = 0, 1, 2, y Ω
1
, Ω
2
⊆ [−π, π] son conjuntos de
frecuencia dados. T
´
ıpicamente Ω
1
es la uni
´
on de la banda
de paso y la banda de transici
´
on de los filtros, y Ω
2
es la
banda de rechazo. Las funciones A
id
i
son las respuestas en
amplitud ideales de cada filtro en la banda de paso, es decir,
A
id
0
(ν) = 1, A
id
1
(ν) =
ν
T
, y A
id
2
(ν) = −
ν
2
T
2
. Adem
´
as, los
par
´
ametros λ
i
∈ [0, 1] se definen de acuerdo a la relaci
´
on de
compromiso deseada entre la distorsi
´
on de la se
˜
nal de inter
´
es
y la reducci
´
on de la energ
´
ıa de las se
˜
nales de interferencia y
el ruido. Notar que las funciones costo f
i
(·) son funciones
cuadr
´
aticas convexas en a
i
, de modo que se pueden escribir
como f
i
(a
i
) = a
T
i
P
i
a
i
+q
T
i
a
i
+r
i
, donde P
i
son matrices
semidefinidas positivas. Esta familia de funciones costo es
muy general. Por ejemplo, si tomamos λ
i
= 0 y Ω
2
=
[−π, π], obtenemos que la funci
´
on costo es proporcional a
la potencia de ruido a la salida de los filtros, de modo que
estar
´
ıamos minimizando la contribuci
´
on del ruido al error.
Las restricciones del problema est
´
an formuladas detal-
ladamente en el trabajo [26]. Aqu
´
ı, con el fin de ilustrar
el enfoque CSIP, presentaremos algunos ejemplos de la
formulaci
´
on de las restricciones semi-infinitas. Por ejemplo,
una gran limitaci
´
on del estimador cl
´
asico de la DFT consiste
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
87
http://elektron.fi.uba.ar
en su comportamiento pobre en condiciones de frecuencia
no nominal, tal como se discuti
´
o ampliamente en la Secci
´
on
III. Cuando se contempla que la se
˜
nal tiene una frecuencia
no nominal, tenemos el modelo de se
˜
nal y
n
= Xe
jδn
,
donde δ pertenece al conjunto de frecuencias de la banda de
paso Ω
pb
. Si forzamos al TVE a estar acotado por una cota
superior TVE
STA
, relacionada con pero no necesariamente
igual al l
´
ımite del Std., obtenemos la siguiente restricci
´
on
convexa semi-infinita:
TVE = |A
0
(δ) − 1| ≤ TVE
STA
, δ ∈ Ω
pb
. (60)
En forma similar, si evaluamos el FE e imponemos que
sea menor a cierto valor, digamos FE
STA
, obtenemos la
restricci
´
on semi-infinita
FE =
1
2π
A
1
(δ)
A
∗
0
(δ)
−
δ
T
≤ FE
STA
, δ ∈ Ω
pb
. (61)
Cabe notar que en esta expresi
´
on, la respuesta en amplitud
A
∗
0
(parametrizada por a
∗
0
) representa la soluci
´
on del prob-
lema de optimizaci
´
on correspondiente al filtro A
0
, que se
debe resolver previamente para poder plantear el problema
de optimizaci
´
on correspondiente al filtro A
1
. La restricci
´
on
anterior es del tipo convexa semi-infinita en los coeficientes
a
1
. Por
´
ultimo, tal como se mencion
´
o en la Secci
´
on IV, el
RFE es exactamente nulo para este modelo de se
˜
nal, por el
hecho de que los filtros H
i
son de fase lineal. Por lo tanto,
no es necesario imponer ninguna restricci
´
on sobre A
2
para
este test.
Otra restricci
´
on que es interesante presentar es la que se
obtiene cuando se consideran los tests de saltos en fase y
amplitud que aparecen en el Std. [4], pues esta es la gran
limitaci
´
on que tienen los algoritmos basados en TFT y es
una de las grandes motivaciones para utilizar el enfoque
CSIP. Comenzamos considerando un salto en amplitud y
luego uno en fase. El modelo de se
˜
nal para un salto en
amplitud es
y
n
= X[n] = X(1 + ∆au[n]), n ∈ N , (62)
donde u[n] es el escal
´
on unitario, y ∆a es el salto en am-
plitud. En este caso, la estimaci
´
on del fasor es simplemente
b
X[n] = X(A
0
(0) + ∆as
0
[n]), n ∈ N , (63)
donde s
0
[n] es la respuesta al escal
´
on del filtro A
0
. Supong-
amos, por simplicidad, que el salto es positivo: ∆a > 0. El
caso ∆a < 0 se puede tratar en forma similar. Luego, el
sobrepico en amplitud (AO) en este caso se puede escribir
como
AO = max
n=0,...,R
|
b
X[n]| − |X[n]|
|X[n]|
= max
n=0,...,R
(A
0
(0) + ∆as
0
[n]) − (1 + ∆a)
1 + ∆a
, (64)
donde A
0
(0) = g
T
0
(0)a
0
y s
0
[n] = u
T
0
[n]a
0
, con
(u
0
[n])
k
0
=
u[n+k
0
]+u[n−k
0
]
2
para k
0
= 0, . . . , R. Luego,
la restricci
´
on AO ≤ AO
max
se puede expresar f
´
acilmente
como R+1 restricciones lineales en a
0
, que son por supuesto
compatibles con el enfoque CSIP. Del mismo modo, para un
salto en fase tenemos el siguiente modelo de se
˜
nal
y
n
= X[n] = Xe
j∆φu[n]
, n ∈ N , (65)
donde ∆φ es el salto en fase. En este caso, el estimador
dentro del bloque de datos se puede escribir como
b
X[n] = Xv
H
0
[n]a
0
, n ∈ N , (66)
donde (v
H
0
[n])
k
0
=
e
j∆φu[n+k
0
]
+e
j∆φu[n−k
0
]
2
,
(v
H
1
[n])
k
1
=
e
j∆φu[n+k
1
]
−e
j∆φu[n−k
1
]
2
, (v
H
2
[n])
k
2
=
e
j∆φu[n+k
2
]
+e
j∆φu[n−k
2
]
2
, k
0
, k
2
= 0, . . . , R, k
1
= 1, . . . , R.
Consideremos nuevamente un salto positivo, ∆φ > 0,
para simplificar la presentaci
´
on. El sobrepico en fase (PO)
entonces se puede escribir como
PO = max
n=0,...,R
∠
b
X [n] − ∠X [n]
∆φ
= max
n=0,...,R
∠v
H
0
[n]a
0
− ∆φ
∆φ
.
(67)
Desafortunadamente, la restricci
´
on PO ≤ PO
max
pro-
duce restricciones no convexas. Sin embargo, notando que
∠v
H
0
[n]a
0
oscilar
´
a en un peque
˜
no entorno alrededor de
∆φ luego del salto, podemos relajar la restricci
´
on, reem-
plaz
´
andola por
Im{v
H
0
[n]}a
0
≤ tan[∆φ(1 + PO
max
)] Re{v
H
0
[n]}a
0
,
(68)
para n = 0, . . . , R, que son R + 1 restricciones lineales
ordinarias en a
0
. Estas formulaciones permiten entonces
controlar precisamente, mediante los par
´
ametros AO
max
y
PO
max
, los sobrepicos en amplitud y fase cuando ocurren
cambios abruptos en dichos par
´
ametros, lo cual es una
soluci
´
on al fen
´
omeno de Gibbs que se presenta en los es-
timadores basados en TF. Para demostrar esto, presentamos
en la Fig. 5 la respuesta a un salto en fase de un estimador
CSIP y un estimador TFT de dos ciclos. El estimador
CSIP fue dise
˜
nado para satisfacer los requisitos del est
´
andar,
mientras que para el estimador TFT se utiliz
´
o una ventana
rectangular (o ventana de Kaiser con β = 0) pues con
β > 0 el estimador presenta un ancho de banda excesivo. La
diferencia en el sobrepico de ambos estimadores es notable.
De hecho, el sobrepico del estimador CSIP es despreciable
(est
´
a en el orden de 10
−10
%), mientras que el estimador
de TFT presenta un sobrepico igual a 8, 02%, que excede
considerablemente al l
´
ımite del 5% que establece el Std. para
el PMU clase P. Este es un claro ejemplo de los beneficios
de incluir expl
´
ıcitamente las restricciones del problema en
la etapa de dise
˜
no.
t [s]
4.96 4.98 5 5.02 5.04
Fase estimada [rad]
0
0.1
0.2
CSIP
TF
Escalón de fase
Fig. 5: Respuesta a un salto en fase de π/18 rad. de los
estimadores CSIP y TFT.
Es interesante observar que los estimadores que se ob-
tienen con el banco de filtros diferenciador dise
˜
nado con
un enfoque CSIP son generalizaciones de los estimadores
DFT, WDFT, TF, y WTF. A continuaci
´
on mostramos
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
88
http://elektron.fi.uba.ar
expl
´
ıcitamente esta conexi
´
on. Comencemos notando que
todos los estimadores se pueden interpretar como la re-
spuesta de filtros FIR, de modo que en general se pueden
parametrizar por uno o m
´
as vectores h
i
. Los estimadores
DFT y WDFT se pueden parametrizar por un
´
unico vector
h
0
. Recordemos de la Secci
´
on III que el estimador DFT es
un estimador lineal insesgado de m
´
ınima varianza para el
modelo de fasor constante. Cuando se tiene ruido blanco,
la potencia de ruido a la salida es proporcional a kh
0
k
2
.
Luego, el estimador de la DFT se puede obtener como la
soluci
´
on del siguiente problema de optimizaci
´
on:
min
h
0
kh
0
k
2
sujeto a H
0
(0) = h
T
0
1 = 1, (69)
como puede verificarse f
´
acilmente [13]. Del mismo modo,
el estimador WDFT se puede interpretar como un estimador
lineal insesgado de m
´
ınima varianza cuando la matriz de
covarianza Σ del ruido v es modelada como una matriz
diagonal y donde las varianzas satisfacen
10
σ
2
n
= 1/w
n
. De
este modo, la interpretaci
´
on estad
´
ıstica de la WDFT es que
las muestras que son pesadas m
´
as fuertemente son muestras
de varianza peque
˜
na y vicecersa. Luego, la WDFT es la
soluci
´
on del siguiente problema de optimizaci
´
on
min
h
0
h
T
0
Σh
0
sujeto a H
0
(0) = h
T
0
1 = 1, (70)
donde Σ = W
−1
= diag({w
−1
n
}). Cabe notar que en este
planteo uno podr
´
ıa usar incluso informaci
´
on de correlaci
´
on
entre distintas muestras del ruido si se conociera o estimara
previamente, lo cual generaliza la WDFT. La soluci
´
on en
general ser
´
a
11
h
0
=
Σ
−1
1
1
T
Σ
−1
1
, (71)
de modo que el estimador es
b
X = h
T
0
y =
1
T
Σ
−1
y
1
T
Σ
−1
1
, (72)
que generaliza a la Ec. (33).
Consideremos ahora el estimador TF, que tambi
´
en es un
estimador lineal insesgado y de m
´
ınima varianza, pero para
un modelo de fasor polin
´
omico. En este caso, el estimador
queda determinado por L+1 filtros, de modo que la variable
de optimizaci
´
on es una matriz H = [h
0
, . . . , h
L
]. Para que
el estimador
b
X = H
T
y sea insesgado, debemos pedir que
E[
b
X] = X para cualquier X ∈ C
L+1
. Entonces, obtenemos
la condici
´
on
H
T
N = I. (73)
Por otro lado, notar que la potencia total de ruido a la
salida del banco de filtros diferenciador es proporcional
a kh
0
k
2
+ . . . + kh
L
k
2
= Tr(H
T
H). Luego, los filtros
10
Esta correspondencia tiene sentido pues generalmente se utilizan
ventanas positivas, es decir, tales que w
n
> 0 para todo n ∈ N.
11
La misma puede hallarse f
´
acilmente mediante el m
´
etodo de multipli-
cadores de Lagrange.
de Taylor-Fourier son soluciones al siguiente problema de
optimizaci
´
on:
min
H
Tr(H
T
H)
sujeto a H
T
N = I. (74)
An
´
alogamente a la generalizaci
´
on que se plante
´
o para inter-
pretar al estimador WDFT, el estimador WTFT se obtiene
resolviendo el siguiente problema de optimizaci
´
on, que
generaliza el problema anterior
min
H
Tr(H
T
ΣH)
sujeto a H
T
N = I, (75)
donde en este caso hacemos la asociaci
´
on Σ = (W
T
W )
−1
.
De hecho, la soluci
´
on a este problema es
12
H = Σ
−1
N (N
T
Σ
−1
N )
−1
, (76)
lo cual produce los estimadores
b
X = (N
T
Σ
−1
N )
−1
N
T
Σ
−1
y, (77)
que generaliza la Ec. (48). Por
´
ultimo, es interesante no-
tar que los problemas (74) y (75) se pueden desacoplar
f
´
acilmente en L + 1 problemas independientes, lo cual
permite optimizar la resoluci
´
on del problema en una forma
similar a la que se plantea en [30].
VI. CONCLUSIONES
En este trabajo hemos presentado un resumen de algunos
de los estimadores de sincrofasores m
´
as importantes que se
han desarrollado desde que el problema cobr
´
o inter
´
es en
los a
˜
nos 80 y, m
´
as importante a
´
un, la conexi
´
on entre ellos.
En primer lugar, analizamos los algoritmos m
´
as simples
basados en un modelo de fasor constante y, en forma
m
´
as general, en un modelo de fasor polin
´
omico. Luego,
presentamos los estimadores CSIP, un enfoque novedoso
basado en optimizaci
´
on semi-infinita convexa que ofrece
una gran flexibilidad de dise
˜
no. Concretamente, permite
incluir precisamente las restricciones que debe satisfacer el
sistema en el dominio del tiempo y de la frecuencia. En
particular, esto permite controlar tanto el comportamiento
del estimador en frecuencia no nominal, como los sobrepi-
cos, de modo que resuelve los problemas que afectan a los
estimadores basados en la DFT y la TFT. Hemos finalizado
el trabajo mostrando que existe una conexi
´
on, a trav
´
es de
la formulaci
´
on de problemas de estimaci
´
on de cuadrados
m
´
ınimos con distintas restricciones, entre los estimadores
DFT, WDFT, TF, WTF, y CSIP. Se espera que en el
desarrollo futuro de los estimadores de sincrofasores, se
propongan diversos esquemas adaptivos partiendo de filtros
prototipos dise
˜
nados con el enfoque CSIP, con el fin de
relajar las restricciones del problema original en pos de
un mejor desempe
˜
no en general. En el presente estamos
trabajando en esta tem
´
atica interesante.
12
Nuevamente el problema se puede resolver mediante multiplicadores
de Lagrange, planteando en este caso una matriz Λ de multiplicadores de
Lagrange. La funci
´
on de Lagrange en este caso ser
´
ıa, por ejemplo,
L(H, Λ) = Tr(H
T
ΣH) + Tr(Λ(I − H
T
N )).
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
89
http://elektron.fi.uba.ar
AGRADECIMIENTOS
El trabajo de F. Messina es financiado por una beca
de doctorado Peruilh de la Facultad de Ingenier
´
ıa de la
Universidad de Buenos Aires. Este art
´
ıculo fue posible
por la financiaci
´
on del proyecto FONARSEC UREE 4
”Sistema de Medici
´
on fasorial orientado al desarrollo de
redes inteligentes” del FONARSEC, Ministerio de Cien-
cia, Tecnolog
´
ıa e Innovaci
´
on Productiva; y el proyecto
PIP11220150100578CO del CONICET.
REFERENCIAS
[1] G. Giannakis, V. Kekatos, N. Gatsis, S.-J. Kim, H. Zhu, and B. Wol-
lenberg, “Monitoring and Optimization for Power Grids: A Signal
Processing Perspective,” IEEE Signal Process. Mag, vol. 30, no. 5,
pp. 107–128, Sept. 2013.
[2] F. Aminifar, M. Fotuhi-Firuzabad, A. Safdarian, A. Davoudi, and
M. Shahidehpour, “Synchrophasor measurement technology in power
systems: Panorama and state-of-the-art,” IEEE Access, vol. 2, pp.
1607–1628, 2014.
[3] J. D. L. Ree, V. Centeno, J. S. Thorp, and A. G. Phadke, “Syn-
chronized phasor measurement applications in power systems,” IEEE
Transactions on Smart Grid, vol. 1, no. 1, pp. 20–27, June 2010.
[4] IEEE Standard for Synchrophasor Measurements for Power Systems,
IEEE Std. C37.118.1-2011 (Revision of IEEE Std C37.118-2005),
Dec. 2011.
[5] IEEE Standard for Synchrophasor Measurements for Power Systems –
Amendment 1: Modification of Selected Performance Requirements,
IEEE Std. C37.118.1a-2014 (Amendment to IEEE Std C37.118.1-
2011), Apr. 2014.
[6] P. Castello, M. Lixia, C. Muscas, and P. Pegoraro, “Impact of the
Model on the Accuracy of Synchrophasor Measurement,” IEEE Trans.
Instrum. Meas., vol. 61, no. 8, pp. 2179–2188, Aug. 2012.
[7] G. Barchi, D. Macii, and D. Petri, “Synchrophasor estimators accu-
racy: A comparative analysis,” IEEE Trans. Instrum. Meas., vol. 62,
no. 5, pp. 963–973, May 2013.
[8] G. Barchi, D. Macii, D. Belega, and D. Petri, “Performance of
synchrophasor estimators in transient conditions: A comparative anal-
ysis,” IEEE Trans. Instrum. Meas., vol. 62, no. 9, pp. 2410–2418, Sept
2013.
[9] M. Karimi-Ghartemani and H. Karimi, “Processing of symmetrical
components in time-domain,” IEEE Trans. Power Syst., vol. 22, no. 2,
pp. 572–579, May 2007.
[10] J. Grainger and W. Stevenson, Power System Analysis. McGraw-Hill,
1994.
[11] F. Messina, P. Marchi, L. Rey Vega, C. G. Galarza, and H. Laiz,
“A Novel Modular Positive-Sequence Synchrophasor Estimation Al-
gorithm for PMUs,” IEEE Transactions on Instrumentation and
Measurement, vol. 66, no. 6, pp. 1164–1175, June 2017.
[12] A. Phadke and J. Thorp, Synchronized Phasor Measurements and
Their Applications, ser. Power Electronics and Power Systems.
Springer US, 2008.
[13] S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation
Theory. Upper Saddle River, NJ, USA: Prentice-Hall, Inc., 1993.
[14] G. Turin, “An introduction to matched filters,” IRE Transactions on
Information Theory, vol. 6, no. 3, pp. 311–329, June 1960.
[15] A. Oppenheim and R. Schafer, Discrete-Time Signal Processing, ser.
Prentice-Hall signal processing series. Pearson Education, 2011.
[16] D. Macii, D. Petri, and A. Zorat, “Accuracy of dft-based synchropha-
sor estimators at off-nominal frequencies,” in 2011 IEEE International
Workshop on Applied Measurements for Power Systems (AMPS), Sept
2011, pp. 19–24.
[17] G. Andria, M. Savino, and A. Trotta, “Windows and interpolation
algorithms to improve electrical measurement accuracy,” IEEE Trans-
actions on Instrumentation and Measurement, vol. 38, no. 4, pp. 856–
863, Aug 1989.
[18] C. Offelli and D. Petri, “The influence of windowing on the accuracy
of multifrequency signal parameter estimation,” IEEE Transactions
on Instrumentation and Measurement, vol. 41, no. 2, pp. 256–261,
Apr 1992.
[19] P. Stoica and R. Moses, Spectral Analysis of Sig-
nals. Pearson Prentice Hall, 2005. [Online]. Available:
https://books.google.com.ar/books?id=h78ZAQAAIAAJ
[20] D. Macii, D. Petri, and A. Zorat, “Accuracy Analysis and En-
hancement of DFT-Based Synchrophasor Estimators in Off-Nominal
Conditions,” IEEE Trans. Instrum. Meas., vol. 61, no. 10, pp. 2653–
2664, Oct. 2012.
[21] D. Belega and D. Petri, “Accuracy Analysis of the Multicycle Syn-
chrophasor Estimator Provided by the Interpolated DFT Algorithm,”
IEEE Trans. Instrum. Meas., vol. 62, no. 5, pp. 942–953, May 2013.
[22] D. Belega, D. Dallet, and D. Petri, “Statistical description of the sine-
wave frequency estimator provided by the interpolated dft method,”
Measurement, vol. 45, no. 1, pp. 109 – 117, 2012. [Online]. Available:
http://www.sciencedirect.com/science/article/pii/S0263224111003058
[23] J. A. de la O Serna, “Dynamic Phasor Estimates for Power System
Oscillations,” IEEE Trans. Instrum. Meas., vol. 56, no. 5, pp. 1648–
1657, Oct. 2007.
[24] M. A. Platas-Garza and J. A. de la O Serna, “Dynamic Phasor and
Frequency Estimates Through Maximally Flat Differentiators,” IEEE
Transactions on Instrumentation and Measurement, vol. 59, no. 7, pp.
1803–1811, July 2010.
[25] D. Petri, D. Fontanelli, and D. Macii, “A frequency-domain algorithm
for dynamic synchrophasor and frequency estimation,” IEEE Trans.
Instrum. Meas., vol. 63, no. 10, pp. 2330–2340, Oct. 2014.
[26] F. Messina, L. Rey Vega, P. Marchi, and C. G. Galarza, “Optimal
Differentiator Filter Banks for PMUs and their Feasibility Limits,”
IEEE Transactions on Instrumentation and Measurement, 2017, en
prensa.
[27] D. Belega, D. Fontanelli, and D. Petri, “Dynamic phasor and fre-
quency measurements by an improved taylor weighted least squares
algorithm,” IEEE Transactions on Instrumentation and Measurement,
vol. 64, no. 8, pp. 2165–2178, Aug 2015.
[28] J. A. de la O Serna, “Synchrophasor measurement with polynomial
phase-locked-loop taylor-fourier filters,” IEEE Transactions on Instru-
mentation and Measurement, vol. 64, no. 2, pp. 328–337, Feb 2015.
[29] M. A. Platas-Garza and J. A. de la O Serna, “Dynamic harmonic
analysis through taylor-fourier transform,” IEEE Transactions on
Instrumentation and Measurement, vol. 60, no. 3, pp. 804–813, March
2011.
[30] ——, “Polynomial implementation of the taylor-fourier transform
for harmonic analysis,” IEEE Transactions on Instrumentation and
Measurement, vol. 63, no. 12, pp. 2846–2854, Dec 2014.
[31] A. Jerri, The Gibbs Phenomenon in Fourier Analysis,
Splines and Wavelet Approximations, ser. Mathematics and
Its Applications. Springer US, 1998. [Online]. Available:
https://books.google.com.ar/books?id=y2c00LAMnFUC
[32] D. Macii, G. Barchi, and D. Petri, “Design criteria of digital filters for
synchrophasor estimation,” in 2013 IEEE International Instrumenta-
tion and Measurement Technology Conference (I2MTC), May 2013,
pp. 1579–1584.
[33] A. J. Roscoe, B. Dickerson, and K. E. Martin, “Filter design masks
for c37.118.1a-compliant frequency-tracking and fixed-filter m-class
phasor measurement units,” IEEE Trans. Instrum. Meas., vol. 64,
no. 8, pp. 2096–2107, Aug 2015.
[34] F. Messina, P. Marchi, L. Rey Vega, and C. G. Galarza, “Design
of Synchrophasor Estimation Systems with Convex Semi-Infinite
Programming,” in Electrical Power and Energy Conference (EPEC),
Ottawa, Canada, Oct. 2016.
[35] A. W. Potchinkov, “Design of optimal linear phase fir filters by a semi-
infinite programming technique,” Signal Processing, vol. 58, no. 2,
pp. 165–180, 1997.
[36] R. Reemtsen and J. R
¨
uckmann, Semi-Infinite Programming, ser.
Nonconvex Optimization and Its Applications. Springer US, 1998.
[37] S. Boyd and L. Vandenberghe, Convex Optimization, ser. Berichte
¨
uber verteilte messysteme. Cambridge University Press, 2004.
[38] MATLAB, Optimization Toolbox. Natick, Massachusetts: The Math-
Works Inc., 2016.
[39] “Semi-infinite programming solvers: CSIP and NSIPS,”
http://www.norg.uminho.pt/aivaz/software.html.
Revista elektron, Vol. 1, No. 2, pp. 79-90 (2017)
ISSN 2525-0159
90
http://elektron.fi.uba.ar