Metodos Estocasticos de Otimizacaono Calculo de Azeotropos Reativos

Gustavo M. Platt,Depto de Modelagem Computacional, IPRJ, UERJ,

28630-050, Nova Friburgo, RJ

E-mail: gmplatt@iprj.uerj.br

Resumo: Neste trabalho e apresentada uma analise da utilizacao de metodos estocasticos deotimizacao na solucao de sistemas algebricos de equacoes nao-lineares com aplicacoes em Ter-modinamica: o calculo de azeotropos reativos. A existencia (ou nao) de azeotropos em misturasimplica em diferentes estrategias de separacao dos componentes da mistura; portanto, o calculorobusto de azeotropos e importante do ponto de vista industrial. O problema de azeotropia reativae apresentado no sistema de coordenadas composicionais de Ung-Doherty. Sao apresentados re-sultados numericos para calculos de azeotropos reativos em um sistema quaternario reativo ideale com duas estrategias de otimizacao estocastica: o metodo de Luus-Jaakola e o algoritmo deEvolucao Diferencial. Os resultados indicam que – ao menos neste problema – o metodo deEvolucao Diferencial oferece uma alternativa de simples implementacao e extremamente efici-ente e robusta do ponto de vista computacional.

Palavras-chave: Otimizacao, Evolucao Diferencial, Algoritmo de Luus-Jaakola, Calculo deAzeotropos.

1 Introducao

Um azeotropo (em seu conceito tradicional) e uma condicao termodinamica onde uma mistura(em fase lıquida) produz um vapor com composicao identica (ou seja, as quantidades relativasdos componentes nas fases lıquida e vapor sao iguais). Estados azeotropicos podem introduzirregioes de confinamento de produtos para separacoes que exploram diferencas nos pontos deebulicao dos componentes da mistura. Diversos trabalhos abordam, tanto do ponto de vistaexperimental quanto computacional, o problema de azeotropia (ver, por exemplo, [5, 10, 12]).

Em alguns casos, azeotropos podem ocorrem em sistemas reativos sob condicoes de equilıbrioquımico. O conceito de azeotropo reativo (em sistemas com uma unica reacao quımica emequilıbrio) foi introduzido por Barbosa e Doherty [1]. Neste caso, as composicoes (fracoes mo-lares) nao sao iguais na condicao azeotropica, mas sao relacionadas pela equacao:

y1 − x1ν1 − νTx1

=yi − xi

νi − νTxii = 2, . . . , C − 1 (1)

onde yi e xi referem-se as fracoes molares do componente i nas fases vapor e lıquida, respecti-vamente; νi e o coeficiente estequiometrico de i na reacao quımica; νT e a soma dos coeficientesestequiometricos na reacao em equilıbrio e C e o numero de componentes.

Empregando uma transformacao nao-linear nos vetores de fracoes molares, Ung e Doherty[14] produziram um novo conjunto de coordenadas composicionais transformadas para as faseslıquida e vapor, Xi e Yi, respectivamente, que sao identicos na condicao azeotropica:

Xi =

(xi − νt

iN−1xref

1− νttotN

−1xref

);Yi =

(yi − νt

iN−1yref

1− νttotN

−1yref

)(2)

118

ISSN 1984-8218

onde νti refere-se ao vetor de coeficientes estequiometricos do componente i nas reacoes quımicas

das quais ele participa (t indica transposicao), νttot e um vetor com as somas dos coeficien-

tes estequiometricos para cada reacao, N e uma matriz de coeficientes estequiometricos para oscomponentes de referencia (o numero de componentes de referencia e igual ao numero de reacoesquımicas em equilıbrio), xref e o vetor de fracoes molares (na fase lıquida) para os componen-tes de referencia. yref e definido de maneira analoga, para a fase vapor. Os componentes dereferencia sao escolhidos dentre as especies participantes da reacao quımica, evitando-se compo-nentes inertes, que introduziriam linhas nulas na matriz N e impediriam a inversao da matriz.

Azeotropos homogeneos reativos podem ser observados mesmo em misturas ideais, enquantoazeotropos homogeneos ordinarios sao consequencias de nao-idealidades da fase lıquida. Alemdisto, diversos sistemas de interesse apresentam azeotropos reativos e, em alguns casos, mais deum azeotropo reativo [9].

O problema de calculo de azeotropos reativos e caracterizado por um sistema de equacoesalgebricas nao-lineares, representando as condicoes de coexistencia entre as fases fluidas, ascondicoes de equilıbrio quımico e as relacoes de azeotropia reativa. Este sistema algebricopode ser convertido em um problema de minimizacao e resolvido por metodologias classicas(determinısticas) ou por metodos estocasticos de otimizacao.

O calculo de azeotropos reativos foi previamente estudado por Maier et al. [8] utilizandouma metodologia do tipo Newton intervalar/bisecao generalizada. Este algoritmo tem garantiasmatematicas e computacionais, de acordo com Maier et al. [8], de robustez para a localizacaode todos os azeotropos reativos; por outro lado, e extremamente custoso do ponto de vistacomputacional. Bonilla-Petriciolet et al. [3] apresentaram uma nova formulacao para o calculode azeotropos reativos e nao-reativos e obtiveram resultados com uma versao do algoritmoSimulated Annealing (SA) (desenvolvido por Goffe et al. [6] e implementado por Corana et al.[4]) e um hıbrido do tipo Simulated Annealing-Newton. Os resultados obtidos demonstram queo algoritmo hıbrido apresenta caracterısticas de robustez; contudo, deve-se ter em conta que aversao do algoritmo SA utilizada nao e a versao canonica (a implementacao de Corana et al. [4]incorpora uma busca unidimensional no metodo) e, portanto, ha uma lacuna no que diz respeitoa capacidade de metaheurısticas classicas em resolver o problema de azeotropia reativa. Reis etal. [11] apresentaram um algoritmo robusto para a identificacao de azeotropos reativos, baseadona construcao de uma malha de estimativas iniciais e no uso de um metodo determinıstico paraa solucao do sistema algebrico nao-linear representado pelas relacoes de equilıbrios de fases equımicos. Quando o algoritmo identifica a vizinhanca de um azeotropo (||X − Y || proxima azero), um metodo de otimizacao e iniciado de modo a localizar o azeotropo reativo.

Neste trabalho e abordado o calculo de azeotropia reativa empregando-se dois metodos deotimizacao estocasticos: o algoritmo de Luus-Jaakola [7] e o metodo de Evolucao Diferencial[13]. Os algoritmos serao comparados no que diz respeito a robustez, tempo de computacao,numero de avaliacoes da funcao-objetivo e facilidade de implementacao.

2 Formulacao do Problema de Azeotropia Reativa

Nesta secao serao apresentados o calculo de azeotropia reativa e o problema de otimizacaoassociado.

2.1 Equilıbrio Simultaneo de Fases e Quımico e o Azeotropo Reativo

Um azeotropo reativo homogeneo e uma condicao de coexistencia de fases fluidas (lıquida e va-por) onde as coordenadas composicionais, no sentido Ung-Doherty [14], sao identicas. Portanto,tal condicao pode ser representada como:

Xi = Yi i = 1, 2, . . . , C −R− 1 (3)

119

ISSN 1984-8218

onde C e o numero de componentes na mistura e R representa o numero de reacoes quımicasem equilıbrio.

O equilıbrio simultaneo de fases e quımico e expresso pelo seguinte sistema algebrico nao-linear (para sistemas a baixas pressoes e sem associacoes em fase vapor):

Pyi − xiγiPsati = 0 i = 1, 2, . . . , C∑C

i=1 xi − 1 = 0∑Ci=1 yi − 1 = 0

ln (K)− νt ln (a) = 0 i = 1, 2, . . . , R

(4)

O primeiro conjunto de equacoes representa a condicao de coexistencia lıquido-vapor (Lei deRaoult modificada). As duas equacoes seguintes sao as restricoes de soma unitaria das fracoesmolares. O ultimo conjunto de equacoes contempla as condicoes de equilıbrio quımico. Nasequacoes (4), P e a pressao total do sistema, γi e o coeficiente de atividade para o componentei e P sat

i e a pressao de saturacao para i (tambem conhecida como pressao de vapor). Ainda nasEqs. (4), K e o vetor de constantes de equilıbrio quımico. Em sistemas monoreativos, ν e ovetor com os coeficientes estequiometricos de cada componente participante da reacao quımicae, naturalmente, K e um escalar. O vetor a e o vetor de atividades (em fase lıquida), sendo aatividade de um componente i calculada por ai = γixi.

Uma analise de graus de liberdade para este sistema algebrico nao-linear indica C + 2 + Requacoes. As coordenadas termodinamicas sao x, y, T e P , em um total de 2C + 2 variaveis.Portanto, o numero de graus de liberdade para o problema de coexistencia de fases e (2C+2)−(C + R + 2) = C − R. Ora, ao considerarmos a pressao especificada (ou a temperatura, semperda de generalidade), o numero de graus de liberdade cai para C − R − 1. Este numero deequacoes pode ser obtido utilizando-se as relacoes de azeotropia reativa (Eqs. 3), de modo aresolver-se o sistema algebrico nao-linear.

As pressoes de vapor sao calculadas via Equacao de Antoine (mm Hg):

log10(P sati

)= Ai −

Bi

T (oC) + Ci(5)

2.2 O problema de otimizacao

O sistema algebrico nao-linear que representa as condicoes simultaneas de equilıbrios de fases equımicos (Eqs. 4) e as condicoes de azeotropia reativa (Eqs. 3) podem ser convertidos em umproblema de otimizacao, a ser resolvido via metodos estocasticos.

O vetor de incognitas pode ser escrito como: θt = [x1 x2 . . . xc y1 y2 . . . yc T ], sob pressaoespecificada. Portanto, o sistema nao-linear escreve-se como:

Fl (θ) = 0 l = 1, 2, . . . , 2C + 1 (6)

onde Fl refere-se a cada equacao nao-linear formada pelas equacoes (4) e (3).Portanto, um problema de otimizacao relacionado pode ser construıdo, atraves da seguinte

funcao escalar nao-negativa f(θ):

f (θ) =

2C+1∑l=1

(Fl (θ))2 (7)

Um minimizador θ da funcao f(θ) que produz um mınimo nulo e tambem uma solucao dosistema nao-linear representado pelas Eqs. (3-4).

120

ISSN 1984-8218

3 Metodos Numericos

3.1 O algoritmo de busca aleatoria adaptativa (algoritmo de Luus-Jaakola)

A ideia basica do algoritmo de Luus-Jaakola [7] e simples: uma busca alaetoria e feita nodomınio da funcao e, periodicamente, o raio da regiao de busca e diminuıdo, de modo a produziruma solucao acurada. O metodo, assim como outras metaheurısticas, tem alguns parametros decontrole: uma estimativa inicial aleatoria (θ0), o numero de iteracoes no laco interno (jmax,antesde reduzir a regiao de busca), o numero de iteracoes do laco externo (imax), o raio de buscar (um vetor com as mesmas dimensoes de θ), o fator de contracao ϵ (usado para diminuir, aofinal de um conjunto de iteracoes internas, o raio de busca). A partir de um otimo corrente,um novo ponto e gerado por uma expressao do tipo: θnovo = θ0 + Rr, onde R e uma matrizdiagonal com elementos nao-nulos aleatorios com distribuicao uniforme no intervalo [−0, 5; 0, 5].Caso o valor da funcao no ponto θnovo seja menor do que em θ0, faz-se θ0 ← θnovo. Osexperimentos numericos neste trabalho foram feitos com imax = 400, jmax = 200, ϵ = 0, 98,rt = [1 1 1 1 1 1 100]. O raio de busca e atualizado por uma expressao do tipo: r ← (1− ϵ)r.

3.2 O algoritmo de Evolucao Diferencial

O algoritmo de Evolucao Diferencial foi proposto por Storn e Price [13] e baseia-se em umaestrutura populacional que evolui atraves de 3 mecanismos: mutacao, recombinacao e selecao.Para cada ponto da populacao (vetor), o vetor mutante e calculado por uma expressao do tipovi,G+1 = θr1,G+F (θr2,G−θr3G), onde os pontos r1, r2 e r3 sao distintos entre si e do ponto i. Oındice G refere-se a geracao. O parametro F e um numero no intervalo [0, 2]. Apos a mutacao,o vetor-tentativa e obtido comparando-se, elemento a elemento, os vetores θi,G e vi,G+1, em umprocesso conhecido como recombinacao. Um elemento uji do vetor-tentativa sera vji,G+1 casorandb(j) < CR ou j = rnbr(i). Por outro lado, uji sera θji,G caso randb(j) < CR e j = rnbr(i).CR e um numero entre 0 e 1, randb(j) e um numero aleatorio com distribuicao uniforme 0-1 ernbr(i) e um ındice aleatorio dentre os elementos do vetor. Finalmente, na selecao, o valor dafuncao-objetivo entre os pontos ui e θi,G e comparado e o ponto com menor valor e escolhido paraa proxima geracao (G+1). O ponto base para a mutacao pode ser aleatorio (representado comorand) ou o melhor dos pontos (best), dentre outras alternativas. A mutacao pode ser aplicadaa partir de um unico par de pontos ou a diversos pares de pontos. Neste trabalho, escolhemos aopcao DE/rand/1, indicando que o ponto escolhido e aleatorio e que a mutacao ocorre sobre umunico par de pontos. Neste trabalho, consideramos N = 75 (tamanho da populacao), G = 350,CR = 0, 8803, F = 0, 4717.

4 Resultados

Conforme mencionado anteriormente, azeotropos reativos podem aparecer mesmo em sistemasideais (fase vapor gas ideal e fase lıquida solucao ideal). Neste trabalho, consideraremos umsistema quaternario ideal (proposto por Barbosa e Doherty [2]) sob pressao de 1 atm, onde umareacao quımica em equilıbrio ocorre: A + B⇔ C + D. Portanto, os coeficientes estequiometricossao νA = −1; νB = −1; νC = 1; νD = 1, e, portanto, νT = 0 ou, equivalentemente, νt

tot = 0. Umavez que a fase lıquida e considerada solucao ideal, consideraremos γi = 1, para i ∈ {A,B,C,D}.Os valores para as constantes da equacao de Antoine sao apresentados na Tabela 1.

Considerando D como componente de referencia, as condicoes de azeotropia reativa nestesistemas sao expressas como XA = YA e XB = YB. Portanto, em termos de fracoes molares,aplicando-se a definicao da transformacao Ung-Doherty [14], tem-se: xA + xD = yA + yD exB + xD = yB + yD. Estas equacoes, juntamente com as condicoes de coexistencia de fasese as relacoes de equilıbrio quımico (Eqs. 4) formam o sistema de equacoes algebricas a seremresolvidas para a determinacao do azeotropo reativo.

121

ISSN 1984-8218

Componente Ai Bi Ci

A 7,38781 1533,313 222,299B 8,07130 1730,630 233,416C 8,11219 1592,864 226,174D 7,10178 1244,951 217,871

Tabela 1: Constantes para Equacao de Antoine

A Tabela 2 apresenta as coordenadas composicionais e temperatura para o azeotropo reativono sistema quaternario em questao. Estes valores foram obtidos com estrategias do tipo Newton,com estimativas iniciais proximas a solucao do problema, e servem apenas para comparacoescom os resultados obtidos pelos algoritmos estocasticos.

Componente xi yi T (K)

A 0,1883651 0,0732845 362,54111B 0,3583621 0,2432814 -C 0,2150891 0,3301697 -D 0,2381837 0,3532644 -

Tabela 2: Azeotropo reativo em sistema ideal sob pressao de 1 atm.

A Figura 1 representa uma execucao tıpica do algoritmo de Luus-Jaakola na busca peloazeotropo reativo quaternario no sistema ideal. O comportamento e apresentado em coordenadas(XA, XB, T ), mas o problema esta, de fato, em ℜ7, uma vez que tem-se C − 1 coordenadaspara a fase lıquida, C − 1 coordenadas para a fase vapor e a temperatura. A evolucao doalgoritmo e apresentada simultaneamente com o diagrama de temperaturas de bolha reativas(superfıcie na figura). O azeotropo reativo, neste caso, corresponde ao ponto de sela, claramenteindicado na figura e encontrado pelo metodo LJ. Valores das coordenadas de Ung-Dohertyacima de 1 ou negativos, embora nao facam parte do domınio composicional Ung-Doherty,podem ser produzidas a partir dos valores de fracoes molares durante o processo de busca.Um aspecto interessante nesta Figura diz respeito aos valores de temperaturas ao longo dasiteracoes: os valores tornam-se extremamente baixos e, finalmente, retornam para as vizinhancasda temperatura azeotropica.

Figura 1: Execucao tıpica do algoritmo de Luus-Jaakola.

122

ISSN 1984-8218

Analogamente, a Figura 2 ilustra os pontos iniciais e finais para a abordagem com DE/rand/1.Assim como no caso do LJ, o metodo DE/rand/1 localizou a condicao de azeotropia reativa paraquaisquer estimativas iniciais. Os pontos finais estao colapsados, por questoes de escala, em umunico triangulo, exatamente sobre o ponto de sela.

Figura 2: Execucao tıpica do DE/rand/1. ×: pontos iniciais. △ :pontos finais.

A Tabela 3 apresenta uma comparacao (para 10 execucoes de cada algoritmo) entre osmetodos Luus-Jaakola e Evolucao Diferencial (DE/rand/1) no calculo de azeotropos reativos.Os parametros de cada algoritmo foram ajustados de modo a obter-se, a cada execucao, valoresde funcao-objetivo proximos a zero. De fato, nao ha – para os dois algoritmos – valores dafuncao-objetivo acima de 1 × 10−7. Os dois algoritmos foram implementados com o mesmonıvel de otimizacao no que diz respeito a sua programacao e executados no mesmo equipamento.Portanto, as diferencas nao podem ser explicadas por qualidades de implementacao distintas oupor maquinas distintas. Claramente, o algoritmo de Evolucao Diferencial – neste problema – esuperior, tanto do ponto de vista de qualidade da solucao quanto ao que diz respeito ao tempode computacao, ao metodo proposto por Luus-Jaakola. Todos os calculos foram obtidos emlinguagem Scilab 5.3.3 (64 bits) [15], em um computador com processador Intel I7 com 4GB deRAM.

Luus-Jaakola Evolucao Diferencial

Funcao-objetivo (media) 2, 772× 10−8 2, 301× 10−14

Funcao-objetivo (melhor) 1, 130× 10−8 4, 107× 10−17

Tempo (medio, s) 167,15 44,38Tempo (melhor, s) 144,23 40,98

Tabela 3: Comparacao entre os metodos LJ e DE (10 instancias para cada algoritmo).

5 Conclusoes

Neste trabalho foram obtidos resultados de calculo de azeotropos reativos em um sistema qua-ternario ideal, empregando metodos de otimizacao estocasticos. O problema de azeotropia rea-tivo foi apresentado, bem como os parametros necessarios ao calculo do azeotropo. O sistemaalgebrico nao-linear que representa o calculo de azeotropos em sistemas reativos foi obtido econvertido em um problema de minimizacao de uma funcao nao-negativa. A estrutura de oti-mizacao foi entao resolvida pelo metodos de Luus-Jaakola [7] e Evolucao Diferencial [13]. Os

123

ISSN 1984-8218

resultados indicaram uma clara superioridade, em termos da qualidade do otimo e do tempo decomputacao, do algoritmo de Evolucao Diferencial frente ao metodo de Luus-Jaakola.

Agradecimentos

Gustavo M. Platt agradece pelo apoio financeiro da Fundacao Carlos Chagas Filho de Apoio aPesquisa do Estado do Rio de Janeiro (FAPERJ).

Referencias

[1] Barbosa, D. and Doherty, M.F. (1987) ‘Theory of phase diagrams and azeotropic conditionsfor two-phase reactive systems’, Proc. R. Soc. London, Vol. A413, pp.443–458.

[2] Barbosa, D. and Doherty, M.F. (1988) ‘The influence of equilibrium chemical reactions onvapor-liquid phase diagrams’, Chemical Engineering Science, Vol. 43, pp.529–540.

[3] Bonilla-Petriciolet, A., Iglesias-Silva, G.A. and Hall, K.R. (2009) ‘Calculation of homo-geneous azeotropes in reactive and non-reactive mixtures using a stochastic optimizationapproach’, Fluid Phase Equilibria, Vol. 281, pp.22–31.

[4] Corana, A., Marchesi, M., Martini, C., Ridella, S. (1987) ‘Minimizing multimodal functionsof continuous variables with the simulated annealing algorithm’, ACM Transactions onMathematical Software, Vol. 13, pp.262–280.

[5] Chapman, R.G. and Goodwin, S.P. (1993) ‘A general algorithm for the calculation of aze-otropes in fluid mixtures’, Fluid Phase Equilibria, Vol. 85, pp.55–69.

[6] Goffe, W.L., Ferrier, G.D. and Rogers, J. (1994) ‘Global optimization of statistical functionswith simulated annealing’, Journal of Econometrics, Vol. 60, pp.65–99.

[7] Luus, R. and Jaakola, T.H.I. (1973) ‘Optimization by direct search and systematic reductionof the size of search region’, AIChE Journal, Vol.19, pp.760.

[8] Maier, R.W., Brennecke, J.F. and Stadtherr, M.A. (2000) ‘Reliable computation of reactiveazeotropes’, Computers and Chemical Engineering, Vol. 24, pp.1851–1858.

[9] Okasinski, M. J. and Doherty, M. F. (1997) ‘Thermodynamic behavior of reactive azeotro-pes’, AIChE Journal, Vol. 43, pp.2227–2238.

[10] Qi, Z. and Sundmacher, K. (2005) ‘Geometrically Locating Azeotropes in Ternary Systems’,Ind. Eng. Chem. Res., Vol. 44(10), pp.3709–3719.

[11] Reis, M.H.M., Mascolo, L.F.S and Wolf-Maciel, M.R. (2006) ‘Development of a robustalgorithm to compute reactive azeotropes’, Brazilian Journal of Chemical Engineering, Vol.23(3), pp.395-403.

[12] Song, W., Huss, R.S., Doherty, M.F. and Malone, M.F. (1997) ‘Discovery of a reactiveazeotrope’, Nature, Vol. 388, pp.561–563.

[13] R. Storn, K. Price, Differential evolution – a simple and efficient heuristic for global opti-mization over continuous spaces, J. Global Optimization, 11 (1997) 341-359.

[14] Ung, S. and Doherty, M.F. (1995) ‘Necessary and sufficient conditions for reactive azeotro-pes in multireaction mixtures’, AIChE Journal, Vol. 41, pp.2383–2392.

[15] Consortium Scilab - Digiteo (2011). Scilab: Free and open source software for numericalcomputation (Windows, Version 5.3.3) [Software]. Disponıvel em: http://www.scilab.org.

124

ISSN 1984-8218