Continuación Homotópica

Partiendo de lo anterior, la continuación homotópica es un método numérico de convergencia global, es decir, que es capaz de encontrar múltiples soluciones a sistemas de ecuaciones algebraicas no lineales mediante una aproximación inicial. El método de continuación homotópica tiene como premisa el resolver sistemas de ecuaciones algebraicas no lineales a través de la deformación continua de una función G(x) hasta llegar a convertirse en otra función continua f(x) perteneciente al sistema de ecuaciones no lineales. Este mapeo se lleva a cabo con la ayuda del parámetro homotópico τ para todas las ecuaciones del sistema (^{Diaz-Arango et al., 2018}). De esta manera el sistema se vuelve un sistema homotópico expresado como:

donde Hfx=Hx,τ:Rn+1→Rn, x ∈Rn τ∈0,1.

Dada la ecuación 1, y tomando en cuenta las condiciones impuestas para Hx,τ, es conveniente definir a la ecuación Gx de manera que tenga una solución propuesta o conocida; de este modo, se cumplen las condiciones propias de la continuación homotópica:

Si τ=0 y H(x, τ)= Gx, entonces la solución propuesta es obtenida de manera muy sencilla (dada su definición).

En contraste, si τ=1 y H(x, τ)=f(x), entonces la solución del sistema original es obtenida (^{Oliveros-Muñoz & Jiménez-Islas, 2013}).

Con respecto a estos tres primeros conceptos (ecuación 1, condición 1 y condición 2), resulta útil visualizar que el proceso de variación, llamado en lo sucesivo SHH, consiste en deformar ese sistema intencionalmente construido para que sea fácil resolver Gx hasta que se convierta en f(x), todo a través de la variación paulatina y sistemática de τ, desde 0 hasta 1. Un ejemplo geométrico sencillo para integrar este proceso de manera intuitiva es considerar la clásica deformación homotópica necesaria para obtener un toroide a partir de una taza (Figura 1).

Fuente: Elaboración propia.

Figura 1 Interpretación esquemática de la deformación de una taza para obtener un toroide al llevar τ de 0 a 1. 

Los sistemas algebraicos asociados al toroide y la taza son bastante complejos; el ejemplo es usado únicamente para establecer el concepto intuitivo de deformación geométrica asociada a la ecuación homotópica (1).

Una de las consecuencias de la condición descrita para el dominio de las variables de la ecuación 1 (Rn+1→Rn, x ∈Rn τ∈0,1) es la noción de hiperespacio, que denota la generalización de las propiedades del espacio euclídeo (donde operan las geometrías clásica y analítica) de tres a cuatro o más dimensiones. Téngase en consideración que todo sistema con tres o más variables (y ecuaciones) cumple con x ∈Rn (en este caso n sería 3), y al considerar el parámetro homotópico τ∈0,1, se tiene entonces Rn+1→Rn (que en este ejemplo es R4→R3). Esta aclaración es importante porque el SHH hace uso de propiedades euclídeas de las figuras regulares, como la definición R3 de la esfera en sistemas con tres o más ecuaciones.

En sistemas algebraicos, las formas más conocidas y sencillas para Gx (^{Diaz-Arango et al., 2018}; ^{Oliveros-Muñoz & Jiménez-Islas, 2013}) son denominadas:

Aunque existen más formas homotópicas Gx disponibles en la literatura científica, dada la intención didáctica del presente trabajo, solo se analizará la homotopía de Newton, la cual se sustituye en la ecuación 1 para así obtener:

Deformación homotópica y Parametrización del sistema

Para generar los vectores de predicción se necesitan las pendientes de dichos vectores, estas se obtienen de las ecuaciones que componen el camino homotópico. Por otra parte, se mencionó anteriormente que el sistema de ecuaciones algebraicas no lineales se deforma por medio de la adición de una ecuación G(x) y su respectivo parámetro homotópico τ. Dicha deformación da como resultado la ecuación 1, de la cual se puede deducir:

Si se considera que H(x, τ) es también función de un parámetro p, con el cual se estima la distancia recorrida a lo largo del espacio homotópico, entonces todas las variables implicadas (x∈Rn y τ∈ 0,1) tienen funcionalidad de p.

De esta manera se parametriza el sistema.

Definición de matriz Jacobiana

Una vez que se tiene el sistema homotópico, se procede a agrupar el vector de x y el parámetro homotópico τ en un vector auxiliar y.

Diferenciando la ecuación anterior, y empleando la regla de la cadena, se obtiene:

donde ∂H∂yi es una columna del Jacobiano H' del sistema homotópico, por lo que se tiene:

De esta manera se obtiene un sistema de n ecuaciones con n+1 variables. Al tratarse de un sistema homogéneo donde el número de variables es mayor al número de ecuaciones, existen soluciones diferentes a la trivial.

Cálculo de vector predictor

Para realizar el cálculo de los vectores de predicción se toman en cuenta las pendientes del camino homotópico mencionadas en el punto anterior y el algoritmo de Euler, que en este caso corresponde a:

donde xip y τp son los valores de predicción obtenidos por el método de Euler que se encontrarán en la superficie de la hiperesfera de radio r. Por otro lado, para calcular el avance horizontal ∆p se hace uso de las ecuaciones 14 y 15 en conjunto con la ecuación de la hiperesfera:

Por lo tanto, ∆p se calcula con la siguiente fórmula:

donde el signo dependerá de si el seguimiento se realiza a la derecha o a la izquierda, es decir, en sentido positivo o negativo. De esta manera es posible calcular los vectores de predicción para x y τ de acuerdo con las ecuaciones 14 y 15.

Para calcular el punto p correspondiente al avance del camino homotópico se utilizará la ecuación:

Corrección hiperesférica

Una vez calculados los vectores de predicción, estos se vuelven el valor inicial para empezar con la corrección hiperesférica, dicho procedimiento tiene como objetivo el encontrar la intersección entre la hiperesfera y el camino homotópico.

Este procedimiento se realiza empleando el método de Newton-Raphson para resolver las ecuaciones homotópicas del sistema, descritas como la ecuación 2 en conjunto con la ecuación de la hiperesfera, ecuación 16.

Recordemos que el método de Newton-Raphson requiere una matriz de derivadas parciales, denominado Jacobiano (J), y un vector con la función evaluada en un punto inicial (B) para comenzar la corrección. Esto puede verse de forma más clara con la expresión siguiente:

donde los valores del vector X se sumarán a los valores del vector predictor para obtener un vector corrector:

Este proceso también es iterativo, se recomienda un rango de tres a cinco iteraciones para obtener un buen resultado.

Por lo tanto, para un sistema de dos ecuaciones, el Jacobiano (J) quedaría expresado de la siguiente forma:

donde c es la ecuación de la hiperesfera derivada parcialmente en cada uno de los componentes del sistema. De manera similar, el vector denominado B, que contiene las funciones homotópicas evaluadas en los valores del vector de predicción y la función de la hiperesfera, quedará expresado como:

De esta manera el punto corregido se convertirá en el centro de la siguiente hiperesfera y, a su vez, en el punto inicial del siguiente vector de predicción.

Detección y refinamiento de cruces con τ=1

De acuerdo con los puntos anteriores, para que el camino homotópico pueda trazarse, es necesario un alto número de iteraciones. Estas iteraciones tienen que monitorearse para detectar el momento en el que el camino homotópico cruce el hiperplano τ=1, ya que, de acuerdo con la teoría, este será el punto en el que se encuentre el vector solución del sistema.

Sin embargo, es casi imposible que el cruce del camino se localice en el punto exacto de τ=1. Por lo tanto, cuando se detecte el paso de un número inferior (τ<1) a un número superior a (τ>1), o viceversa, es cuando hay que realizar el refinamiento de cruces con τ=1.

Este proceso se llevará a cabo con el mismo método que la corrección hiperesférica, es decir, se empleará el método de Newton-Raphson. En esta ocasión se contará solamente con las funciones del sistema de ecuaciones algebraicas no lineales originales y el último vector corregido que haya sobrepasado a τ=1, esto aplica tanto en el sentido +τ como en la dirección contraria -τ.

Para un sistema de dos ecuaciones se expresaría como:

donde el subíndice R indica que se trata de elementos usados en el refinamiento de cruces con τ=1. Finalmente, se recomienda un número mínimo de cinco iteraciones para obtener el vector solución lo más exacto posible.

El pseudo-código del algoritmo de SHH puede encontrarse en el Apéndice B de este mismo artículo.

Ejemplo práctico

Definición de matriz Jacobiana

Tomando en cuenta lo anterior, se procede a obtener el Jacobiano (H') del sistema derivando parcialmente las ecuaciones homotópicas. De esta manera se obtiene:

Sustituyendo el Jacobiano H' del sistema se obtiene:

Sustituyendo los valores iniciales:

Cálculo de pendientes aplicando la Regla de Cramer

En este paso se tiene por objetivo obtener los determinantes de cada uno de los elementos del sistema, (xi,τ). Este proceso se puede resumir a la ecuación 13, la cual expresa el cálculo de los determinantes de matrices a las cuales se les ha retirado la columna i perteneciente a la variable en turno. Esto es:

Aplicando la ecuación 13 para el cambio de signo, los determinantes quedan expresadas como:

Cálculo de vector predictor

El vector predictor se calculará usando las pendientes que hemos calculado en el punto anterior, es decir, y˙i:

Las ecuaciones necesarias para calcular el vector de predicción son las ecuaciones 14 y 15, las cuales requieren un elemento ∆p que se calcula con la ecuación 17. Sustituyendo los valores de ∂xi∂p, ∂τ∂p y r en la formula, la ecuación se reduce a:

Posteriormente, se calcula el vector de predicción con las ecuaciones 14 y 15. Esto da como resultado:

Por último, se calcula el avance de p.

Corrección hiperesférica

Para llevar a cabo la corrección hiperesférica se calculará un Jacobiano (J) de acuerdo con la ecuación 21, detallada previamente. En dicho Jacobiano se sustituirán los valores pertenecientes al vector predictor y al vector inicial de la siguiente manera:

donde xi_0 permanecerá constante, mientras que xi0 cambiará en cada iteración de la corrección con el punto anterior. Por otro lado, el vector B, que contiene las funciones homotópicas evaluadas en el punto predictor y la ecuación de la hiperesfera, se definirá de acuerdo con la forma de la ecuación 22.

Una vez definida la matriz del Jacobiano (J) y el vector B, se sustituyen los valores de manera que xi0, xi_0, τ0, xi_p, xip y τp contengan los valores del vector inicial y el vector de predicción.

El siguiente paso es obtener la matriz inversa de J, es decir J-1:

Dicha matriz se multiplica por el vector B para obtener un vector X que se sumará al vector de predicción para obtener un nuevo vector X.

En la siguiente iteración el vector X ocupará el lugar del vector de predicción hasta que se cumplan las iteraciones del proceso de corrección hiperesférica, se recomienda que sea de tres a cinco iteraciones. Finalmente, el último vector X será el vector corregido Xc que pasará a ser el punto inicial xi0 para la siguiente iteración global del proceso de continuación homotópica (no confundir con xi_0).

Detección y refinamiento de cruces con τ=1

El proceso iterativo de continuación homotópica anteriormente explicado fue realizado durante cuatro ciclos cuando se registró un cruce en τ=1 donde τ cambia de 0.9636 a 1.0540. Dado que no se trata del punto exacto τ=1, se procede a refinar la solución utilizando el último punto corregido xi_c como punto inicial para el método de Newton-Raphson descrito por la ecuación 23, donde la matriz Jacobiana correspondiente al proceso de refinado JR es igual a:

y el vector BR corresponde a:

Posteriormente, se calculó la inversa de JR y se multiplicó por el vector BR para obtener el vector XR de acuerdo con la ecuación 23 que se mencionó anteriormente.

Dicho vector XR se sumó al vector inicial para obtener el nuevo vector XR.

El nuevo vector XR pasará a ser el vector inicial para la segunda iteración en el proceso de refinamiento de cruces en τ=1. Este proceso se repitió hasta que los valores del vector BR, que contiene las funciones del sistema original evaluadas con el vector inicial en turno, se volvieron prácticamente 0 (Es decir, los valores están por debajo de la tolerancia fijada de antemano, por ejemplo 10^-5.) De esta manera es como se obtuvo el primer vector solución del sistema, para este caso particular es de:

En este caso de estudio se realizó un seguimiento de camino homotópico corto con 12 hiperesferas, las cuales fueron suficientes para encontrar los dos vectores solución conocidos para este sistema. Dichos vectores se muestran en la Tabla 1.

Todos los vectores corregidos obtenidos durante el seguimiento hiperesférico se muestran en la Tabla 2, donde los cruces en τ=1 están resaltados en color gris oscuro. Se registra un total de tres cruces que fueron refinados con el método de Newton-Raphson, sin embargo, los dos primeros cruces resultaron ser el mismo vector solución, por lo que se omitió en la Tabla 1. En el Apéndice A, la Figura A1 muestra las gráficas 3D de las superficies H1(x1,x2,τ) y H2(x1,x2,τ) en un momento inicial (τ=0), y la Figura A2 muestra las seis primeras deformaciones y su influencia sobre las intersecciones de las funciones homotópicas con H=0. En este caso, la deformación homotópica se representa por el cruce simultaneo de las dos superficies con el plano H=0, y la deformación homotópica implica la expansión de ambas superficies. En sistemas con más ecuaciones y variables, este tipo de visualización geométrica se dificulta; sin embargo, dada la generalización de las propiedades euclídeas del espacio R3 a Rn, numéricamente pasa algo que tiene propiedades geométricas similares.

Finalmente, el trazo del camino homotópico se muestra en la Figura 4, donde se graficó p en el eje horizontal y τ en el eje vertical. El camino homotópico se muestra en color azul y una línea ilustrativa de τ=1 en color rojo, esta línea sirve para detectar a simple vista los tres cruces que se mencionaron anteriormente.

Fuente: Elaboración propia.

Figura 4 Gráfico del camino homotópico para el caso ejemplo de dos ecuaciones. 

Caso de estudio 1

Consideremos el caso de estudio utilizado por ^{Malinen & Tanskanen (2010)}. Se trata de un sistema de dos ecuaciones algebraicas no lineales, ecuación 30, cuyas soluciones se encuentran situadas entre -4≤x1,x2≤4. Para este sistema se encuentran reportados un total de nueve vectores solución encontrados con el método de homotopía con un límite en su parámetro homotópico.

El programa de seguimiento homotópico hiperesférico (SHH) fue capaz de encontrar el mismo número de vectores solución (V_n), los cuales se muestran en la Tabla 3. Este resultado se obtuvo con un vector inicial x10,x20,τ0,r0=1,1,0,0.1, el tiempo de cómputo utilizado fue de 16 s.

Gracias a la parametrización del sistema homotópico con respecto a p, es posible graficar el camino homotópico, tal como se mostró en el caso ejemplo. Así pues, para este caso de estudio el camino homotópico se muestra en la Figura 5, donde ∆p recorrió el espacio de p tanto en sentido negativo como en positivo. El sentido negativo se encuentra marcado de color rojo, mientras que el sentido positivo se muestra de color azul. La línea de color verde representa τ=1, de manera que se puedan detectar a simple vista los nueve cruces que representan los nueve vectores solución.

Fuente: Elaboración propia.

Figura 5 Camino homotópico del caso de estudio 1. 

De acuerdo con la definición del método de SHH como un método de convergencia global, que tiene como premisa encontrar la mayor cantidad de vectores solución con una sola aproximación inicial, los resultados presentados anteriormente son los esperados comparados con otros métodos de convergencia local.

A continuación, se muestra un breve análisis de este caso de estudio realizado por un método numérico basado en el concepto de continuación homotópica (BCH) (^{Martín, 2013}). Dicho método deforma el sistema original en un sistema homotópico de tal forma que se transforma en un problema de valor inicial; posteriormente, se utiliza el algoritmo de Runge-Kutta de 4° orden para hacer la búsqueda del vector solución en el intervalo [0<τ<1] (Apéndice C). Para complementar la evaluación del caso de estudio, también se utilizó el método de Newton-Raphson para sistemas de ecuaciones no lineales.

Los resultados de los dos métodos anteriores, en conjunto con los resultados del método de SHH, se muestran en la Tabla 4. Los tres métodos fueron evaluados con diferentes vectores iniciales en un equipo computacional con procesador Intel® Xeon® E5-2620 V3 @2.4GHz con 64GB de RAM. Los resultados incluyen los vectores solución encontrados por cada método, numerados como se muestra en la Tabla 3, así como el tiempo de cómputo utilizado durante la evaluación y la proximidad del vector solución obtenido con la solución en fx=0.

Caso de estudio 2

El siguiente caso de estudio es un sistema de tres ecuaciones algebraicas no lineales con cuatro vectores solución reportados en la literatura (^{Kuno & Seader, 1988}; ^{Rahimian et al., 2011b}). Dichos vectores fueron encontrados por los autores utilizando un método de homotopía global de punto fijo.

En el método de SHH se utilizó un vector inicial de: x10,x20,x30,τ0,r0=[0,0,0,0,0.1]. El camino homotópico se muestra en la Figura 6; sin embargo, debido a la cercanía de los cruces, es imposible detectarlos a simple vista, por lo que se recurrió a varios acercamientos.

Fuente: Elaboración propia.

Figura 6 Camino homotópico del caso de estudio 2. 

La Figura 7a muestra un acercamiento a la sección del camino homotópico en sentido negativo que atraviesa la línea transversal de τ=1, mientras que en la Figura 7b se muestra un acercamiento a la sección del sentido positivo que atraviesa la misma línea. Cada una de estas líneas representan dos cruces, uno de forma ascendente y otro de forma descendente.

Fuente: Elaboración propia.

Figura 7 a) Sección del camino homotópico en sentido negativo; b) Sección del camino homotópico en sentido positivo. 

Finalmente, la Figura 8 muestra un acercamiento mucho más cercano a los cruces en la sección en sentido positivo del camino homotópico. Esto es para demostrar cómo el camino homotópico cruza ascendentemente τ=1 para luego descender generando otro cruce. Tomando en cuenta esto, hay un total de ocho cruces del camino homotópico con τ=1, cuatro más que los reportados en la literatura.

Fuente: Elaboración propia.

Figura 8 Sección en sentido positivo del camino homotópico que muestra dos cruces en τ=1. 

Los vectores solución correspondientes a cada uno de los cruces antes mencionados se muestran en la Tabla 5, señalando cuáles fueron los vectores solución adicionales encontrados. Este proceso llevó un tiempo de cómputo de 25 s.

Caso de estudio 3

En un caso de estudio aplicado en el área de ingeniería química, considere el evaporador de triple efecto mostrado en la Figura 9. Dicho evaporador tiene por objetivo realizar la concentración de una solución de azúcar del 10% al 50% en peso. La alimentación de la primera etapa usa un vapor saturado a 280 °F y produce líquido a 100 °F a un flujo másico de 50000 lb/h. La tercera etapa debe operarse a una presión absoluta, de manera que el punto de ebullición de la solución sea de 125 °F. Descarte el efecto del incremento del punto de ebullición en las variaciones del calor específico y el calor latente con respecto a la temperatura y la concentración. Determine el área de transferencia para cada etapa (se consideran iguales), las temperaturas (T1,T2), los flujos (L1,L2,L3), las componentes (x1,x2) y el flujo V0 del vapor de alimentación. Tome Cp=1 BTU/lb°F y un calor latente λ=1000 BTU/lb con coeficientes de transferencia de calor U1=500 BTU/ft2h °F, U2=300 BTU/ft2h °F y U3=200 BTU/ft2h °F (^{Oliveros-Muñoz & Jiménez-Islas, 2013}) .

Fuente: ^{Oliveros-Muñoz & Jiménez-Islas (2013)}.

Figura 9 Esquema de un evaporador de triple efecto. 

De un balance de materia global se obtiene que L₃ = 10000 lb/h, por lo que, planteando balances macroscópicos de materia, energía y la transferencia de calor regida por la ley de enfriamiento de Newton y sustituyendo los valores conocidos, resulta el siguiente sistema de ecuaciones no lineales:

Para realizar el SHH se utilizó el vector inicial x10,x20,x30,x40,x50,x60,τ0,r0=0,0,10,10,100,100,0,0.1 con un tiempo de cómputo de 20 s.

El camino homotópico se muestra en la Figura 10. Se puede observar un total de cinco cruces, sin embargo, cuando se refinan, se muestra que solo tres de ellos son vectores solución que no se repiten. Dichos vectores solución se muestran en la Tabla 6. Para una mejor visualización del camino homotópico se consideró un factor de multiplicación en p de 1X1025. El vector solución factible termodinámicamente es V₃, que no es fácil de encontrar con otros métodos como el de Newton-Raphson.

Fuente: Elaboración propia.

Figura 10 Camino homotópico para el caso de estudio 3. 

Al igual que el caso de estudio 1, este caso de estudio se evaluó con el método BCH y Newton-Raphson multivariable en igualdad de condiciones, los resultados se muestran en la Tabla 7. Dichos resultados muestran los vectores iniciales utilizados, los vectores solución encontrados por cada método, numerados de acuerdo con la Tabla 6, mientras que la precisión de cada vector solución con respecto a fx=0 se muestra en la parte inferior de la Tabla.

Caso de estudio 4

El siguiente caso describe el problema cinemático que presenta un manipulador paralelo 3-RRS. Este sistema consta de dos plataformas triangulares, la plataforma inferior es fija mientras que la superior es móvil; ambas están unidas por tres brazos con tres articulaciones, uno en cada una de las plataformas y otro en medio del brazo de unión, el sistema se muestra en la Figura 11.

Fuente: (^{Gallardo-Alvarado, 2019}).

Figura 11 Esquema de manipulador paralelo 3-RRS. 

El problema consiste en encontrar las posiciones alcanzables de la plataforma móvil que cumple con un conjunto prescrito de coordenadas generalizadas, para este caso son q1=120°,  q2=135° y q3=95°. Dichas coordenadas están relacionadas a la plataforma fija del manipulador paralelo, por lo tanto, el sistema se encuentra descrito por el siguiente sistema de ecuaciones algebraicas no lineales (^{Gallardo-Alvarado, 2019}):

Para resolver el sistema por el método de SHH se utilizó un vector inicial xi0,τ0,r0=0.1,0,2.5,0.5,1.5,1,0.5,-0.1,1.5,0,0.1. Los vectores solución encontrados tomaron un tiempo de cómputo de 1.5156 s, fueron un total de seis vectores, los cuales se presentan en la Tabla 8 y coinciden con el número reportado en la literatura.

En la Figura 12a se muestra el camino homotópico del problema. Al igual que en problemas anteriores, se presenta el camino en sentido positivo en color azul y en sentido negativo en color rojo, además de una línea τ=1 donde se puede ver una gran cantidad de cruces. La Figura 12b muestra un acercamiento en una sección del camino homotópico en sentido positivo, de esta manera se puede ver con más claridad que, a pesar de tratarse de una gran cantidad de cruces, todo se reduce a una repetición periódica de seis cruces, cuya longitud de onda va en aumento conforme el camino se desplaza a la derecha o a la izquierda. Dicho de otra manera, se trata de los mismos vectores solución presentados con anterioridad y repetidos con frecuencia uno detrás de otro.

Fuente: Elaboración propia.

Figura 12 a) Camino homotópico del caso de estudio 4; b) Acercamiento a una sección del camino homotópico en sentido positivo. 

Al igual que en los casos de estudio 1 y 3, este caso de estudio se evaluó con el método BCH y Newton-Raphson multivariable en igualdad de condiciones, los resultados se muestran en la Tabla 9. En dicha Tabla se muestra el vector inicial, los vectores solución encontrados por cada método, numerados de acuerdo con la Tabla 8, mientras que la precisión de cada vector solución con respecto a fx=0 se muestra en la parte inferior de la Tabla.

Caso de estudio 5

A continuación, se presenta una reacción de esterificación catalizada por ácido:

la cual se lleva a cabo continuamente en dos reactores de tanque de agitado continuo (CSTR) conectados en serie (Figura 13) (^{Seader et al., 1990}).

Fuente: (^{Seader et al., 1990}).

Figura 13 Diagrama de flujo de producción de ftalato de dibutilo. 

El equilibrio de vapor / líquido se establece en el reactor con ventilación de vapor. El primer CSTR es alimentado con una cantidad conocida de ácido, y una cantidad adicional de ácido ingresa al segundo CSTR para lograr una concentración de ácido específica en su efluente. Las ecuaciones de modelado consisten en cuatro ecuaciones de equilibrio de fase y tres ecuaciones de equilibrio de materia, la versión polinómica simplificada de dichas ecuaciones se muestra a continuación (^{Seader et al., 1990}):