4. METODOS NUMERICOS UTILIZADOS EN LA RESOLUCION DEL SISTEMA by xpr28091

VIEWS: 66 PAGES: 43

									4.   METODOS NUMERICOS UTILIZADOS EN LA RESOLUCION DEL SISTEMA DE
     ECUACIONES NO LINEAL.


     4.1.    Introducción.


     El estudio del flujo de cargas lleva a plantear un sistema de ecuaciones no lineales que
se expresará como un sistema no lineal igualado a cero de la forma,

                                                                                         (4.1)

     Este sistema puede presentar múltiples soluciones matemáticamente posibles y su
resolución numérica debe proporcionar la solución físicamente correcta, x(s). Existen distintos
métodos para la resolución del sistema de ecuaciones (4.1).


     El método de Newton es el más utilizado en la bibliografía, [7], [19] y [22], para la
resolución del flujo armónico de cargas. Este método tiene como ventaja una relativa
simplicidad y una rápida convergencia pero tiene el inconveniente que su convergencia
únicamente se puede asegurar si los valores iniciales de las incógnitas están lo suficientemente
cerca de la solución (dominio de atracción de la solución).


     Así en el estudio del flujo de cargas convencional, [20] y [21], y del flujo armónico de
cargas, [22], para sistemas mal condicionados el método de Newton puede presentar
problemas de convergencia. Estos problemas son más frecuentes en el flujo armónico de
cargas dada la dificultad en la inicialización de la fase de las tensiones armónicas.


     Con el fin de intentar mejorar la convergencia del método de Newton existen en la
bibliografía los métodos de Newton modificados [20], [26], [27], [28] y [29], basados algunos
de ellos en criterios de optimización. Estos métodos aún suponiendo una mejora respecto al
método de Newton no aseguran la convergencia global a la solución correcta del problema
porque pueden converger hacía mínimos locales del sistema tratado.


     Por otra parte, frente a la resolución del flujo armónico de cargas mediante el método
de Newton, limitado por la elección de las condiciones iniciales, existe como alternativa para


                                           - 40 -
los casos de no convergencia los métodos de continuación, [22], [30], [31], y en concreto el
método de Davidenko en los que no existen problemas debido a la elección de las condiciones
iniciales aunque, por contra, el proceso de cálculo es más lento. En este método la elección
del valor inicial no es crítica, siempre que corresponda a una situación físicamente real, ya
que la homotopía conduce de dicho valor hacia la solución.


     4.2.    Método de Newton.


     El método de Newton partiendo de unas condiciones iniciales x(0) genera una secuencia
de valores x(1), .., x(i), .., x(m) que convergen a la solución x(s) si las condiciones iniciales,
x(0), están suficientemente próximas a ella. En particular, tal como se puede ver en la
referencia [22], consiste en la aplicación iterativa del algoritmo,

                                                                                           (4.2)

desde el valor inicial, x(0). Donde DF(x) es el jacobiano de la función F(x). El proceso
iterativo finaliza cuando F(x) 2<∈ (un valor usual de ∈ es 10-4).




                      Figura 4.1. Convergencia del método de Newton.


     La convergencia hacia una solución del algoritmo está asegurada por el teorema de
Newton-Kantorovich, [6], si se cumplen ciertos requisitos. De estos, los más importantes son
que el punto inicial debe estar suficientemente cerca de la solución y dentro de un


                                            - 41 -
determinado recinto, donde F(x) sea continuamente diferenciable. Esto queda ilustrado en la
figura 4.1. donde el proceso numérico diverge o tiende a otra solución si la inicialización es
x2(0) por estar fuera del dominio de atracción de la solución, x(s), mientras que el proceso
converge si la inicialización es x1(0).


     Así, este algoritmo tiene la ventaja de una rápida convergencia si las condiciones
iniciales, x(0), están próximas a la solución buscada. Es por ello que una mala elección de las
tensiones armónicas iniciales puede dar problemas de convergencia, los cuales son,


-    El método es divergente, [7].
-    El método converge hacia una solución falsa, [6] y [23].
-    El método converge después de muchas iteraciones, [19].


tal como se estudiará en ejemplos posteriores.


     Los problemas de convergencia son más habituales en el estudio del flujo armónico de
cargas porque los valores iniciales de los ángulos de las tensiones armónicas pueden ser muy
diferentes de los ángulos solución. Así, los sistemas mal condicionados son más frecuentes
en el flujo armónico de cargas que en el flujo de cargas convencional.


     4.3.    Métodos de Newton modificados


     Los métodos de Newton modificados incorporan al algoritmo de Newton, (4.2), un
factor amortiguante, λ, reescribiéndose el algoritmo de la forma,

                                                                                        (4.3)

     El objetivo de este factor es amortiguar el efecto del término,




para evitar la divergencia del método de Newton debido a "saltos" excesivamente grandes en
alguna iteración.



                                          - 42 -
     Existen distintos criterios en la bibliografía para su determinación,


     En la referencia [29] y [32] el factor amortiguante tiene la expresión




así, para cada iteración i se buscará el valor del parámetro ni (partiendo de ni=0) que lleve
a cumplir la condición F(x(i+1)) ∞< F(x(i)) ∞.


     En la referencia [26] se presentan 3 variantes de los métodos de Newton modificados
atendiendo al criterio del cálculo de λ(i).


     La primera consiste en utilizar dos factores amortiguantes fijos. Un primer factor, S1,
inferior a la unidad, p.ej. λ(i) = S1 = 0.5, para las primeras iteraciones hasta que la función
error fuese inferior a un determinado valor ∈, F(x(i)) ∞< ∈. Y un segundo factor, S2, mayor
que el primero, p. ej. λ(i) = S1 = 1.0, para finalizar el proceso iterativo. Con ello se controla
la convergencia del método en las primeras iteraciones, que suelen ser las más críticas, para
luego acelerar el proceso aumentando el valor de λ(i).


     La segunda consiste en que en cada iteración el factor de la variante anterior, λ(i) = S1
ó S2, es reducido siempre que el término ∆x(i) = DF(x(i))-1 F(x(i)) supere un cierto límite L1,
p. ej. ∆x(i) > L1 = .25 p.u. Y su determinación se realiza imponiendo que las tensiones
armónicas no superen un límite L2, p. ej. x(i+1) < L2 = .4 p.u. Es decir,




     La tercera variante es similar a la presentada en la referencia [29] y [32] pero utiliza un
factor amortiguante distinto para cada armónico, por lo que el algoritmo (4.3) se reescribe de
la forma,




donde [λ(i)] es un vector cuyos términos son,



                                              - 43 -
     Idénticamente que en [29] y [32] para cada iteración y, en esta variante, para cada
armónico se buscará el valor de nik que lleve a cumplir la condición general
 F(x(i+1)) ∞< F(x(i))   ∞   y las condiciones particulares   Fk(x(i+1)) ∞< Fk(x(i)) ∞.


     En la referencia [20], formulando el problema en coordenadas rectangulares, el
multiplicador óptimo para una iteración i, λ(i) se determina minimizando la función H( λ),




donde, considerando el sistema (4.1) como [F(x)] = [G(x)] - [Gdat] = 0, se tiene,




     Los métodos de Newton modificados mejoran la convergencia del método de Newton
al controlar los "saltos" del proceso iterativo con el factor amortiguador pero no aseguran
dicha convergencia pudiéndose alcanza mínimos locales de la función estudiada.


     Así, en estos casos el factor multiplicador λ(i) obtenido en una iteración es muy próximo
a cero, [29], por lo que el proceso iterativo queda detenido y no converge a la solución del
problema. Esto ocurre al alcanzar un punto de jacobiano singular, DF(x(i)) =0, ya que la
reducción del término DF(x(i))-1 F(x(i)) impone un factor λ(i) cercano a cero.


     4.4.     Métodos de continuación. Método de Davidenko.


     En los métodos de continuación a partir del sistema no lineal (4.1), por medio de una
parametrización, se fabrica una homotopía o camino,




que lleva de un valor arbitrario x(0) a la solución del problema, x(s). La homotopía es
seleccionada tal que,


                                               - 44 -
por lo que los valores iniciales y finales son solución del problema original.


     Diferenciando la función implícita respecto a t se obtiene la ecuación diferencial,




cuya resolución numérica con las condiciones iniciales proporcionará la curva x(t) que finaliza
para tm en la solución.


     El procedimiento que se detalla aquí es el método de Davidenko, [31], que pertenece
a la amplia gama de los métodos de continuación [25].


     En el método de Davidenko es definida la homotopía exponencial,

                                                                                        (4.4)

donde el parámetro t va de 0 a ∞ y x(0) es un valor inicial arbitrario.


     Así, la homotopía va desde un valor inicial arbitrario,




a la solución del problema,




de forma continua.


     Derivando la homotopía (4.4) respecto a t es posible obtener el problema de ecuaciones
diferenciales,

                                                                                        (4.5)

     La resolución numérica del problema, [22], nos permitirá obtener el camino que recorre


                                          - 45 -
la homotopía al variar t. Este camino finaliza para t ⇒ ∞ en la solución, x(s). En la práctica,
es suficiente alcanzar t=5 para obtener una buena aproximación a la solución.


     La resolución numérica de la ecuación de Davidenko, (4.5), puede ser obtenida con
cualquier algoritmo de ecuaciones diferenciales [25] (p.ej. Runge-Kutta, predictor-corrector,
..). Usando el método de Euler se tiene el algoritmo iterativo,

                                                                                        (4.6)

donde el paso de integración h es constante.


     Como se puede comprobar el algoritmo obtenido a través de la resolución de la ecuación
diferencial por el método de Euler coincide con el propuesto en los métodos de Newton
modificados (4.3) y si se utiliza un paso de integración unidad, h=1, se obtiene el algoritmo
de cálculo del método de Newton (4.2).




                    Figura 4.2. Convergencia del método de Davidenko.

     La aplicación del algoritmo (4.6) da lugar, aproximadamente, a un conjunto de
aplicaciones sucesivas del método de Newton, en el cual podemos estar todo lo cerca que se
quiera de la solución local de cada paso del proceso de resolución. Así, debido a la especial

                                          - 46 -
forma de evolucionar el método para asegurar la convergencia a una solución x(s) basta con
que se cumplan las soluciones del teorema local de Newton, a diferencia del teorema de
Kantorovich, que es un teorema de convergencia global. Esto queda ilustrado en la figura 4.2.


     La condición crítica que impone el teorema local de Newton es que el jacobiano de F(x)
sea invertible.


     Así, la utilización de un paso de integración grande puede dar problemas de estabilidad
en el proceso si las condiciones iniciales no están próximas a la solución. Es por ello que se
deberá utilizar un paso de integración suficientemente pequeño con el fin de asegurar la
convergencia pero ello supondrá que el proceso de cálculo sea muy lento.


     La ventaja del método de Davidenko, es que, partiendo de una situación físicamente real
se asegura la convergencia de forma continua a la solución del problema sin que la elección
de las condiciones iniciales deba ser realizada por el usuario. Así tomando como condiciones
iniciales las de reposo del sistema (inicialización trivial) el proceso seguido por el método
puede ser interpretada como una lenta y continua puesta en marcha del sistema hasta las
condiciones de carga establecidas. Esta interpretación física ofrece interesantes posibilidades
al estudio del flujo armónico de cargas como comprobaremos en los ejemplos posteriores.


     Como desventaja está la lentitud del proceso hasta que se obtienen los resultados finales.


     4.5.     Elección de las condiciones iniciales.


     En la implementación del método de Newton es el usuario quien decide los valores de
las condiciones iniciales y, como ya se ha comentado, una mala elección de dichos valores
puede llevar a la divergencia del método.


     En este sentido cuando se utiliza para la resolución del flujo de cargas convencional, sin
armónicos, no suele presentar problemas porque las tensiones de los nudos de la red son
próximas a la del nudo slack,
y éstas suelen ser condiciones iniciales suficientemente buenas aunque tal como se menciona


                                            - 47 -
en la referencia [25] existen casos de redes mal condicionadas donde se pueden presentar
problemas de convergencia.


     En el estudio del flujo de cargas con armónicos existirá el problema de la elección de
las condiciones iniciales para las tensiones armónicas, Vi(k) para i=1, .., n y k=5, 7, .., ya que
no serán comparables a la fundamental del slack, la única conocida, y no se tendrá ningún
criterio de elección en especial para la fase de dichas tensiones (el modulo es predecible que
sea de valor reducido si el índice de distorsión armónica es bajo), lo que puede llevar a
importantes problemas de convergencia tal como se menciona en [19].


     La inicialización más habitual utilizada en la bibliografía, [19], para las tensiones es,




existiendo distintos criterios para la inicialización de los parámetros del convertidor según sea
su modelización. De acuerdo con [6],




donde α2(0) se deduce de (2.2) sin considerar la presencia de armónicos.


     Los métodos de Newton modificados siguen el mismo criterio que el método de Newton
en la elección de las condiciones iniciales y la utilización del factor multiplicador que se
incorpora al algoritmo perseguirá reducir los problemas de convergencia que aparecen debido
a la dificultad en la inicialización de las tensiones armónicas.


     El planteamiento utilizando los métodos de continuación, y en concreto el método de
Davidenko, permite un enfoque distinto del problema de las condiciones iniciales dado que
su elección, siempre que corresponda a una situación real, no será crítica por establecerse un

                                            - 48 -
camino continuo a través de la homotopía entre dichos valores y la solución del problema.


     El hecho de poder elegir el punto de inicio invita a hacerlo como el punto de "reposo"
del sistema, el cual suele tener una solución trivial. En el flujo armónico de cargas tal
situación corresponderá a la de un sistema que no consume potencia, por tanto no existe
inyección de armónicos por parte de los convertidores. Esta circunstancia permite representar
la evolución del método de continuación como una "lenta" puesta en marcha del sistema lo
que facilita el análisis de los resultados intermedios y apoya la idea de que la solución hallada
sea la físicamente correcta. El método de Davidenko permite dar la interpretación física de
que la homotopía se inicia en el estado del sistema relajado y progresivamente éste evoluciona
hasta alcanzar para t ⇒ ∞ el régimen de funcionamiento planteado en el problema,
obteniéndose para dicho régimen las soluciones buscadas. La idea de comenzar el proceso en
el estado de reposo de la red lleva a realizar las siguientes inicializaciones.


     El consumo de potencia será nulo en las cargas y los convertidores. Así,
aproximadamente, la siguiente inicialización es posible,




     No existe inyección inicial de armónicos, al no existir consumo de potencia por parte
de los convertidores, por lo que se deberían inicializar las tensiones armónicas de todos los
nudos a cero. Esto provoca problemas de dependencia lineal en el jacobiano en la primera
iteración por lo que se inicializarán las tensiones como las producidas por una inyección
armónica correspondiente al consumo de los convertidores de una potencia con un valor
próximo a cero (pequeña fracción de su potencia nominal, Pd = 10-3 Pdato) con Ud = UdDato.


     Así, los valores iniciales de las tensiones armónicas son calculados con,




donde Ibus(k) es la corriente inyectada cuando los convertidores AC/DC consumen un pequeño
porcentaje de su potencia nominal. Con esta pequeña modificación el sistema continua estando
próximo a la condición de reposo.



                                           - 49 -
     Para los convertidores se ha considerado como condición de reposo un consumo de
potencia nulo mientras que la tensión en el lado de continua no es cero debido a que la
tensión en el lado de alterna es la del slack por lo que estas condiciones deberán ser
trasladadas a las variables que definen su comportamiento. Así, la inicialización de los
parámetros del convertidor según [6],




     Al ser nulo el consumo de potencia las admitancias fundamentales de las cargas tendrán
un valor también nulo y las correspondientes admitancias armónicas adoptarán el valor que
se deduzca de su función imponiendo la condición anterior.




     4.6.    Comentario a los métodos presentados.


     Se han presentado el método de Newton, las modificaciones existentes sobre dicho
método y el método de Davidenko para la resolución de los sistemas de ecuaciones no
lineales que se obtienen en el flujo armónico de cargas.


     Como ya se ha comentado la incorporación del factor amortiguante al método de
Newton, método de Newton modificado, mejora su convergencia aunque no la asegura
pudiendo aparecer en el proceso iterativo problemas de convergencia hacia mínimos locales
de F(x) .


     El método de continuación, por una parte ofrece un enfoque nuevo de la resolución
numérica a través de las ecuaciones diferenciales lo que permitirá utilizar todas sus
herramientas para el tratamiento del problema. Así, resolviendo la ecuación diferencial (4.5)
con el método de Euler se obtiene de forma justificada al algoritmo propuesto por los
métodos de Newton modificados dando una justificación teórica a la utilización del factor
multiplicador λ y su influencia en la convergencia del proceso al corresponder al paso de
integración h.

                                          - 50 -
     Por otro lado ofrece un enfoque nuevo a la inicialización del problema proponiendo unas
condiciones iniciales, situación físicamente real de la red a estudio, que junto a la integración
del sistema diferencial derivado de la homotopía (4.4) con un paso suficientemente pequeño
asegurará la convergencia a la solución del problema de forma continua, aunque presenta el
inconveniente de la lentitud del proceso de cálculo.


     Así comparando el método de Newton con el método de continuación se podrían
considerar ambos procedimientos complementarios pareciendo deseable combinar sus ventajas
de manera que, frente a un problema, primero se pruebe la convergencia del método de
Newton y si éste falla se aproxime a la solución con el método de Davidenko de paso
variable, estableciéndose un compromiso entre la rapidez y la seguridad a la hora de
seleccionar dicho paso.


     Todo esto se puede comprobar en el ejemplo de la figura 4.3,




                                Figura 4.3. Sistema a estudio.


el cual ha sido estudiado en el punto 6.1 con los datos de la tabla 4.1 considerando 6
armónicos,
     Así para este ejemplo, según se ve en la figura 4.4, el método de Newton diverge en 5
iteraciones, el método de Davidenko converge de forma continua a la solución correcta pero
su tiempo de ejecución es excesivamente largo y el método de Newton modificado presentado
en las referencias [29] y [32] alcanza un mínimo local de la función F(x)         ∞   (el proceso

                                           - 51 -
       PARAMETRO        VALOR (p.u.)       PARAMETRO         VALOR (p.u.)
               V1           1.05                 Yr               .00
               X            .341                 Yc               .0455
               R            .02                  U2D              1.3
                              Tabla 4.1 Datos del sistema.

numérico se detiene al obtenerse un factor amortiguante λ próximo a cero).




                         Figura 4.4. Evolución de los procesos.

     En la figura 4.4 se ha representado la norma infinito de la función error, F(x) ∞, en
trazo discontinuo, y el ángulo de disparo del convertidor, α2, en trazo continuo frente al
parámetro t.




                                        - 52 -
5.   NUEVO PLANTEAMIENTO PARA EL ESTUDIO DEL FLUJO ARMONICO DE
     CARGAS.


     5.1.    Introducción.


     Se ha visto que en los trabajos existentes en la bibliografía con relación al flujo
armónico de cargas se pueden diferenciar tres puntos básicamente,


- La modelización de los elementos del sistema.
- La formulación del flujo de cargas.
- Y la resolución numérica del sistema de ecuaciones no lineales.


     En el presente capítulo se pretende abordar el análisis del flujo armónico de cargas con
relación a los últimos dos puntos.


     Respecto a la formulación del flujo armónico de cargas se ha visto que una de las
simplificaciones más comunes que se realizan en la bibliografía es la de considerar que la
potencia consumida es debida únicamente a la onda fundamental con lo que se simplifica el
problema al reducirse el número de incógnitas a tratar (desaparecen como variables no
conocidas las admitancias fundamental y armónica).


     Esta simplificación sólo es correcta, como se vera en un ejemplo posterior, si la
aportación armónica de potencia frente a la fundamental es reducida. Es por ello que se
propone realizar la formulación completa del problema, método 1 (apartado 3.2), considerando
la potencia armónico-fundamental e incorporando a dicha formulación los nudos PV.


     Respecto a la resolución numérica del sistema de ecuaciones no lineales planteado en
la formulación del flujo de cargas con armónicos se han presentado tres métodos,


- Método de Newton.
- Método de Newton modificado.
- Métodos de continuación. Método de Davidenko.


                                         - 53 -
     El primero presenta las ventajas de una relativa simplicidad y de una rápida
convergencia si las condiciones iniciales están próximas a la solución, pero tiene el
inconveniente de la dificultad en la elección de dichas condiciones iniciales.


     El segundo es una modificación del método de Newton que reduce sus problemas de
convergencia pero no asegura dicha convergencia.


     Y el tercero tiene la ventaja de asegurar la convergencia del proceso a la solución
correcta siempre que las condiciones iniciales correspondan a una situación físicamente real,
pero tiene el inconveniente de la lentitud del proceso de cálculo.


     Así se presentará un nuevo método de resolución numérica para el problema del flujo
de cargas, llamado método de h-Newton, que busca combinar las ventajas de la rapidez en
la convergencia del método de Newton con la seguridad en dicha convergencia del método
de Davidenko.


     Así, el nuevo planteamiento para el estudio del flujo armónico de cargas pretende
abordar el problema a través de una formulación completa del mismo, a la que se han
incorporado los nudos PV, y una resolución numérica del sistema no lineal convergente y
rápida.


     5.2.    Formulación completa del flujo armónico de cargas.


     Se desarrolla la formulación planteada en el método 1 incorporando los nudos PV al
problema tal como se muestra en la tabla 5.1.


     En la nueva formulación se distinguirán dos tipos de nudos PV,


- Nudos PV de generador (generador PV). Son datos la tensión y la potencia entregada en
  el nudo por venir impuestas por el generador.


- Nudos PV de carga (carga PV). Son datos la tensión, por venir fijada mediante una batería


                                          - 54 -
  de condensadores de valor desconocido, y el consumo de potencia activa en el nudo.

        1                   :      Slack o referencia (S)
        2   ag              :      Generador PQ (GQ)
        g+1 a h             :      Generador PV (GV)
        h+1 a c             :      Carga PV (CV)
        c+1 a l             :      Carga PQ (CQ)
        l+1 a n             :      Dispositivo no lineal (D)

                                Tabla 5.1 Tipos de nudos de la red.


     Así el planteamiento del problema considerando que a representa el número de
armónicos más el fundamental vendrá definido por,


    NUDO          DATOS               INCOGNITAS                NUM. DE INCOGNITAS

      (S)          V1(1)             V1(k) (k=5, 7, ..)                   2 (a - 1)

     (GQ)           Ss               Vs(k) (k=1, 5, ..)                 2 a (g - 1)

     (GV)         Ph y Vh            Vh(k) (k=1, 5, ..)                 2 a (h - g)

     (CV)         Pc y Vc        Vc(k) y Yc(k) (k=1, 5, ..)              4 a (c - h)

     (CQ)           St           Vt(k) y Yt(k) (k=1, 5, ..)              4 a (l - c)

     (D)         Pd y UdD       Vd(k) (k=1, 5, ..) y αd, βd           2 (n - l) (a + 1)

                                                   TOTAL:      2 (a (n + l - h) + n - l - 1)


  Tabla 5.2 Planteamiento completo del flujo armónico de cargas con incorporación de
                                             nudos PV.

siendo Pd la potencia activa en el lado de alterna y UdD la tensión en el lado de continua del
convertidor.


     Para obtener las 2 (a (n+l−h)+n−l−1) incógnitas del problema se deben plantear el
mismo número de ecuaciones.



                                              - 55 -
      Así se tienen las ecuaciones presentadas en el método 1 (apartado 3.2),


      Balance de potencias: Las 2 (g-1)+ 2 (n-c) ecuaciones de la potencia activa y reactiva,
expresiones (3.3), para los nudos de generador PQ, excepto el slack, de carga PQ y de
convertidor.


      Balance de corrientes: Para el nudo slack y los nudos de generador PQ y carga PQ, las
2 g (a-1) + 2 (l-c) (a-1) ecuaciones correspondientes a la expresión (3.4).


      Corrientes inyectadas: Las 2 (n-l) (a-1) ecuaciones correspondientes a la expresión (3.5)
en los nudos no lineales.


      Ecuaciones de los convertidores: Para todos los convertidores de la red (i=l+1, .., n) se
plantean las 2 (n-l) expresiones correspondientes a (3.6).


      Para los nudos de carga PQ (i=c+1, .., l) se tendrán las 2 (l-c) restricciones
correspondientes a sus potencias consumidas y las 2 (l-c) (a-1) ecuaciones que definen el
valor de las admitancias para cada orden de armónico k=5, 7, .., expresiones (3.7) y (3.8)
respectivamente.


      Y se deben añadir las necesarias para cubrir las 2 a (h-g)+4 a (c-h) incógnitas
añadidas por los nudos PV.


      Estas ecuaciones son las siguientes,


      Balance de potencia activa: (c-g) expresiones de la potencia activa en todos los nudos
PV de la red (i=g+1, .., c). Utilizando las expresiones (3.1) y (3.2) resultan la ecuaciones no
lineales,


                                                                                        (5.1)


donde k representa el orden del armónico perteneciendo dicho valor, siempre que las cargas


                                             - 56 -
sean simétricas, al conjunto K = (k∈N, k=6n+1 ó k=6n-1 con n=0, 1, ..) secuencias directa
e inversa respectivamente.


     Valor eficaz de la tensión: (c-g) expresiones del valor eficaz de la tensión en todos los
nudos PV (i=g+1, .., c),

                                                                                       (5.2)


     Balance de corrientes: Teniendo en cuenta que en los nudos PV no se inyectan
armónicos en la red o sea se comportan como una carga con admitancia Yh = gh + jbh
conocida correspondiente a la admitancia del generador (nudos PV generador) en secuencia
directa e inversa ó Yc = gc + jbc correspondiente a la admitancia de la carga y de la batería
de condensadores que regula la tensión (nudos PV carga). Para un orden de armónico dado
(k = 5, 7, ..) al aplicar el método de los nudos (3.1) resulta la expresión,


                                                                                       (5.3)


y al separar las partes real e imaginaria resultan 2 (c-g) (a-1) expresiones,




     Además, para los nudos PV de carga (i=h+1, .., c) se tendrán las 2 (c-h) restricciones
correspondientes a sus potencias consumidas,


                                                                                       (5.4)



donde la potencia reactiva consumida por la carga PV no es un dato, ya que la batería de

                                           - 57 -
condensadores utilizada para regular la tensión en dichos nudos es de valor desconocido, por
lo que será calculada desde el lado de la red con la expresión,




obtenida de las expresiones (3.1) y (3.2).


     Y las 2 (c-h) (a-1) expresiones que definen el valor de las admitancias para cada orden
de armónico k = 5, 7, ..,

                                                                                      (5.5)

     Así el número de ecuaciones es,


                   ECUACIONES                 NUMERO DE ECUACIONES

                       (3.3) ..............       2 (g - 1) + 2 (n - c)

                   (5.1) y (5.2) ........                2 (c - g)

                   (3.4) y (5.3) ........               2 l (a - 1)

                       (3.5) ..............          2 (n - l) (a - 1)

                       (3.6) ..............              2 (n - l)

                       (5.4) ..............              2 (c - h)

                       (3.7) ..............              2 (l - c)

                       (5.5) ..............          2 (c - h) (a - 1)

                       (3.8) ..............          2 (l - c) (a - 1)

                              TOTAL:           2 (a (n + l - h) + n - l - 1)



                               Tabla 5.3 Balance de ecuaciones.

que corresponden, según lo presentado, a las siguientes expresiones,




                                              - 58 -
- 59 -
     5.3.    Método de h-Newton.


     Frente a los problemas de convergencia que presenta el método de Newton en los
sistemas mal condicionados, y que pueden llevar al método a divergir o a converger a una
solución falsa, se ha visto que existe en la bibliografía los métodos de Newton modificados
que incorporan un factor amortiguante λ para intentar asegurar dicha convergencia. Sin
embargo estos métodos no pueden asegurar la convergencia global del proceso porque pueden
converger hacia mínimos locales tal como veremos en ejemplos posteriores.


     Por otra parte se ha presentado el método de continuación, método de Davidenko, que
trabajando con un algoritmo de ecuaciones diferenciales asegura la convergencia del proceso
a la solución correcta siempre que las condiciones iniciales correspondan a una situación
físicamente real y el paso de integración h sea lo suficientemente pequeño pero tiene el
inconveniente de que su tiempo de ejecución es excesivamente largo.


     El método propuesto, cuyo nombre es método de h-Newton, es un método numérico de
convergencia global, no divergente y comparable en tiempo de ejecución al método de
Newton cuando éste converge.


     Así, el nuevo método estará basado en dos ideas,


- La sustitución del algoritmo de Newton por un algoritmo para la resolución numérica de
  ecuaciones diferenciales.


- La elección de la red descargada (consumo de potencia nulo en las cargas y los
  convertidores) como condiciones iniciales.


     La primera idea permite usar los algoritmos numéricos típicos de ecuaciones
diferenciales para la resolución del problema. Así, el método de Euler (4.6) de paso variable,
en cada iteración i se puede utilizar un paso h(i) diferente, parece ser el apropiado, ya que el
único objetivo de la integración es alcanzar la solución correcta, (t⇒∞), tan rápido como sea
posible.


                                           - 60 -
     Por ello, el algoritmo (4.6) se reformula de la forma,

                                                                                        (5.6)

     La velocidad de integración (h(i) grande) irá en contra de la estabilidad del proceso
numérico. Por ello será necesario un detector de estabilidad como criterio de elección del paso
óptimo, h(i) (compromiso entre la estabilidad y la rapidez). El detector de pérdida de
estabilidad utilizado es la norma infinito de F(x), F(x) ∞. Cuando, F(x(i+1)) ∞> F(x(i)) ∞,
h(i) debe ser reducida hasta obtener, F(x(i+1)) ∞< F(x(i) ∞.




                    Figura 5.1. Convergencia del método de h-Newton.

     La elección de la red descargada como condición inicial (segunda idea presentada) y el
control de la estabilidad en la integración, evitará obtener soluciones matemáticamente
posibles pero físicamente no reales, [6].


     La necesidad de comenzar el proceso numérico con la red descargada es claro en el
método de Davidenko (el cual utiliza el método de Euler con un paso de integración fijo y
pequeño, p. ej. h=.025). En este caso, la evolución del estado de la red puede ser interpretada
como una "lenta y continua" puesta en marcha del sistema, desde el consumo de potencia
nulo hasta el consumo nominal. Así, esta suave evolución siempre lleva a la solución correcta

                                            - 61 -
(evitando otras soluciones matemáticamente posibles) si la inversa del jacobiano existe. Pero
el hecho de trabajar con un paso reducido tiene el inconveniente de hacer que el proceso
numérico sea excesivamente largo. Si este proceso alcanzase un punto de jacobiano singular,
esto significaría que la red está en una situación físicamente imposible (colapso de tensiones).
Las condiciones iniciales podrían ser cualquier situación físicamente real de la red, pero la
situación de descarga de la red es una solución trivial del sistema.


     Así, el proceso numérico sería similar al planteado para el método de Davidenko, figura
5.1, (sucesiva resolución de problemas intermedios hasta alcanzar la solución) pero
seleccionando la zona de convergencia mayor a través del criterio de la norma infinito,
F(x(i+1)) ∞< F(x(i) ∞.


     El algoritmo del nuevo método es,


     1)   Leer la información y los datos del sistema.
     2)   Formar las matrices de admitancia, YBUS(k) con k=1,5, ..
     3)   Inicializar los parámetros, n(i)=n0=0 y i=0.
     4)   Inicialización, x(0) y cálculo de F(x(0)) y DF(x(0)).
     5)   Calcular, h(i)n=1/2n.
     6)   Calcular, x(i+1) = x(i) - h(i)n DF(x(i))-1 F(x(i)).
     7)   Calcular, F(x(i+1)).
     8)   Si F(x(i+1)) 2<∈ (Convergencia obtenida) ir a 11).
     9)   Si F(x(i+1)) ∞< F(x(i))   ∞   (criterio de estabilidad correcto) entonces
          9.1) i=i+1 y n(i)=n0.
          9.2) Actualizar, DF(x(i)).
          sino
          9.3) n(i)=n(i)+1.
          Fin Si
     10) Ir a 5)
     11) Imprimir resultados.
     12) Fin.



                                              - 62 -
trabajando con un valor de ∈ igual a 10-4.


      Así, el control de la estabilidad del algoritmo h-Newton permite seguir aproximadamente
la misma evolución que el método de Davidenko pero de una forma más rápida. Por tanto,
el método propuesto consiste en integrar la ecuación diferencial (4.5) con el método de Euler
(4.6) usando un paso de integración variable, h(i), desde las condiciones del sistema
descargado. La elección del paso h(i) se realiza buscando el compromiso entre la rapidez y
la estabilidad.


      Hasta este punto el método propuesto sería similar a los métodos de Newton
modificados pero la utilización de los algoritmos de ecuaciones diferenciales para su
caracterización permite justificar el uso del factor del factor amortiguador h(i)=λ(i) así como
su determinación, y permite incorporar la idea de inicializar el proceso numérico con la red
descargada. Por otra parte también se justifica la reducción del paso para evitar la divergencia
del método pues ésta sería debida a un problema de la resolución numérica (paso de
integración excesivamente grande).


      El nuevo método, al igual que los métodos de Newton modificados, podría presentar
problemas de convergencia hacia mínimos locales por lo que para evitar dichos problemas se
incluye la siguiente sentencia,


      5’) Si h(i)n<hmin (Problemas de convergencia detectados) entonces
           5.1’) n(i)=n0=n0+1 y i=0.
           5.2’) Ir a 4)
           Fin Si


después de 5). Los problemas de convergencia son detectados controlando que el paso no sea
inferior a un determinado valor (p.ej., hmin = 0.025). Así la progresiva reducción del valor
inicial del paso, hn0, convertirá el método de h-Newton en el método de Davidenko. Como
ya ha sido comentado anteriormente el método de Davidenko es un proceso continuo basado
en la resolución de una ecuación diferencial y esto asegura la convergencia del proceso a la
solución correcta evitando los problemas de mínimos locales que aparecerían en los


                                           - 63 -
algoritmos discretos, figura 5.2.




           Figura 5.2. Problemas de mínimos locales en la resolución numérica.

     Así el método de h-Newton es muy simple porque únicamente incorpora dos
modificaciones respecto a los métodos de Newton amortiguados.


-    La reducción del multiplicador amortiguante inicial, hn0, cuando el método alcanza un
     mínimo local.


-    La elección como valores iniciales del proceso numérico iterativo, x(0), los
     correspondientes a la red descargada (consumo de potencia nulo).


ambas condiciones aseguran la convergencia del método de h-Newton a la solución correcta.


     Por último, el método de h-Newton obliga a elegir una situación de consumo de
potencia nulo como punto inicial del proceso. Este punto de inicio (punto de reposo del
sistema de potencia) es una solución trivial del sistema e impone la inicialización presentada
en el apartado 4.5.



                                          - 64 -
     5.4.   Estructura del programa de resolución del flujo de cargas con armónicos.


     El programa ha sido realizado en un ordenador PC-486 y el lenguaje utilizado ha sido
FORTRAN 5.0 compilado para ser ejecutado bajo entorno WINDOWS.


     El flujograma general del programa se muestra en la figura 5.3,




                      Figura 5.3. Flujograma general del programa.

en el destacan los siguientes módulos.


     5.4.1. Módulo 1. Lectura de datos.


     En este módulo se realizará la lectura de todos los datos necesarios para poder
desarrollar el estudio del flujo armónico de cargas correspondiente. Destacan dos tipos de
datos,


     > Los datos del proceso; que serán introducidos por el usuario a petición del programa.
         Estos son,



                                          - 65 -
         -   Planteamiento del flujo de cargas.


             -   Consideración del consumo de potencia armónico y fundamental.
             -   Consideración del consumo de potencia sólo fundamental.


         -   Método numérico de resolución.


             -   Método de Newton.
             -   Método de Davidenko.
             -   Método de h-Newton.


         -   Archivo de datos del sistema.
         -   Número de armónicos a tratar.


     Así, el programa permite resolver el problema considerando las potencias armónico-
fundamentales (método 1) o sólo fundamentales (método 2) con la incorporación de los nudos
PV. Y también permite abordar la resolución de las ecuaciones no lineales por los tres
métodos numéricos expuestos anteriormente. De esta forma se podrá realizar un análisis
comparativo de todos estos aspectos.


     > Los datos del sistema; serán leídos por el programa del archivo de datos que se ha
         especificado. Estos son,


         -   Datos de los nudos.
         -   Datos de las líneas y transformadores.


     5.4.2. Módulo 2. Inicialización de variables.


     El proceso de inicialización dependerá del método de resolución numérica seleccionado.


     En la implementación del método de Newton es el usuario quien decide los valores de
las condiciones iniciales. Esto se realizará a través del archivo de datos del sistema de forma


                                          - 66 -
que el programa leerá de dicho archivo la inicialización elegida por el usuario.


       El planteamiento utilizando el método de Davidenko y por el método h-Newton permite
un enfoque distinto del problema de las condiciones iniciales dado que su elección, siempre
que corresponda a una situación real, no será crítica. Así, el hecho de poder elegir el punto
de inicio invita a hacerlo como el punto de "reposo" del sistema, el cual suele tener una
solución trivial. También sería posible partir de un punto correspondiente a un estado del
sistema conocido. Es por ello que la inicialización para estos métodos ya estará establecida
y se realizará a través de una subrutina del programa de acuerdo a lo expuesto en el punto
4.5.


       5.4.3. Módulo 3. Determinación de la matriz Ybus(k).


       La matriz de admitancias de la red para un armónico k dado será,




siendo Yim(k) = Gim(k) + jBim(k) el elemento de la matriz de admitancias entre los nudos i
y m para un armónico k dado.


       Así, el elemento Yim(k) está definido, en redes sin acoplamientos, por,




donde n es el número de nudos de la red en contacto con el nudo j, sin contar el de
referencia, es decir, tierra.


       Debido al posible tamaño de las matrices Ybus(k) y al alto contenido en elementos nulos
que pueden tener se ha trabajado con matrices dispersas con el fin de evitar la ocupación de
memoria de forma innecesaria.




                                           - 67 -
     Así la matriz dispersa sólo trabaja con los elementos no nulos y los distribuye, para cada
armónico k, de la siguiente forma,


GIA(kmax, nv) y BIA(kmax, nv)        Elementos no nulos de las matrices Gbus(k) y Bbus(k)
                                     escritos por filas.
GJA(kmax, nv) y BJA(kmax, nv)        Columna de los elementos no nulos de las matrices
                                     Gbus(k) y Bbus(k) escritos por filas.
GA(kmax, ny) y BA(kmax, ny)          Posición de las matrices YIA y YJA donde se inicia la
                                     fila "ny".


donde nv es el número de elementos no nulos de la matriz Ybus(k), ny es el número de filas
más 1 de la matriz Ybus(k) y kmax es el número de armónicos a tratar más la onda
fundamental.


     5.4.4. Módulo 4. Resolución del flujo de cargas.


     El programa permitirá al usuario elegir la resolución numérica del flujo de cargas,


     - Método de Newton.
     - Método de Davidenko.
     - Método de h-Newton.


     El programa también permitirá al usuario seleccionar la formulación del flujo de cargas,


     -         Consideración de las potencias armónicos y fundamentales.
     -         Consideración de las potencias sólo fundamentales.


     El proceso para la resolución del flujo de cargas será el mismo para ambas
formulaciones y lo único que variará será el cálculo de la función F(x) y de su jacobiano
DF(x) o sea la formulación matricial del problema, así, en función del método numérico
seleccionado se tendrá el siguiente algoritmo según lo presentado en el punto 5.3,



                                          - 68 -
1)   Leer la información y los datos del sistema.
2)   Formar las matrices de admitancia, YBUS(k) con k=1,5, ..
3)   Inicializar los parámetros, n(i)=n0=0 y i=0.
4)   Inicialización, x(0) y cálculo de F(x(0)) y DF(x(0)).
5)   Si
     5.1)     Método Newton entonces h(i)n=1.0.
     5.2)     Método Davidenko entonces h(i)n=0.025          (a)
                                                                   .
     5.3)     Método h-Newton entonces h(i)n=1/2n.
     Fin Si
6)   Si h(i)n<hmin y Método h-Newton entonces
     6.1)     n(i)=n0=n0+1 y i=0.
     6.2)     Ir a 4)
     Fin Si
7)   Calcular, x(i+1) = x(i) - h(i)n DF(x(i))-1 F(x(i)).
8)   Calcular, F(x(i+1)).
9)   Si F(x(i+1)) 2<∈ ir a 11).
10) Si Método h-Newton
     10.1) Si F(x(i+1)) ∞< F(x(i))     ∞   (criterio de estabilidad correcto) entonces
              10.1.1) i=i+1 y n(i)=n0.
              10.1.2) Actualizar, DF(x(i)).
              sino
              10.1.3) n(i)=n(i)+1.
              Fin Si
     10.2) Ir a 5)
     sino
     10.3) i=i+1.
     10.4) Actualizar, DF(x(i)).
     10.5) Ir a 7)
     Fin Si
11) Imprimir resultados.
12) Fin.
              (a) El valor de h(i)n puede ser seleccionado por el usuario.


                                           - 69 -
         5.4.5. Módulo 5. Salida de resultados.


         El programa presentará los resultados en dos pantallas distintas atendiendo al tipo de
estos.


         La primera indicará si el proceso ha convergido o no y en caso de convergencia ofrece
al usuario el tiempo empleado en la ejecución del flujo de cargas, el número de iteraciones
y el error, F(x)) 2, en el que se ha finalizado el proceso.


         La segunda pantalla corresponde al valor obtenido de las incógnitas del problema tal
como se verá en los ejemplos posteriores.


         5.5.   Formulación matricial del problema.


         La resolución del problema del flujo armónico de cargas por el método de h-Newton
consiste en la aplicación iterativa del algoritmo (5.6) desde las condiciones iniciales, x(0). Así
matricialmente dicho algoritmo quedará planteado de la forma,




que se reformula de la forma,

                                                                                           (5.7)


donde,


[∆x(i)] = [x(i+1)] - [x(i)]     Vector de incógnitas.
[F(x(i))]                       Vector de funciones de error.
[DF(x(i))]                      Matriz jacobiana de las funciones de error.


         Cada uno de los elementos anteriores y su correspondiente cálculo dependerá del tipo
de formulación del flujo armónico de cargas que se esté utilizando. En el programa
desarrollado existe la posibilidad de trabajar con dos formulaciones distintas.


                                            - 70 -
     5.5.1 Formulación         del flujo de cargas considerando las potencias armónico-
              fundamentales.


     El vector de incógnitas, según lo planteado en el punto 5.2, tendrá la siguiente
estructura,




donde,




     Como ya se ha visto el estudio del flujo de cargas lleva a plantear un sistema de
ecuaciones no lineales que se expresará como un sistema no lineal igualado a cero, vector de

                                          - 71 -
funciones de error, de la forma,




donde,




                                   - 72 -
     Para ello el sistema de ecuaciones del flujo de cargas se reformulan de la forma,




     Y por último se deberá determinar el jacobiano de las funciones error, [DF(x)]. Así la
expresión (5.7) quedará,

                                        - 73 -
     Debido al posible tamaño de la matriz jacobiana [DF(x)] y al alto contenido en
elementos nulos que pueden tener se ha trabajado con matrices dispersas con el fin de evitar
la ocupación de memoria de forma innecesaria.


     Así la matriz dispersa sólo trabaja con los elementos no nulos de la siguiente forma,


JIA(ne) Elementos no nulos de las matriz dispersa del jacobiano, [DF(x)], escritos por filas.
JJA(ne) Columna de los elementos no nulos de la matriz dispersa del jacobiano, [DF(x)],
         escritos por filas.
JA(nf)   Posición de las matrices JIA y JJA donde se inicia la fila "ne".


donde ne es el número de elementos no nulos de la matriz [DF(x)], nf es el número de filas
más 1 de la matriz [DF(x)].


     El cálculo de la matriz del jacobiano se realizará numéricamente utilizando un paquete
informático con el algoritmo 618, [38], que estima el jacobiano de una función F(x)
trabajando con matrices dispersas.


     En problemas de gran escala el jacobiano suele estar constituido por una matriz dispersa
lo que hace interesante la estimación que se propone.


     Para estimar el jacobiano se plantea,




y tomando un valor de d suficientemente pequeño se tiene,

                                          - 74 -
las subrutinas comprendidas en el algoritmo 618, [38], se basan en el cálculo anterior y en
la dispersidad del jacobiano.


     Definida la dispersidad del jacobiano se pueden agrupar sus columnas en grupos que
cumplen la propiedad de no tener ningún elemento no nulo en la misma fila, por ejemplo para
una matriz,




se realizarían dos grupos, el primero con las columnas 1 y 4, y el segundo con 2 y 3, es decir,




donde se ve que todas las filas sólo tienen un elemento no nulo.


     Una vez realizada esta partición se procede a calcular el jacobiano de cada grupo e
incorporarlo a la matriz jacobiana total. Así por ejemplo para el primer grupo,




     El proceso de estimación del jacobiano utilizando el algoritmo 618, [38], será,



                                          - 75 -
      - Definición de la dispersidad del jacobiano. Se deben de contabilizar los elementos
          no nulos que constituirán el jacobiano y para cada elemento se indicará la fila y la
          columna donde se encuentra.
      - Llamada a la subrutina DSM( ) para realizar la partición en los grupos cuyas
          columnas no tienen ningún elemento no nulo en la misma fila.


Para cada grupo,


      - Cálculo de F(x) y de F(x + d).
      - Cálculo del jacobiano con la subrutina FDJS( ).


      Así, únicamente se deberá definir la dispersidad del jacobiano [DF(x)] localizando los
elementos no nulos. Esta dispersidad queda reflejada en la fig. 5.4 donde las incógnitas son,


θs1 = θs(1) y Vs1 = Vs(1)              (s = 2, .., h)
θt1 = θt(1) y Vt1 = Vt(1)              (t = h+1, .., l)
θd1 = θd(1) y Vd1 = Vd(1)              (d = l+1, .., n)
θsk = θs(k) y Vsk = Vs(k)              (s = 1, .., h ; k = 5, 7, ..)
θtk = θt(k) y Vtk = Vt(k)              (t = h+1, .., l ; k = 5, 7, ..)
θdk = θd(k) y Vdk = Vd(k)              (d = l+1, .., n ; k = 5, 7, ..)
αd y βd                                (d = l+1, .., n)
gt1 = gt(1) y bt1 = bt(1)              (t = h+1, .., l)
gtk = gt(k) y btk = bt(k)              (t = h+1, .., l ; k = 5, 7, ..)


y las funciones error son,


FPi             (i=2, .., c)                 FItk = FIt(k)               (t=h+1,.., l ; k=5, ..)
FPd             (d=c+1, .. n)          FIdk = FId(k)            (d=l+1, .., n ; k=5, ..)
FQi             (i=2, .., g)                 FUPDd ⇒ FUDd, FPDd (d=l+1, .., n)
FVi             (i=g+1, .., c)         FSt                      (t=h+1, .., l)
FQd             (d=c+1, .., n)         FYtk = FYt(k)                     (t=h+1, .., l ; k=5, ..)
FIsk = FIs(k) (s=1, .., h ; k=5, ..)


                                             - 76 -
Figura 5.4. Dispersidad del jacobiano considerando las potencias armónico-fundamentales.

                                       - 77 -
y las distintas zonas definen la dispersidad de la matriz,


i)       Zona con todos los elementos nulos.
ii)      Zona correspondiente a una matriz diagonal.
iii)     Zona con todos los elementos no nulos.
iv)      Zona con la dispersidad caracterizada por la figura 5.5.




                        Figura 5.5. Caracterización de la dispersidad iv)

       5.5.2. Formulación del flujo de cargas considerando las potencias fundamentales.


       El vector de incógnitas,




                                           - 78 -
donde,




     El sistema de ecuaciones no lineales expresado como un sistema no lineal igualado a
cero corresponde al vector de funciones de error de la forma,




donde,




                                         - 79 -
       El sistema de ecuaciones no lineales del capítulo 3 se reformula de la forma,




       Y por último se deberá determinar el jacobiano de las funciones error, [DF(x)]. Así la
expresión (5.7) quedará,




y la caracterización de la dispersidad del jacobiano vendrá definida por la figura 5.6,


i)       Zona con todos los elementos nulos.
ii)      Zona correspondiente a una matriz diagonal.
iii)     Zona con todos los elementos no nulos.
iv)      Zona con la dispersidad caracterizada por la figura 5.5.


                                           - 80 -
Figura 5.6. Dispersidad del jacobiano considerando la potencia fundamental.

                                 - 81 -
     Finalmente comentar que para la resolución del sistema lineal disperso caracterizado por
la expresión (5.7) se ha realizado por métodos directos utilizando las subrutinas
correspondientes al algoritmo 533, [39].




                                           - 82 -

								
To top