© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
Rev. Acad. Canar. Cienc., XIII (Núms. 1-2-3), 137-152 (2001) (publicado en Julio de 2002)
Implementación de un Algoritmo Genético en
los Cálculos de Equilibrio Líquido-Vapor
Fernando Espiau, Juan Ortega, Pedro Saavedra, Carmelo González
Universidad de Las Palmas de Gran Canaria
Resumen
En este trabajo se utiliza una expresión polinómica, con coeficientes dependientes de la temperatura
y función de la llamada fracción activa de los constituyentes de una mezcla, z(xJ, que
permite correlacionar simultáneamente las cantidades termodinámicas de un equilibrio de fases,
como son la función de energía de exceso de Gibbs y sus coeficientes de actividad, y la
entalpía de exceso. La correlación se lleva a cabo utilizando el principio de mínimos-cuadrados
e implementando un algoritmo genético para la optimización de la.función objetivo que resulta.
Se realiza una aplicación del procedimiento a un conjunto de sistemas binarios reales, comparando
los resultados con los obtenidos por otro método clásico conocido. La valoración de la
aplicación es excelente y produce buenas estimaciones de las propiedades de mezcla (Gt/RT) y
(ff!R) que caracterizan al equilibrio líquido-vapor de mezclas binarias elegidas.
Se concluye que la expresión que mejor correlaciona la función de energía de exceso de
Gibbs para un sistema binario, en.función de la concentración de uno de los componentes y de
la temperatura es:
Palabras clave: Algoritmo genético, equilibrio líquido-vapor
Summary
A po/ynomia/ expression with temperature-dependent coefficients and a function for the socalled
active fraction of the mixture components, z(xJ, was used to corre/ate such phase equilibrium
thermodynamic quantities as the Gibbs excess energy .function, the activity coefficients
for that .function, and the excess enthalpy simultaneously. The correlation was based on the
principie of /east-squares, and a genetic a/gorithm was imp/emented to optimize the resu/ting
objective .function. The procedure was performed on a set of real binary systems, and the resu/ts
were compared with those obtained using a well-known conventiona/ method. Excellent results
were achieved using the procedure, which yie/ded good estima/es of the mixing properties
(GEIRT) and (ff!R) characterizing the vapor-liquid equilibria for the binary mixtures employed.
In conclusion, the expression that yie/ded the best correlation of the Gibbs excess energy
.function for binary systems, as a function of the concentration of one of the system components
and temperature is:
~;(x;,T}=z(x;J(i-z(x;){( A;, +A02 )+(A;, +A22 } 1 (x;)]
Key words: Genetic algorithm, vapor-liquid equilibrium
137
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
l. Introducción.
En este trabajo se analiza la utilidad de un método o procedimiento que proporcione
valores adecuados de las cantidades termodinámicas macroscópicas más características
de sistemas de más de un componente en equilibrio de dos fases, mediante el
uso de expresiones empíricas y con datos experimentales previamente conseguidos. Se
discute en primer lugar un formalismo de representación de las propiedades termodinámicas
de mezcla utilizando un modelo termodinámico-matemático que no presenta ninguna
relación con la Física Molecular o la Mecánica Estadística. Con ello no se desprecian
los esfuerzos que actualmente se llevan a cabo por los teóricos en esas áreas, para
explicar a través de un modelo el comportamiento de los sistemas fluidos, sino que,
aceptando el hecho de no existir una teoría satisfactoria, con una buena base científica
que proporcione buenas predicciones y/o correlaciones de datos termodinámicos de
mezclas líquidas, se opta también por desarrollar un método empírico como forma de
incentivar el refinamiento de los modelos que son de aplicación en esta área de trabajo.
El interés de los datos experimentales directos de propiedades termodinámicas
y/o de sus correlaciones está en los distintos usos que de ellos puede hacerse. Entre
otros destacan, la posibilidad de su aplicación en los cálculos de diseño de operaciones
de separación en Plantas Químicas, la verificación de modelos teóricos, o bien la de
encauzar la formulación de aquellas teorías que intentan estimar las propiedades de
mezclas a partir de las de los componentes puros, aspecto donde se desea incluir la aportación
de este trabajo. Un aspecto importante y básico es lograr un adecuado tratamiento
de datos experimentales con modelos que permitan estimar el mayor número de propiedades
de manera precisa. Lo ideal sería adoptar en una primera etapa, un modelo simple
con un cierto fundamento teórico del problema y adaptable a distintas situaciones. En
una segunda fase, ya durante el procesamiento de datos, lograr su inserción en un adecuado
algoritmo haciéndolo convergente. Todo ello constituirá una procedimiento que
aquí se trata de establecer incluyendo las dos etapas mencionadas, con la propuesta de
un modelo sencillo y su implementación en un algoritmo genético, de uso común en
ingeniería, y cuyo principio fundamental es la convergencia en solución única. Tales
algoritmos (Holland, 1992; Gen y Cheng, 1997) tiene como objetivo su aproximación a
óptimos absolutos del modelo considerado, permitiendo trabajar con un elevado número
de variables sin que ello suponga alto coste computacional.
Entre los diferentes procesos de separación a los que puede aplicarse el método
que se elabora y se presenta en este trabajo, cabe mencionar el de la destilación y/o rectificación
industrial, que es una operación de la Ingeniería Química que siempre opera a
presión constante. Los estudios sobre equilibrios de fases, y sobre todo los de equilibrio
líquido-líquido y líquido-vapor, son de gran ayuda para lograr un diseño óptimo de dichas
operaciones. Específicamente, los datos de equilibrio isobáricos, o bien los de unas
correlaciones precisas, son importantes para lograr el mejor diseño de los sistemas que
operan en la destilación.
Aquí se llevará a cabo una aplicación a varios casos reales de sistemas binarios
formados por un éster de propilo con los cuatro primeros n-alcanoles, con el fin de verificar
la validez de las expresiones y la metodología que se propone en aras de conseguir
el mejor tratamiento matemático de las cantidades que caracterizan termodinámicamente
al equilibrio líquido-vapor de sistemas binarios. Los resultados de este trabajo se
compararán con los obtenidos por otro método de uso común en esa área.
138
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
2. Planteamiento termodinámico-matemático
Las propiedades termodinámicas de exceso con un mayor interés, tanto bajo un
punto de vista práctico como teórico, en los estudios que se llevan a cabo sobre el comportamiento
de sistemas fluidos, son la función de energía de Gibbs y la entalpía de
exceso; la entropía puede calcularse a partir de aquellas. Uno de los aspectos que ha
recibido mayor atención en el campo de la Termodinámica de soluciones, es el estudio
del equilibrio líquido-vapor, por eso la literatura muestra abundantes datos para multitud
de sistemas binarios y multicomponentes, siendo muchos conjuntos de esos valores
presentados a presión constante y más concretamente a presión atmosférica. Si bien los
datos experimentales a temperatura constante o isotérmicos son preferidos por la Termodinámica
teórica, los datos a presión constante, normalmente más exactos, proporcionan
una excelente fuente de estudio cuando también se conocen las entalpías de exceso,
como se verá en este artículo.
La ecuación fundamental que contiene a las cantidades que juegan un importante
papel en los procesos de equilibrio y que establece las condiciones para el comportamiento
de la fase líquida, surge a partir de un planteamiento sobre la dependencia de la
función de Gibbs con las otras variables que intervienen, como son la presión, la temperatura
y la cantidad o número de moles de las sustancias que constituyen el sistema; así,
si G=G(p, T,nJ, su diferencial produce:
dG=(ªG ) dr+(ªG) dp+ ¿(ªG) dn;
ar p,n ap . j an¡
1 7 ,n1 p ,T,nj.,,_1
(1)
y de aquí surge (como puede verse con más detalle en cualquier obra de Termodinámica
Química, Smith y Van Ness, 1987) una de las ecuaciones más importantes y generales
que pueden plantearse en los estudios de la Termodinámica del equilibrio entre fases:
HE VE
- -dr +-dp-"xdlny =0
Rr2 Rr '7- ' ' (2)
También, a la Termodinámica Química elemental pertenecen otras expresiones que son
utilizadas en este trabajo, como son las que relacionan la/unción de Gibbs con los coeficientes
de actividad,
o bien lny -- -G+E x (-ª(G-E-/ R-r- ))
' Rr 1 &:; 7.
,p
(3)
Un asunto relevante que debe considerarse en el tratamiento termodinámico del
equilibrio de fases es la dependencia con la temperatura de las propiedades de exceso.
Hasta ahora este aspecto ha recibido por los científicos menos atención de la necesaria.
Sin embargo, una vez se haya aceptado su necesidad, se deben intentar conseguir relaciones
termodinámicas más exactas que permitiría reducir considerablemente el esfuerzo
experimental. Dos son las relaciones básicas que pueden utilizarse para conseguir
dicho objetivo:
e: =[ 8(HE / R)]
R ar p,r,
(4)
139
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
-HE =T[accE 1 RT)J
RT 8T p ,x, ..
(5)
si bien, como se observa, ambas ecuaciones pueden fundirse en una sola. Para concretar
la utilidad de las mismas y sobre todo para poder aplicarlas en este trabajo, deben hacerse
algunas consideraciones que permita integrarlas con facilidad. Así, tanto para los
líquidos puros, e p, como para las mezclas, e:, la dependencia de la capacidad térmica
con la temperatura, no suele ser muy diferente de una función lineal de primer o segundo
grado. En el caso de la e: de un sistema fluido de más de un componente, a una
determinada concentración, dicha cantidad podría estar bien representada en un amplio
intervalo de temperaturas por una expresión cuadrática sencilla, como se expresa a continuación,
CE
_P =a+bT +cT2
R
(6)
si ahora en la ecuación (4) se sustituye la (6) y se integra respecto a la temperatura se
produce:
HE =aT+!z.T2 +E..T3 +!
R 2 3 I
(7)
donde la 11 es una constante de integración. Repitiendo el proceso con las ecuaciones (5)
y (7), se tiene,
GE =-alnT-Q_T_E..T2+11 +!
RT 2 6 T 2
donde h es otra constante que surge de esta última integración.
(8)
En principio, con las ecuaciones (6-8) todas las constantes podrían estar determinadas,
eso sí, dependiendo de la disponibilidad de datos experimentales de las cantidades
termodinámicas en juego. Una vía sería la de utilizar valores de equilibrio líquido-
vapor a varias temperaturas para proporcionar datos de GEIRT en función de T, lo que
daría lugar a varias situaciones de la ecuación (8) que permitirían calcular las constantes
a, b, e, 11 e h de una sola vez. Hecho esto y volviendo atrás, las diferentes ecuaciones
(6,7) para las propiedades de exceso como función de Testarían también definidas. Sin
embargo, una desventaja de este procedimiento es que las cantidades primarias, como
las ¡.¡E y las e: se lograrían a partir de, respectivamente, una primera y una segunda
diferenciación de la función de Gibbs, perdiéndose en cada proceso diferencial un orden
de magnitud. Esto se confirma porque, de las cinco constantes que se requieren para
GEIRT, solo cuatro son necesarias para las entalpías HEIR, y tres para las capacidades térmicas
de exceso, e: IR. Otro aspecto negativo de este procedimiento es que, para calcular
valores de la función Gº/RT a partir de datos de equilibrio se requiere de otras cantidades,
a menudo imprecisas, como son los coeficientes de actividad, ver ecuación (3).
Este hecho ocurre a pesar de tener la firme convicción de la alta precisión de los datos
que se consiguen en la experimentación del equilibrio entre fases.
Otra forma de proceder para calcular las constantes consiste en introducir varias
medidas de capacidades térmicas en función de T. De esta forma se tendrían valores
140
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
para a, b y e, y luego con la ecuación (7) y un valor de I-f a una sola temperatura, se
calcularía un valor para la constante de integración / 1• Por último, un valor de GEIRT vs T
permite determinar la segunda constante de integración h El mayor inconveniente de
este segundo procedimiento es que se necesitan datos de tres cantidades termodinámicas
diferentes, lo que conlleva tener dispuestos tres sistemas experimentales diferentes. En
descarga de este método puede decirse, por la evidencia experimental contrastada por
nosotros, que al menos los valores de e; IR determinados por diferenciación de los datos
de entalpías, son tan exactos como los que se obtienen a partir de la experimentación
directa, lo cual, unido al interés secundario de las capacidades térmicas en exceso en
relación a las otras cantidades consideradas, haría reducir la labor experimental a solo
dos de las cantidades. Este comentario revela bastante bien la necesidad de poseer datos
precisos de las entalpías de exceso en función de la temperatura.
Pueden admitirse otras simplificaciones sucesivas de las ecuaciones (6-8) que
dan lugar a que e; IR sea función lineal de T, con c=O, o incluso independiente de T,
con b=c=O. Con estas dos consideraciones particulares, las relaciones anteriores se convierten,
respectivamente, en las siguientes expresiones.
CE
_P =a+bT
R
HE b 2
, -=aT +- T +! . R 2 1
GE b f ¡
- =- alnT--T+- +I
RT 2 T 2
GE /
-=-alnT +-1.+J
RT T 2
(9)
(10)
Resulta conveniente aclarar que cuando se trabaja en un corto intervalo de tempe~aturas,
la consideración de independencia entre e; IR y T es acertada en la práctica, lo cual
equivale a considerar que 1-FIR es lineal en T, como se recoge en la segunda expresión de
(1 O). En los conjuntos de ecuaciones (9) y (1 O), los parámetros pueden calcularse a partir
de valores apropiados de I-f IR y de GEIRT en función de T, como se indicó antes.
Si se continua con este mismo razonamiento sistemático de simplificación, se
llegaría a los casos más simples y extremos, como aquel en el cual se verifica que
a=b=c=O en la ecuación (6), o el más sencillo, donde a=O como caso particular de la
primera expresión de (1 O), ocasionando también que las entalpías sean independientes
de la temperatura (lo que resulta bastante usual en la práctica por la falta de datos reales
que permita establecer dicha dependencia). Esta última consideración hace que la segunda
y tercera expresión de (1 O) reduzcan, respectivamente, a:
HE
- R = / I
GE=!_¡_+/
RT T 2
(11)
(12)
Incluso, los casos planteados en el conjunto de ecuaciones (9) y (1 O) pueden ser
fusionados en (11) y (12), considerando que algunas de las cantidades presentadas en
esas ecuaciones, como la correspondiente a lnT, pueda ser desarrollada en serie de potencias
y sustituida por términos en T, como se verá en la sección 4. Todo esto forma
141
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
parte de un planteamiento formal que se basa en la utilización como Ciencia herramienta
de la Termodinámica, pero que debe adaptarse para la consecución de un modelo
empírico que tenga un sentido práctico y útil para cumplir con el objetivo que se pretende,
simular adecuadamente los valores más característicos de los equilibrios entre fases,
y particularmente del equilibrio líquido-vapor.
El planteamiento realizado hasta ahora, y donde puede captarse un cierto grado
de conveniencia, se ha desarrollado en base a considerar únicamente la dependencia con
la temperatura de las cantidades adimensionales de la entalpía de mezcla y de la función
de Gibbs. Sin embargo, la realidad es que dichas funciones también dependen de otras
variables, que deben ser consideradas para encontrar las expresiones adecuadas de correlación,
estas deben ser también tenidas en cuenta. Así, es trascendental incluir la dependencia
con la composición x; de un sistema multicomponente. Para completar este
estudio, se asume por ejemplo que la función de Gibbs en su forma adimensional GEIRT
para un sistema binario, queda bien representada a una determinada temperatura, por
una ecuación polinómica que fue presentada anteriormente (Ortega y Alcalde, 1992)
para correlacionar las cantidades de exceso con la composición y cuya expresión truncada
en el tercer término, viene dada por
GE 2
-=Z1Z2( Áo + Á¡Z¡ + A2Z1) RT
(13)
siendo z2=l-z1 y z1=x¡/(x1+kx2), son las llamadas fracciones activas de cada uno de los
componentes en la mezcla. El método más directo de integrar (13) teniendo en cuenta el
formalismo termodinámico recogido hasta ahora, es expresar cada parámetro A; como
función de la temperatura T, de acuerdo a como se plantea en la (12) y que para esta
ocasión podría escribirse como: A;=A;c/T+A¡¡, quedando ahora la expresión (13) algo
más compleja por la inclusión de la dependencia de lafanción de Gibbs con la temperatura
y la concentración,
m
<> z(J-zJ"L.A,z; (14)
i=O
Esta expresión, en la que los A¡¡ son constantes, proporciona una aproximación
más real de la/unción de Gibbs de mezcla, siempre que su aplicación se restrinja a los
casos que comprendan intervalos no muy amplios de temperatura. Recordando ahora la
relación (5), podrían reproducirse las entalpías de exceso. Así, derivando la ecuación
(14) se obtiene:
En el modelo que se plantea, el parámetro k influye más bien en la forma de la
curva que en su tamaño, por lo que puede considerarse independiente de la temperatura,
tomando la ecuación (15) la forma simplificada de la página siguiente. Una explicación
más detallada del significado de k puede verse en Pacheco y Ortega (2000).
142
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
(16)
sustituyendo las derivadas de A;, de acuerdo con la definición establecida antes, se tiene:
--H-=E z (1-z ) (A--2 .!._+_A!_ !_z +_Al. !_z 2 )
RT I I T T I T J
(17)
dando lugar a que los coeficientes A;¡/T sean idénticos a los que se conseguirían en un
ajuste o en una correlación directa de datos, con las entalpías presentadas en su forma
adimensional, ff!RT vs z, con una ecuación similar a la (13) pero, donde ahora, la expresión
es,
E
H =zJ 1-z, J(A; +A; z1+A;z:)
RT
siendo: A,1 = ATn (18)
o bien directamente, cuando se correlaciona a las ff!R vs z. Lo esencial de este procedimiento
de cálculo estriba en que una vez estimados los A/ con la (18), puedan determinarse
fácilmente los parámetros Au de la ecuación (14) y (17), debiendo luego calcularse
solo los Ai2 con un adecuado procedimiento de regresión de los datos GERTylo los y;
frente a la composición y la temperatura, lo que supone un motivo de comentario de este
trabajo.
Centrémonos mínimamente en el hecho experimental. Los pasos a seguir antes
del tratamiento de datos y el uso de las ecuaciones planteadas a lo largo de esta sección
son los siguientes: en primer lugar, es habitual en los estudios sobre mezclas líquidas
utilizar las cantidades de if y ~como propiedades de exceso y las de x1, y¡, T, y p en
los equilibrios líquido-vapor, ya sean isotérmicos y/o isobáricos; a continuación, a partir
de los datos de equilibrio, se calculan los coeficientes de actividad y los valores de función
de energía de exceso de Gibbs utilizando la expresión (3).
Con todo ello, se comprende la importancia de la expresión (2), que puede considerarse
como un auténtico modelo termodinámico-matemático para verificar, como si
de un test se tratase, la consistencia del equilibrio entre fases, y que puede aplicarse a
distintas condiciones de operación. Para el caso concreto en el que se estudien los equilibrios
líquido-vapor isobáricos, deben emplearse tres conjuntos de datos (ff!RT vs x1),
(lny¡ y Inri vs x1) y a partir de estos últimos los (GEIRT vs x 1), para conseguir un adecuado
tratamiento de las cantidades termodinámicas del equilibrio. No obstante, la elección
más adecuada dependerá del método de regresión como se verá en el siguiente apartado.
3. Procesamiento de los datos de equilibrio
Del planteamiento llevado a cabo en la sección anterior, se deduce que la Termodinámica
constituye una buena herramienta de la Ingeniería Química para calcular
con precisión las cantidades fisicoquímicas que establecen las condiciones de equilibrio
de sistemas fluidos. Con el fin de obtener interpolaciones y/o extrapolaciones de datos
en condiciones diferentes a aquellas que se consiguen en la experimentación directa, se
necesita un modelo que permita la estimación que se pretende. En dicho proceso de estimación,
que actuaría de forma muy similar al de otros campos o áreas de la Ciencia,
deberían seguirse los siguientes pasos o etapas:
143
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
(a) El desarrollo y la propuesta de un modelo que explique razonablemente el fenómeno
observado.
(b) Disponer de un buen número de medidas experimentales del fenómeno.
( c) Estimar y validar el modelo.
Por la orientación que tiene este trabajo se abordan especialmente los dos últimos
pasos. Además, en la propuesta del modelo, el investigador deberá considerar los errores
aleatorios y sistemáticos, así como aquellos derivados de las hipótesis exigidas al
mismo. Para la estimación de los coeficientes del modelo planteado se necesita optimizar
una función objetivo que puede resultar no-convexa y por tanto tener varios extremos
locales, además de tener un elevado número de variables. Tal función objetivo recogerá
las discrepancias entre las condiciones teóricas impuestas por el investigador,
plasmadas por el modelo, y las experimentales.
En este trabajo deben estimarse los coeficientes del modelo que se propone para la
función de energía de Gibbs, y que resultan de las distintas hipótesis que pueden adoptarse
para la capacidad térmica de exceso, utilizando el principio de mínimos cuadrados,
Draper y Smith (1981), como alternativa a otros métodos ya tratados en la literatura.
Este procedimiento de estimación es óptimo cuando los errores inherentes a la toma de
datos se distribuyen normalmente, son independientes y tienen la misma varianza. Se
utilizan conjuntos de valores, tal y como se han presentado en el último párrafo de la
sección anterior, constituidos por datos experimentales procedentes de la experimentación
directa con sistemas reales, para los que se suponen las condiciones que validan el
procedimiento de estimación utilizado.
Para llevar a cabo la aplicación concreta que se pretende en este trabajo, tal como
se comentó en la §2, para el tratamiento de datos de equilibrio líquido-vapor, las
magnitudes de temperatura, presión y concentración de las fases se toman durante un
proceso experimental, de forma directa las dos primeras e indirecta la última, luego, a
partir de todos esos valores se calculan los coeficientes de actividad y;, y por último, los
valores de la función de Gibbs de exceso. Las entalpías de exceso se miden directamente,
de manera independiente por calorimetría, estando relacionadas estas dos últimas
cantidades a través de la expresión (5).
En resumen, se tienen los siguientes conjuntos de datos que definen los estados de
equilibrio de un sistema fluido binario, por un lado, VT1 , x1.1 , ln(r 1.), ln(r1•1 ), j:l, ... ,n} y
por otro, {x,,,(H E/R); r:l, ... ,q} donde los ln(y¡) y /n(ri.J son los logaritmos de los coeficientes
de actividad obtenidos a la temperatura T¡ y a una concentración del compuesto
1, XJJ, mientras que las cantidades adimensionales (F!R), corresponden a las entalpías
de exceso medidas a una concentración del compuesto de referencia 1, x1,, , generalmente
diferente a la de los datos de equilibrio, y a una temperatura constante de T=298, 15
K. Los segundos subíndices se corresponden con el número de la medición que se lleva
a cabo para las dos magnitudes descritas.
El empleo del principio de mínimos cuadrados plantea una función objetivo que
refleja las discrepancias entre los valores experimentales y los valores estimados por el
modelo en un determinado estado de equilibrio, es decir con igualdad de concentración
y temperatura. Dicha función objetivo, FO, se plantea como:
(19)
144
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
cuyas variables son los coeficientes del modelo que se tome para lafunción de Gibbs.
Los estimadores de dichos coeficientes serán aquellos valores que minimicen FO. Debido
a la gran cantidad de coeficientes que pueda incluir el modelo elegido, lo que conlleva
un elevado número de variables en FO y a la posibilidad que ésta sea no-convexa,
y que por lo tanto presente varios extremos locales, la utilización de métodos numéricos
clásicos para su minimización puede no ser muy adecuada. Además, el empleo de algunos
de esos métodos supone el manejo de complejos sistemas de ecuaciones no-lineales
que complicarían la tarea de la estimación de los coeficientes. Para salvar estos inconvenientes,
se optó en este trabajo por la implementación de un algoritmo genético, cuyos
detalles se recogen en el apéndice.
El procedimiento aquí presentado para estimar los coeficientes del modelo adoptado
para la función de Gibbs presenta la ventaja de no llevarse a cabo en dos etapas,
como ocurre con el procedimiento descrito en la presentación de las ecuaciones (14-18),
en el que, en una primera se estiman algunos coeficientes, entre ellos el de k, considerando
únicamente las !-F!R vs x1, y que pasan como valores fijos a una segunda correlación
de los datos de GEIRT y los y; frente a la composición y a la temperatura; en esta
segunda regresión se obtiene la estimación del resto de los coeficientes y se recalcula
otro valor de k, generalmente con un valor distinto al conseguido en la primera fase. Por
el contrario, la optimización de FO arroja unas estimaciones válidas conjuntamente para
las tres cantidades consideradas, y además, en una sola etapa.
Se observa que en el planteamiento de la FO no se consideran datos la función
adimensional de Gibbs, dado que los valores de dicha función no se miden de manera
directa ni son totalmente independientes respecto a los demás, sino que se determinan a
partir de los coeficientes de actividad y por ello su inclusión en (19) no aportaría una
información de interés adicional al algoritmo y tampoco a la expresión de la FO.
4. Aplicación del modelo propuesto
Como aplicación del modelo termodinámico-matemático que se propone y del
procedimiento de regresión que se ha desarrollado para este trabajo, con el fin de verificar
su utilidad en el tratamiento de datos experimentales de equilibrio líquido-vapor, se
efectúa una aplicación práctica sobre datos experimentales recientes de sistemas binarios
formados por ésteres de propilo con a/cano/es normales, que ya han sido aceptados
en la bibliografia internacional. En concreto se han escogido datos de equilibrio líquidovapor
isobárico y entalpías de exceso para cuatro mezclas binarias del etanoato de propilo
con n-alcanoles (desde el metano/ al butano!). La bondad del procedimiento utilizado
se llevó a cabo comparando los resultados obtenidos con los de otros casos en los
que se utiliza un método/algoritmo clásico, analizando las ventajas e inconvenientes de
cada uno.
En un trabajo reciente, Ortega y col. (2000) utilizaron una ecuación base, similar
a la (13), para la correlación de datos de equilibrio líquido-vapor y propiedades de mezclas
binarias, pero conteniendo el producto x1x2=x(l-x) en lugar de z1z2=.z.(1-z) ya que
este último generaliza al primero, siendo idénticos únicamente si k= 1. Sin embargo, el
modelo matemático planteado en (13) presenta una justificación teórica, ver Pacheco y
Ortega (pendiente) cuyos fundamentos necesitan ser contrastados a partir de aplicaciones
concretas. Por ello nuestro interés en observar aquí su capacidad para correlacionar
las cantidades termodinámicas de equilibrios de fases y en particular las del equilibrio
líquido-vapor, utilizando un método adecuado de regresión.
145
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
La estimación de los coeficientes del modelo más generalizado, expresado en
(14), considerando los comentarios anteriores, es análoga a la utilizada por Ortega y col.
(2000), en el que se utilizó el procedimiento clásico de regresión de mínimos-cuadrados
propuesto por Marquardt (1963) para funciones no-lineales. Los resultados, junto a la
bondad de las estimaciones de las cantidades termodinámicas que se expresan para cada
uno de los sistemas tratados, se recogen en la tabla 1. Para estos casos, los coeficientes
se derivan del ajuste de las Fl!R vs x1, los A;1 se mantienen fijos en la posterior correlación
simultánea de la función de Gibbs, GEIRT y los coeficientes de actividad y; con la
composición de la mezcla, de la que resultarán las estimaciones de los Ai2, reajustándose
en una segunda parte del tratamiento el valor de k que ocasiona la mínima desviación
del conjunto. De esta forma se consiguen unos resultados que constituyen las dos primeras
columnas de la tabla 1, y donde se observa que la ecuación que utiliza el producto
x1x2 (columna 1) mejora la estimación de las FI, pero no la de las otras cantidades. Sin
embargo y como contrapartida, presenta bastante disparidad en los valores del parámetro
k y en los coeficientes, siendo bastante más uniformes los resultados de la columna
11, donde k oscila entre valores muy cercanos a la unidad. La estimación de las entalpías
en esta columna aparenta tener desviaciones elevadas, sin embargo, el error porcentuales
de la correlación es inferior al 8%.
En las figuras 1 a 4 se han representado los resultados de estas dos primeras
aplicaciones, observándose que, al menos cualitativamente, existe una buena representación
de los datos de equilibrio líquido-vapor isobáricos. Incluso, en la figura 1, para la
mezcla etanoato de propilo+metanol, con el modelo de la columna 11 se consigue la
mejor estimación de las GEIRT, aunque la peor de las FI. En los otros tres sistemas binarios,
la estimación de las entalpías de exceso con la forma del modelo considerado en 11
tampoco mejora a los que se consiguen en las columnas 1, III y IV.
En este trabajo se desea verificar la utilidad del método que se ha desarrollado
en las secciones 2 y 3 para el procesamiento de los datos de equilibrio de los sistemas
anteriores, pero empleando el algoritmo genético y la función objetivo descrita en (19),
para estimar los coeficientes del mismo. Las hipótesis que sobre la capacidad térmica de
exceso se han considerado son, la de e; =0, y la más real de e; =a+bT , las cuales llevan
a apareadas las relaciones (11,12) y (9), respectivamente. En la tercera expresión de (9)
se puede tener en cuenta el desarrollo en serie de potencias de (J'-h) para h>O, del logaritmo
neperiano lnT, el cual truncado en el primer término daría lugar a que la ecuación
(9) quedara como sigue:
-G=E - ªT- a1/1,n h- 1) - -bT +1-1 +I = - (-ª+ b-) T+-JI +I -a~nh-1)
RT h 2 T 1 h 2 T 1
(20)
el modelo (13), ajustado a la hipótesis de c ; =a+bT bajo estas consideraciones resulta:
(21)
O sea, dependiendo de cómo se considere la dependencia de la capacidad térmica
de exceso de una mezcla binaria con la temperatura, se tienen dos expresiones para la
función energética de exceso de Gibbs, (14) y (21), que se diferencian en la forma de los
146
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
Tabla 1. Valores obtenidos en la aplicación a las mezclas binarias (etanoato de propilo+n-a/canoles) para
los coeficientes del modelo (14), considerando dos procedimientos de regresión.
Algoritmo regresión no-lineal (Marquardt, 1963) Algoritmo genético
Condición
Ecuación
Mezcla
~"' e
·¡":¡
oc
"o u
Mezcla
~"' e
·¡":¡
oc
"o u
Mezcla
e:/ R=O
1.- x1x2(A0+ A1z) H.- z1z2(A0+A1z+AiJ) IH.- z1z2(A0+AiJ)
Ao1=602.97
A02=-0.625
A 11=-205.42
A12=-0. 189
k=2.317
s(yJ=0.050
s(GEIRT)=0.016
s(JF)=75.5
Ao1=519.4
Aoi=-0.513
A11=166.27
A11=-0.670
k=0,401
s(yJ=0.053
s(GE/RT)=0.008
s(fF)=6.9
Ao1=637.71
Ao1=-I. I 30
A11=224.28
A11=-0. 724
k=2.092
s (yJ=0.009
s(GE IRT) =O. 003
s(fF)=8.9
Au
A.=-+A 2 ' T '
Au A.=-+A 2 ' T '
x1Etanoato de propilo + x2Metanol
A01=387.52
Aoi=-0.093
A11=220.75
A12=-0.424
A21=0.024
An=-0.849
k=0.866
s(yJ=0.026
s(GE IRT) =O. O I 3
s(fF)=76.9
Ao1=436.07
Aoi=-0.244
A11=243.66
An=-1.284
k=0.983
s (yJ=0.029
s(GE/RT)=0.019
s(fF) = 11. I
x 1Etanoato de propilo + x1Etanol
Ao1=449.52
Aoi=-0.362
A11=357.32
A12=-l.137
A21=-0.593
A22=-0. 149
k=0.980
S(yJ=0.023
S(GEIRT) =O. 008
s(fF)=57.0
Ao1=577.32
Aoi=-0. 778
A21=205.26
An=-0.839
k=0,969
s (yJ=0.032
s(GEIRT)=0.009
s(fF)=20.6
x1Etanoato de propilo + x1Propan-l-ol
A01=555.47
Aoi=-0.893
An=/24.57
A11=-0.433
Au=-307.55
An=-0.867
k=I .030
s(yJ=0.008
s(GEIRT)=0.002
s(fF) = 79 8
147
A01=609. 73
Ao1=-I. I06
A11=351.82
An=-0.944
k=0,902
s (yJ=0.01 J
s(GEIRT)=0.002
s(fF)=/ 1.2
c : IR=a+bT
IV.- z1z2(A0+AiJ)
A01=-0.000
Ao1=405.80
Ao1=-0.055
A11=-0.002
A22=63.27
A13=-0.026
k=0.977
s (yJ=0.028
s(GEIRT)=0.018
s(fF) = 13. 3
Ao1=-0.00I
Ao1=442.53
Ao3=0.IOO
Au=-0.003
A22=-I0.26
A13=0.797
k=0.954
s (yJ=0.030
s(GEIRT)=0.009
s(fF)=2 l .9
Ao1=-0.00I
Ao1=420.53
A01=-0.103
A11=-0.003
An=281.37
A13=0.45/
k=0.796
s (yJ=0.013
s(GEIRT)=0.003
s(if)= 14. l
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
Tabla l. Continuación
Mezcla x 1Etanoato de propilo + x2Butan-J-ol
Ao1=710.87 Ao1=835.49 Ao1=729.27 A01=-0.00I
A02=-I.399 A02=-I.697 A02=- l .459 Ao2=282.50
A11=158.51 A 11=-228.87 Ao3=-0. I 37
Á¡¡=-0.468 A12=0.507 Au=-0.006
A21=127.04 A21=98.00 A22=580. 14
An=-0.328 A22=-0.340 A23=/ .372
k=2.609 k=I .120 k=l .016 k=0.486
s (rJ=0.014 s(rJ=0.014 s (rJ=0.014 s (rJ=0.014
s(GEIRT)=0.007 s(GEIRT) =0.007 s(GEIRT) =0.008 s(GEIRT)=0.007
s(¡¡E)=6.5 s(¡¡E)=W l s(¡¡E)=6 7 s(¡¡E)=349
los coeficientes. Por otro lado, si en ambas ecuaciones se desarrollan sus productos y se
expresan como polinomios en "z1" se observa una sobreparametrización innecesaria, al
resultar más parámetros que monomios en z1• Esta circunstancia se evita considerando
en el modelo las potencias pares de z1 en los sumandos de las referidas expresiones,
quedando definitivamente las expresiones finales truncadas en el segundo sumando y
con la siguiente forma.
(22)
(23)
Los resultados de estos modelos para cada una de las mezclas binarias seleccionadas,
se recogen en las columnas III y IV, respectivamente, de la tabla l. La primera
conclusión que se deduce de una detenida observación de la misma, es la poca significación
de los coeficientes de temperatura que se consiguen para el modelo (23), los A01
y A21, lo cual sugiere su eliminación práctica del mismo. Obsérvese que al llevar a cabo
tal supresión el modelo (23) se reduce a (22), lo que implica que el enriquecimiento que
podría ocasionar la relación entre la capacidad térmica de exceso y la temperatura no
redunda en un mejor ajuste. Por otra parte, se aprecia que esta última versión es precisamente
la que mejor estima las entalpías de exceso, siendo aceptables las correlaciones
para las demás cantidades que caracterizan a los equilibrios líquido-vapor de las mezclas
que se ha escogido de (etanoato de propilo+n-alcanoles). También se observa que
las estimaciones de los coeficientes del modelo (22) se comportan de una forma muy
homogénea en los distintos sistemas, no así las de los otros modelos, en los que hay una
mayor fluctuación entre sus coeficientes. Todo esto, unido al menor número de coeficientes,
y por tanto a una menor complejidad, hace que el modelo expresado en (22) sea
el más adecuado.
148
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
4. Conclusión
Se ha desarrollado un procedimiento para mejorar el tratamiento termodinámicomatemático
de los datos que caracterizan a los equilibrios líquido-vapor de las mezclas
binarias, utilizando una nueva forma polinómica para lafanción adimensional de Gibbs,
relacionada con la concentración de los compuestos de la mezcla, a través de la llamada
fracción activa, y de la temperatura. La expresión final más adecuada resultó ser
(24)
Se realizó la estimación de los coeficientes del modelo (24) con los datos de los
coeficientes de actividad y entalpías de mezcla, siguiendo el principio de mínimoscuadrados,
e implementando un algoritmo genético para la optimización del la función
objetivo definida en (19). La aplicación sobre un conjunto de cuatro mezclas binarias de
etanoato de propi/o+n-alcanoles (C1 a C4) dio excelentes resultados, por lo que el procedimiento
aquí desarrollado se propone para el tratamiento de datos de equilibrio entre
fases en futuros trabajos.
S. Bibliografía citada
Draper, N.R.; Smith, H. Applied Regression Analysis. John Wiley & Sons, New York
(1981).
Gen, M.; Cheng, R. Genetic Algorithms and Engineeering Design. John Wiley & Sons,
New York (1997).
Goldberg, D.E. Genetic Algorithms in Search, Optimiza/ion and Machine Learning.
Addison-Wesley, New York (1989).
Holland, J.H. Investigación y Ciencia, Septiembre: 38 (1992).
Holland, J.H. Adapta/ion in Natural and Artificial Systems. University of Michigan
Press, Ann Arbor (1975).
Marquardt, D.W. J Soc. Jndust. Appl. Math., 11, 431 (1963).
Michalewicz, z. Genetic Algorithm+Data Structure=Evolution Programs, 2"d ed.,
Springer-Verlag, New York (1994).
Ortega, J.; Alcalde, R. Fluid Phase Equilibria 71, 49 (1992).
Ortega, J.; González, C.; Peña, J. Galván, S. Fluid Phase Equilibria 170(1): 87 (2000).
Pacheco, J. M.; Ortega, J. Rev. Acad. Canaria de Ciencias, XII: 000 (2000).
Smith, J.M.; Van Ness, H.C. Introduction to Chemical Engineering Thermodynamics,
41h. Ed., McGraw-Hill Book Co., New York (1987).
5. Apéndice. Algoritmo genético
La utilización de métodos numéricos clásicos para minimizar la función objetivo
FO puede no ser adecuada, dado el elevado número de variables de ésta, a la posibilidad
de que existan múltiples extremos locales y a la necesidad de resolver complejos sistemas
de ecuaciones no-lineales. Por ello, en este trabajo se propone utilizar un algoritmo
genético cuyo principios surgen de los trabajos de Goldberg (1989) y Holland (1975),
desarrollados más tarde por otros artículos de Holland ( 1992, 1995), Michalewicz
(1994) y Gen y Cheng (1997), para llevar a cabo tal minimización en la estimación de
149
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
los coeficientes de los modelos (22) y (23). Un algoritmo genético tiene su fundamento
en los mecanismos genéticos y de selección de la naturaleza. Comienzan con un conjunto
aleatorio de vectores, los cromosomas de la primera generación, cuyas componentes,
llamadas genes, se corresponden con cada una de las variables de la FO, y que a su vez
son los coeficientes del modelo que se adopte para la función de Gibbs. Es por ello que
en el modelo planteado los cromosomas son de la forma x=(Ao1, Ao2, Au, A12, k) para el
modelo (22) y x=(Ao1, Ao2, Ao3, A21, A22, A23, k) para el (23). Estos estarán sometidos a
las operaciones "genéticas" de cruzamiento y mutación que dan lugar al conjunto de
cromosomas descendientes. A continuación, el conjunto de cromosomas progenitores y
descendientes será sometido a un proceso de selección para conformar la siguiente generación.
De esta forma, los cromosomas irán evolucionando en sucesivas generaciones,
sobreviviendo los buenos y pasando a la próxima generación, mientras los malos mueren.
Así, los genes convergerán hacia la solución del problema de minimización que se
quiere resolver. Estos valores de convergencia se considerarán como las estimaciones de
los coeficientes antes aludidos.
Los métodos numéricos tradicionales de optimización de funciones consisten en
la sucesión de pasos computacionales deterministas, basados en el gradiente o en derivadas
de orden mayor de la función a optimizar, que convergen en la solución óptima de
la misma. Dichos cálculos se aplican a un solo valor del dominio de la función a optimizar.
Estos métodos "uno a uno" corren el riesgo de fallar en la localización de óptimos
absolutos, dado que la función no tiene porqué ser convexa. Además, en el caso tratado
en este artículo, supone el manejo de complejos sistemas de ecuaciones no-lineales. Por
su parte, los algoritmos genéticos mantienen en cada iteración un conjunto de potenciales
soluciones que, tras las operaciones genéticas ya citadas, dan lugar a la próxima generación
tras un proceso de selección que solo tiene que ver con el valor de la función
en las mismas.
En el problema de minimización de FO no se han considerado restricciones para
sus variables, salvo la imposición de los intervalos donde se quiere que se muevan dichas
variables.
Los siguientes procedimientos genéticos se han mostrado suficientes: el cruzamiento
aritmético, la mutación no-uniforme y la selección determinista.
El cruzamiento aritmético de dos cromosomas pertenecientes a una generación
{x1,x2} producen dos descendientes {x'1,x'2} de la siguiente forma:
x;=x1+( J-J.,)x2
x;=x2+( J-J.,)x1
donde J., E[O, 1} es elegido aleatoriamente. La probabilidad de cruzamiento se ha fijado
en 0,5 así se espera que, en promedio, el 50% de los cromosomas se crucen.
La mutación no-uniforme del gen Yk del cromosoma genérico x=(y1, ... ,yk, ... ,yJ,
que corresponde con la aproximación a la variable k-esima de la función FO, produce
un descendiente x'=(y¡, ... ,y ·k, ... ,yJ, siendo y'k aleatoriamente seleccionado entre las
siguientes posibilidades:
y~=yk +A(t,y~ -yk)
y~=yk-A(t,yk-yf)
150
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
donde yz, y~ son el extremo inferior y superior del intervalo impuesto donde se mueve
la variable k-esima de FO. La función A(t,y)=y·r·(l-t!T/ devuelve un valor en el intervalo
{O,y], siendo r un número aleatorio en {O, 1], T el número máximo de generaciones,
t el número de la generación actual y b un parámetro que determina el grado de nouniformidad.
Se ha considerado que T=J0.000. La probabilidad de mutación se ha fijado
en O, 1, así se espera que, en promedio, el 10% de genes muten.
La selección determinista es un procedimiento genético por el que se conforma
una nueva generación de cromosomas. Consiste en elegir de entre los cromosomas progenitores
y descendientes de una generación aquellos que lleven aparejados los menores
valores de FO.
0.3,....-----------~-.....
(a)
o 0.2 0.4 x, 0.6 0.8
3 120
~~
'c;
.!i¡ o
,.-.:::
:t
'<>
2
-120
o
--. ------------~-\------------/··í'
-5%
0.2
"'- - -./ /
0.4 X 0.6
1
0.8
r, Figura l(a), Representación de valores experimentales
de la función de Gibbs (O) y de los coeficientes
de actividad(•), obtenidos en el ELV isobárico de
la mezcla binaria (x 1etanoato de propi/o+x2metano/)
y correlaciones obtenidas por los métodos de la
Tabla 1: (- - -), I; (- - ), 11; (- · - ), III; (- ·· - ), IV.
1 (b ), Representación de las diferencias de entalpías,
óh' =h:.,,-h:,,, estimadas con cada uno de los métodos.
0.25-r--------------.. 2.5 100..-----------------.
o 0.2 0.4 x, 0.6 0.8
-100+---..----,...---,..----,,----1
r, o 0.2 0.4 X 0.6
1
0.8
Figura 2(a), Representación de valores experimentales
de la función de Gibbs (O) y de los coeficientes
de actividad (• ), obtenidos en el ELV isobárico de
la mezcla binaria (x 1etanoato de propilo+x2etanol)
y correlaciones obtenidas por los métodos de la
Tabla 1: (- - -), I; (- - ), 11; (- · - ), III; (- ·· - ),IV.
2(b ), Representación de las diferencias de entalpías,
óh '" = h,~, - h.~,, estimadas con cada uno de los métodos.
151
© Del documento, de los autores. Digitalización realizada por ULPGC. Biblioteca Universitaria, 2017
0_2.,..-----------------. 2 120~---------------,_..-,--.
(a) r ' (b)
0.15
1.5
0.05
o 0.2 0.4 x, 0.6 0.8
/ \
5% ---------·-··f \
-:;;.:::::~-- / -------.~- \
_.;...- • J ...,-----r-- ""\.
,-%· -----\ Y------,,:---·, ·-..i
~:-... --· ..... :::..-<~ .... ¡ / ',-~ .. ·'
\ . ;---/ .
\ -~~;,-----./-----------------·
\. ----- -120+---~------~--~----
o 0.2 0.4 x, 0.6 0.8
Figura 3(a), Representación de valores experimentales
de la función de Gibbs (0) y de los coeficientes
de actividad (•), obtenidos en el ELV isobárico de
la mezcla binaria (x1etanoato de propilo+x¡propanol)
y correlaciones obtenidas por los métodos de la
Tabla 1: (- - -), I; (--),U; (- · - ), Ill; (- ·· - ),IV.
3(b ), Representación de las diferencias de entalpías,
óh' = h~, - h.:, , estimadas con cada uno de los mé-todos.
0.17 1-7 100~----------------.
o 0.2 0.4 x, 0.6 0.8
~ ~
5%
r,
-5%.... - ~ "
-100+----.-----.---.---........ ---1
o 0.2 0.4 X 0.6 0.8
1
1_2 Figura 4(a), Representación de valores experimentales
de la función de Gibbs (O) y de los coeficientes
de actividad (•), obtenidos en el ELV isobárico de
la mezcla binaria (x1etanoato de propilo+x2bulanol)
y correlaciones obtenidas por los métodos de la
Tabla 1: (- • -), I; (--),U; (- · - ), 111; (- ·· - ),IV.
4(b ), Representación de las diferencias de entalpías,
óh' =h~1 - h.:,, estimadas con cada uno de los métodos.
152