Estimation of kinetic parameters in a chromatographic separation model via Bayesian inference

Resumo

A modelagem de processos de adsorção tem sido empregada com frequência nas indústrias químicas, petroquímicas e refinarias, por exemplo para separação e purificação de misturas em unidade de Leito Móvel Simulado (LMS). Na representação matemática do modelo, a determinação de parâmetros é um passo importante para o projeto de condições cromatográficas para a separação contínua, em processos do tipo LMS. Este trabalho tem por objetivo a análise de estimativa de parâmetros em processos de adsorção, usando um sistema cromatográfico com uma coluna, para a separação das substâncias Glicose e Frutose. Investiga-se o uso da abordagem Bayesiana, através de métodos de Monte Carlo via Cadeias de Markov (MCMC), assim como o uso da abordagem da máxima verossimilhança, utilizando duas técnicas estocásticas diferentes, o Algoritmo de Colisão de Partículas (PCA - Particle Collision Algorithm), e o Algoritmo de Otimização por Enxame de Partículas (PSO - Particle Swarm Optimization) para executar a tarefa de minimização da função objetivo. Diferentes casos são apresentados com o objetivo de analisar a significância estatística das estimativas obtidas para os parâmetros, fazendo-se uma comparação crítica entre a solução via inferência Bayesiana e via minimização da função objetivo com métodos estocásticos. Os resultados obtidos demonstram que o uso da abordagem Bayesiana fornece uma proposta vantajosa para a estimativa de parâmetros em transferência de massa, oferecendo resultados com maior riqueza de informação estatística.

Palavras-chave: Problemas inversos, Estimativa de parâmetros, Métodos estocásticos, Inferência Bayesiana

1. Introdução

O processo de cromatografia em sistema de coluna tem sido empregado como técnica de purificação para isolar compostos desejados a partir de uma certa mistura. É muito útil em química orgânica, principalmente em problemas de separação de dois compostos, por exemplo na obtenção de fármacos enantioméricos [30]. Os métodos de separação baseiam-se na distribuição de componentes entre as fases móvel e estacionária e as técnicas de separação envolvidas estão relacionadas a fenômenos de adsorção e possuem grande importância na resolução de componentes de misturas racêmicas [25].

No processo de caracterização de coluna, usada na concepção de projetos de unidades para a separação contínua, um fator importante está relacionado à estimação e precisão dos parâmetros obtidos através da solução de problemas inversos. Neste contexto é importante caracterizar a significância estatística das estimativas obtidas devido às incertezas dos dados experimentais empregados [12].

Na estimativa de parâmetros em processos cromatográficos, a abordagem da máxima verossimilhança geralmente é empregada [1], resultando em uma função objetivo a ser minimizada, onde algoritmos de otimização estocástica geralmente são empregados na solução deste problema inverso formulado como um problema de otimização [5,6,7,19], devido às características de evitar estagnação em mínimos locais que estes métodos possuem, em comparação com métodos determinísticos. Por outro lado, estes algoritmos tendem a apresentar alto custo computacional em consequência do grande número necessário de avaliações da função objetivo, uma vez que este é um processo custoso computacionalmente, envolvendo esquemas numéricos para a solução de equações diferenciais. Estratégias híbridas que visam aproveitar as características eficientes dos métodos estocásticos e determinísticos podem ser empregadas para localizar rapidamente a melhor solução e com isso reduzir o custo de processamento [27,32]. Cabe ressaltar que, além de obter estimativas pontuais para os parâmetros de interesse, é de suma importância a avaliação da significância estatística deste resultado, por exemplo por meio do cálculo de regiões de confiança para os valores estimados, uma vez que os efeitos das incertezas de medição presentes nos dados de entrada precisam ser avaliados nos resultados obtidos. Ainda, em problemas inversos não-lineares, como é o caso da grande maioria dos problemas práticos, mesmo que os erros de medição possam ser caracterizados por uma distribuição Normal não necessariamente resulta em estimativas com distribuição de probabilidade desta mesma natureza, como demonstrado em [2]. Neste caso, técnicas mais sofisticadas são necessárias para a estimativa da distribuição de probabilidade dos parâmetros de interesse.

Neste cenário, o emprego de métodos de Monte Carlo via Cadeias de Markov se apresentam como uma boa alternativa na solução do problema inverso [11,18]. O objetivo é amostrar a distribuição de probabilidade que descreve os parâmetros de interesse, de modo que a solução do problema inverso não é uma estimativa pontual, mas uma distribuição de probabilidade que permite avaliar diretamente a região de confiança da estimativa, já levando-se em conta as incertezas dos dados experimentais empregados, assim como informações a priori de parâmetros do modelo, caso estas estejam disponíveis.

Para a formulação e solução do problema inverso, este trabalho apresenta então a abordagem Bayesiana como uma proposta alternativa aos métodos de estimativa de parâmetros em processos cromatográficos. É utilizado o algoritmo de Metropolis-Hastings, que é uma técnica de simulação iterativa baseada em Cadeias de Markov. A ideia deste método é simular amostras da distribuição de interesse, uma vez que esta distribuição não pode ser diretamente obtida, e utilizar estas amostras para inferir a densidade de probabilidades [10]. Desta forma, dentre as características importantes na abordagem Bayesiana com o MCMC estão a possibilidade de avaliação de incertezas e a representação do espaço da solução.

De modo a oferecer uma comparação direta entre a solução via inferência Bayesiana apresentada neste trabalho e o procedimento da Máxima Verossimilhança com utilização de métodos estocásticos para a minimização da função objetivo definida, dois algoritmos de otimização são também implementados, o Algoritmo de Enxame de Partículas (PSO - Particle Swarm Optimization) e o Algoritmo de Colisão de Partículas (PCA - Particle Collision Algorithm). O teste de Kruskal-Wallis é utilizado para identificar a representatividade dos dados e fazer comparações.

2. Modelagem matemática e parâmetros do modelo

Na modelagem dos fenômenos de transferência de massa são usados modelos cinéticos que simulam o processo de adsorção e dessorção através de equações diferenciais simples de primeira ordem, obtidas empregando a metodologia denominada Velocidade da Frente de Convecção (Front Velocity) [6,34]. Para modelar a transferência de massa no sistema cromatográfico, dois modelos de transferência de massa concentrada são assumidos neste trabalho, descritos pelas equações a seguir [5,6,8,19]:

${\frac {dC_{j}}{dt}}=-k_{1j}C_{j}+k_{2j}q_{j}$

(1)

${\frac {dq_{j}}{dt}}=+k_{1j}C_{j}-k_{2j}q_{j}$

(2)

onde ${\textstyle C_{j}}$ e ${\textstyle q_{j}}$ são, respectivamente, as concentrações nas fases líquida e sólida, o índice ${\textstyle j}$ corresponde às espécies das misturas a serem separadas (Glicose e Frutose), ${\textstyle k_{1j}}$ e ${\textstyle k_{2j}}$ representam os parâmetros cinéticos de adsorção e dessorção respectivamente, aqui representados por ${\textstyle {\bf {K}}=\{k_{1},k_{2}\}}$.

A taxa de fluxo que se move para o interior da coluna é determinada por um sistema de bombeamento externo que influencia na velocidade de adsorção/dessorção [6]. O tempo que a fase líquida necessita para escoar ao longo do comprimento da coluna pode ser determinado se a taxa de fluxo volumétrica, a porosidade e o volume da coluna são conhecidos experimentalmente. A velocidade é obtida por meio do fluxo da fase móvel, ou seja, ${\textstyle v=Q/\varepsilon A_{T}}$, onde ${\textstyle Q}$ é a vazão volumétrica, ${\textstyle A_{T}}$ é a área total da coluna, e ${\textstyle \varepsilon }$ é a porosidade. A complexidade do sistema de equações diferenciais abordado no presente trabalho depende da relação entre as concentrações de soluto existentes nas fases líquida e sólida, respectivamente $C(t)$ e $q(t)$. O acoplamento entre as duas equações é feito por intermédio de relações denominadas isotermas que diferem em termos de complexidade dependendo da substância de interesse e da faixa de concentrações do experimento. Na literatura podem ser encontradas relações lineares, bem como isotermas de Langmuir ou bi-camadas de maior complexidade [7].

Na configuração experimental usada neste trabalho, tem-se a vazão volumétrica ${\textstyle Q=30ml/min}$, área total da coluna ${\textstyle A_{T}=2,12cm^{2}}$, volume ${\textstyle V=300ml}$, altura da coluna ${\textstyle h=319cm}$, porosidade ${\textstyle \varepsilon =0,4}$, concentração inicial ${\textstyle C(0)=C_{0}=15mg/ml}$ e diâmetro ${\textstyle D=2,6cm}$ conforme descrito em [19]. Para o cálculo dos perfis individuais das concentrações, as equações diferenciais, equações (1)-(2), foram resolvidas através de um esquema numérico usando o Método de Euler Melhorado.

Na figura 1 estão representados os perfis das concentrações experimentais correspondentes à separação enantiomérica das substâncias Glicose e Frutose usados neste trabalho, extraídos da literatura [3].

Figura 1. Dados experimentais das substâncias Glicose e Frutose obtidos em [3].

3. Problema inverso

A primeira abordagem descrita para a formulação e solução do problema inverso consiste no uso da função de máxima verossimilhança, que resulta em um problema de otimização, onde a função objetivo, que caracteriza os desvios entre dados experimentais e as previsões do modelo deve ser minimizada. A segunda abordagem descrita consiste na utilização da abordagem Bayesiana, que permite a inclusão de informações a priori sobre os parâmetros [26] e tem como objetivo a inferência da distribuição a posteriori dos parâmetros de interesse. Na formulação do problema, com a transferência de massa modelada pelas equações (1-2), considera-se que alguns dados experimentais de concentração de cada substância estão disponíveis, ${\textstyle {\bf {C}}^{exp}}$, como indicado na figura 1, e que ${\textstyle {\bf {K}}}$ é o vetor de parâmetros que se deseja estimar.

3.1 Formulação do problema inverso com o procedimento da Máxima Verossimilhança

Na abordagem clássica, problemas inversos podem ser formulados da seguinte forma: dada uma observação ${\textstyle y}$, deseja-se encontrar a variável ${\textstyle x}$, relacionada com ${\textstyle y}$ por ${\textstyle y=f(x,\epsilon )}$, onde ${\textstyle y\in \Re ^{M}}$ representa a quantidade medida, ${\textstyle x\in \Re ^{n}}$ é a quantidade para a qual se busca obter informações, ${\textstyle \epsilon }$ representa o erro associado às medições experimentais e imprecisões contidas no modelo adotado e ${\textstyle f:\Re ^{n}\times \Re ^{m}\longrightarrow \Re ^{M}}$ é o modelo.

O critério adotado na abordagem consiste em estimar ${\textstyle x}$ de modo a minimizar os desvios entre as predições do modelo e os dados experimentais, ${\textstyle (y-f(x))}$. Para se obter então as estimativas de máxima verossimilhança, o problema inverso é formulado implicitamente, resultando em um problema de otimização, onde deseja-se minimizar o funcional de resíduos quadrados entre as quantidades medidas experimentalmente e os valores calculados pelo modelo. Para o problema proposto tem-se:

$f_{obj}({\bf {K}})=\sum _{i=0}^{N_{c}}\left[C_{i}^{exp}-C_{i}^{calc}({\bf {K}})\right]^{2}$

(3)

onde ${\textstyle C_{i}^{exp}}$ e ${\textstyle C_{i}^{calc}}$ são as concentrações experimentais e calculadas, respectivamente, ${\textstyle {\bf {K}}=(k_{1},...,k_{N_{p}})}$ é o vetor de interesse com ${\textstyle N_{p}}$ parâmetros desconhecidos e ${\textstyle N_{c}}$ é o número de dados experimentais, ou seja, a dimensão do vetor de medidas experimentais, ${\textstyle \mathbf {C} ^{exp}}$. Para o presente trabalho tem-se ${\textstyle N_{p}=2}$ e ${\textstyle N_{c}=38}$ (24 para a Glicose e 14 para a Frutose). A minimização da função objetivo dada pela equação (3) é realizada com os métodos estocásticos PSO e PCA, descritos com mais detalhes a seguir.

3.1.1 Particle Swarm Optimization - PSO

O PSO é um método heurístico de otimização global proposto por Eberhart e Kennedy [13]. É uma técnica de otimização estocástica de base populacional inspirado no comportamento dos pássaros em bando ou cardumes de peixes. No que diz respeito ao algoritmo de otimização, cada indivíduo é chamado de partícula e o processo de busca pela melhor solução se assemelha ao comportamento desses seres, que se movem em bando, ou cardumes, em busca de alimentos. A melhor solução no algoritmo de otimização de partículas pode ser resolvida pela colaboração entre os indivíduos.

No algoritmo por enxame de partículas, os candidatos a solução correspondem a um ponto no espaço de busca. O enxame é iniciado de forma aleatória com uma população de soluções candidatas. Cada partícula ${\textstyle i}$ é inicializada com uma posição ${\textstyle x_{i}(t)}$ e uma velocidade ${\textstyle v_{i}(t)}$ no tempo ${\textstyle t}$, que é modificada levando em conta a melhor posição da partícula (${\textstyle p_{best}}$) e a melhor posição do grupo (${\textstyle g_{best}}$) de acordo com a equação (4)[28]. Dessa forma, ao longo do tempo o grupo percorre o espaço em busca do alimento (solução), de acordo com o seguinte procedimento iterativo:

$v_{i}^{t+1}=w.v_{i}^{t}+c_{1}\times r_{1}(p_{best}-x_{i}^{t})+c_{2}\times r_{2}(g_{best}-x_{i}^{t})$

(4)

Na equação (4), ${\textstyle c_{1}}$ e ${\textstyle c_{2}}$ representam os parâmetros social e cognitivo, respectivamente, que ajustam o balanço entre a influência social e a aprendizagem da partícula individual e ${\textstyle r_{1}}$ e ${\textstyle r_{2}}$ são números aleatórios entre ${\textstyle 0}$ e ${\textstyle 1}$ e ${\textstyle w}$ é o fator de inércia. O Algoritmo 1 apresenta um pseudocódigo simplificado do PSO. Mais detalhes deste algoritmo podem ser encontrados em [14].

Algorithm 1. Algoritmo do Enxame de Partículas - PSO

3.1.2 Particle Collision Algorithm - PCA

O PCA é um método estocástico proposto por Sacco e Oliveira [23], inspirado na interação entre partículas em um reator nuclear, particularmente os fenômenos de absorção e espalhamento. A estrutura deste método se aproxima ao algoritmo de recozimento simulado (Simulated Annealing). Inicialmente, uma configuração inicial (${\textstyle Old\_Conf}$) é selecionada e é definida como melhor solução (${\textstyle Best\_Conf}$). A configuração é modificada por uma perturbação estocástica e leva à construção de uma nova solução (${\textstyle New\_Conf}$). Em seguida essa nova solução é comparada por uma função de aptidão (${\textstyle Fitness}$), dada pela função objetivo, e a nova solução pode ou não ser aceita. Uma função de espalhamento que utiliza o paradigma do tipo Metropolis é usada, onde uma solução candidata pode ser aceita com uma determinada probabilidade, mesmo que leve a uma piora no valor da função objetivo. A exploração em posições mais próximas é garantida através de pequenas perturbações [15]. O pseudocódigo simplificado do PCA é apresentado pelo Algoritmo 2 e a descrição completa pode ser encontrada em [15,22,23,24].

Algorithm 2. Algoritmo de Colisão de Partículas - PCA

3.2 Abordagem Bayesiana em problemas inversos

Na abordagem Bayesiana o problema inverso é reformulado como um problema de inferência e leva em conta as informações pré-existentes e as obtidas através dos dados disponíveis. A abordagem Bayesiana baseia-se na ideia de que os parâmetros do modelo podem ser representados por uma função de probabilidade ${\textstyle p({\bf {K}})}$, chamada de distribuição a priori. Quando fornecida com novas informações, tais como dados experimentais, que têm uma densidade de probabilidade correspondente ${\textstyle p({\bf {C}}^{exp}|{\bf {K}})}$, os parâmetros do modelo podem ser atualizados, resultando numa distribuição a posteriori ${\textstyle p\left({\bf {K}}|{\bf {C}}^{exp}\right)}$ dos parâmetros desconhecidos. O resultado fundamental é expresso pelo teorema de Bayes em problemas inversos, que pode ser expresso por [29]:

$p\left({\bf {K}}|{\bf {C}}^{exp}\right)={\frac {p\left({\bf {C}}^{exp}|{\bf {K}}\right)p\left({\bf {K}}\right)}{p\left({\bf {C}}^{exp}\right)}}$

(5)

onde ${\textstyle p\left({\bf {K}}\right)}$ é a função de densidade a priori, ${\textstyle p\left({\bf {K}}|{\bf {C}}^{exp}\right)}$ é a função de densidade de probabilidade a posteriori, a partir da qual é possível obter informações sobre as incógnitas de ${\textstyle {\bf {K}}}$, e ${\textstyle p\left({\bf {C}}^{exp}|{\bf {K}}\right)}$ é a probabilidade de que a evidência ${\textstyle {\mathbf {\bf {C}} ^{exp}}}$ seja observada para os valores dos parâmetros ${\textstyle {\bf {K}}}$, e é conhecida como função de verossimilhança. O denominador da equação (5) assume, geralmente, o papel de uma constante de normalização e assim adota-se a função de densidade de probabilidade a posteriori como sendo proporcional ao produto da verossimilhança e da distribuição a priori:

$p\left({\bf {K}}|{\bf {C}}^{exp}\right)\propto p\left({\bf {C}}^{exp}|{\bf {K}}\right)p\left({\bf {K}}\right)$

(6)

Na maioria dos casos é razoável considerar o ruído como aditivo branco, com média nula e variância constante e sem correlação entre os erros experimentais, de modo que a matrix de variância é simplificada para $\sigma _{{\bf {C}}^{exp}}^{2}{\bf {I}}$. Neste caso, a função de verossimilhança é dada por:

$p\left({\bf {C}}^{exp}|{\bf {K}}\right)={\left(2\pi \sigma _{{\bf {C}}^{exp}}^{2}\right)^{-{\frac {1}{2}}}}\exp \left\{-{\frac {\left({\bf {C}}^{exp}-C({\bf {K}})\right)^{T}\left({\bf {C}}^{exp}-C({\bf {K}})\right)}{2\sigma _{{\bf {C}}^{exp}}^{2}}}\right\}$

(7)

onde ${\textstyle {\mathbf {C} ^{exp}}}$ refere-se às concentrações experimentais e ${\textstyle {\sigma _{{\bf {C}}^{exp}}^{2}}}$ é a variância relacionada às incertezas das medidas experimentais.

Uma das formas de se determinar a densidade à posteriori, ${\textstyle p\left({\bf {K}}|{{\bf {C}}^{exp}}\right)}$, é através de técnicas de amostragem, como os métodos de Monte Carlo via Cadeias de Markov (MCMC). As cadeias de Markov podem ser geradas a partir de algoritmos específicos, como o Metropolis-Hastings [18]. Neste algoritmo são gerados estados candidatos para a cadeia a partir de uma função densidade de probabilidade de transição ${\textstyle Q({\bf {K}}|{\bf {K'}})}$, que neste trabalho é uma distribuição normal centrada em ${\textstyle {\bf {K}}}$, que fornece a probabilidade de se escolher um novo candidato ${\textstyle {\bf {K}}'}$, supondo que a cadeia esteja em um estado ${\textstyle {\bf {K}}^{(i)}}$, ${\textstyle i=0,1,\dots ,N}$. O novo valor de ${\textstyle {\bf {K}}'}$ é aceito com a seguinte probabilidade:

$\alpha =\min \left[1,{\frac {p\left({\bf {K}}'|{C^{exp}}\right)Q\left({\bf {K}}',{\bf {K}}^{(i)}\right)}{p\left({\bf {K}}^{(i)}|{C^{exp}}\right)Q\left({\bf {K}}^{(i)},{\bf {K}}'\right)}}\right]$

(8)

onde ${\textstyle p}$ é a distribuição a posteriori de interesse, calculada de acordo com a equação (5).

Uma vez que ${\textstyle Q}$ é uma distribuição simétrica, ou seja, ${\textstyle Q({\bf {K}}|{\bf {K}}')=Q({\bf {K}}'|{\bf {K}})}$ para todos os valores de ${\textstyle {\bf {K}}}$ e ${\textstyle {\bf {K}}'}$, a probabilidade de aceitação se reduz à razão entre as densidades calculadas e independem de ${\textstyle Q}$, conforme apresentado a seguir:

$\alpha =\min \left[1,{\frac {p\left({\bf {K}}'|{\bf {C}}^{exp}\right)}{p\left({\bf {K}}^{(i)}|{\bf {C}}^{exp}\right)}}\right]$

(9)

Comparam-se então as probabilidades à posteriori do candidato com o estado atual, ${\textstyle p\left({\bf {K}}'|{{\bf {C}}^{exp}}\right)}$ e ${\textstyle p\left({\bf {K}}|{C^{exp}}\right)}$. Se esta probabilidade for maior para o candidato, ${\textstyle {\bf {K}}'}$, aceita-se automaticamente, ${\textstyle {\bf {K}}^{(i+1)}={\bf {K}}'}$. No entanto, mesmo que ${\textstyle p\left({\bf {K}}'|{\bf {C}}^{exp}\right)}$ seja menor que ${\textstyle p\left({\bf {K}}|{\bf {C}}^{exp}\right)}$, tal que ${\textstyle \alpha {<1}}$, o valor candidato pode ainda ser aceito. Isso envolve uma amostragem aleatória de uma densidade uniforme ${\textstyle U^{(i)}}$ e o valor candidato é aceito sujeito a ${\textstyle \alpha \geq U^{(i)}}$. O Algoritmo 3 então é esquematizado de acordo com os passos descritos em [10,16].

O resultado fornecido pelo algoritmo é a construção de uma cadeia de Markov com os estados ${\textstyle ({\bf {K}}^{0},{\bf {K}}^{1},\dots ,{\bf {K}}^{N_{MCMC}})}$. Estes valores são amostras da distribuição a posteriori e podem, portanto ser usados na inferência de informações desta distribuição de interesse [18].

Algorithm 3. Algoritmo de Metropolis-Hasting

Na distribuição a priori deve-se representar o conhecimento que se tem sobre ${\textstyle {\bf {K}}}$ antes da realização do experimento. Neste trabalho considera-se que a única informação a priori disponível sobre ${\textstyle \mathbf {K} }$ é o domínio a qual pertencem os parâmetros. Portanto, a informação a priori é modelada como uma distribuição uniforme entre os valores mínimos e máximos estabelecidos para cada parâmetro.

Devido à necessidade de constatar a convergência da cadeia resultante do método de Monte Carlo, utilizou-se o Critério de Raftery-Lewis na avaliação de convergência além das análises gráficas obtidas. O critério de Raftery-Lewis permite calcular o número de iterações necessárias para uma cadeia atingir a distribuição estacionária, ou seja, permite estimar o número de iterações iniciais que devem ser descartadas e da distância mínima de uma iteração à outra para se obter uma amostra independente [20].

4. Análise de sensibilidade

Um fator importante na estimação dos parâmetros com a aplicação do problema inverso está relacionado à análise de sensibilidade de parâmetros [12,33]. Neste caso, deseja-se estudar o comportamento da resposta do modelo para variações nos parâmetros cinéticos de adsorção e dessorção. O objetivo é observar os parâmetros de transferência de massa contidos na modelagem do sistema cromatográfico e ao fazer a análise de sensibilidade para cada um deles, deve-se verificar quais os parâmetros possuem maior influência sobre o modelo.

Visando comparar os valores dos coeficientes de sensibilidade para parâmetros com diferentes unidades dimensionais, utilizam-se os coeficientes de sensibilidade modificados. Estes são obtidos através do produto do coeficiente de sensibilidade pelo valor do respectivo parâmetro. Dessa forma obtém-se todos os coeficientes com as mesmas unidades. Para a análise de sensibilidade dos coeficientes, considerou-se então a seguinte equação:

${\bf {X}}_{k_{i}}=k_{i}\times {\frac {\partial C_{j}({\bf {K}})}{\partial k_{i}}},\quad i=1,2,\dots ,N_{p}$

(10)

onde ${\textstyle k_{i}}$ representa as incógnitas do problema inverso e o termo diferencial é resolvido usando diferenças finitas centradas conforme a equação (11), considerando ${\textstyle \varepsilon _{i}=0,01\times k_{i}}$:

${\frac {\partial C_{j}({\bf {K}})}{\partial k_{i}}}\approx {\frac {C_{j}\left(k_{1},\dots ,k_{i}+\varepsilon _{i},\dots ,k_{N}\right)-C_{j}\left(k_{1},\dots ,k_{i}-\varepsilon _{i},\dots ,k_{N}\right)}{2\varepsilon _{k_{i}}}}$

(11)

A obtenção de estimativa de boa qualidade para um problema inverso está relacionada à sensibilidade da resposta do modelo aos parâmetros de interesse. Neste sentido, é desejável que os coeficiente de sensibilidade dados pela equação (10) possuam alta magnitude. Além disso, se dois ou mais parâmetros estiverem sendo estimados simultaneamente no problema inverso, seus coeficiente de sensibilidade não devem apresentar dependência linear entre si, caso contrário não há unicidade na solução do problema inverso [4,33].

5. Resultados e discussões

5.1 Análise de sensibilidade e validação dos métodos estocásticos de otimização

Os gráficos apresentados na figura 2 mostram o comportamento de sensibilidade do modelo aos parâmetros que se deseja determinar, onde observa-se que variações nos parâmetros causam essencialmente efeitos sobre a forma da curva das componentes. Ou seja, as curvas de sensibilidade para os parâmetros cinéticos possuem variações distintas em relação à variável observável do problema. Nos dois casos, verifica-se que os parâmetros das componentes apresentam alguma correlação, mas não são linearmente dependentes, indicando que os parâmetros ainda podem ser simultaneamente estimados.

Figura 2. Resultado para a análise de sensibilidade com respeito aos parâmetros das componentes. (a) Glicose. (b) Frutose, respectivamente.

Na tabela 1 estão apresentadas a média (${\textstyle \mu }$), desvio padrão (${\textstyle \sigma }$) e a dispersão (${\textstyle d}$) dos parâmetros de transferência de massa da Glicose e Frutose, obtidas com o procedimento da máxima verossimilhança. Foram realizadas 15 execuções dos algoritmos estocásticos PSO e PCA estabelecendo-se 2700 avaliações da função objetivo em cada execução. Estas simulações utilizam o mesmo conjunto de dados experimentais ilustrados na figura 1. Os valores de referência obtidos em [3] também estão apresentados nesta tabela.

Tabela 1. Representação da média, desvio padrão, coeficiente de variação das estimativas e a dispersão dos parâmetros das substâncias Glicose e Frutose, usando o conjunto de dados experimentais da figura (1).

Alg.	Subst.	$k_{1}$($ml/mg.min$)			$k_{2}$($ml/mg.min$)
		$\mu$	$\sigma$	$\sigma /\mu \times 100\%$	$\mu$	$\sigma$	$\sigma /\mu \times 100\%$	$d$
PSO	Glic.	$0,014300$	$1,8970\times {10}^{-6}$	$0,013$	$0,023581$	$3,016\times {10}^{-6}$	$0,012$	$0,00013$
	Frut.	$0,014934$	$7,931\times {10}^{-4}$	$5,41$	$0,011299$	$1,349\times {10}^{-3}$	$5,55$	$0,05487$
PCA	Glic.	$0,014658$	$1,1\times {10}^{-4}$	$0,73$	$0,024265$	$8,161\times {10}^{-9}$	$0,00$	$0,00769$
	Frut.	$0,014619$	$6,4576\times {10}^{-4}$	$4,41$	$0,010981$	$4,981\times {10}^{-4}$	$4,53$	$0,04476$
Ref.	Glic.	$k_{1}=0,01511$ $ml/mg.min$		$k_{2}=0,02810$ $ml/mg.min$
	Frut.	$k_{1}=0,01454$ $ml/mg.min$		$k_{2}=0,01332$ $ml/mg.min$

Cabe ressaltar que o desvio padrão apresentado refere-se apenas à estocasticidade dos métodos e não considera incertezas provenientes de ruído experimental. A dispersão em relação à média das estimativas apresentada na tabela 1 foi calculada empregando a equação (12) [9], apresentada abaixo, e mostra a discrepância entre os métodos estocásticos. Pode-se observar que a maior variação dos parâmetros encontra-se nas execuções do PCA.

$d={\sqrt {\frac {(1+(\sigma _{k_{2}}/\mu _{k_{2}}))^{2}+(1+(\sigma _{k_{1}}/\mu _{k_{1}}))^{2}}{2}}}-1$

(12)

5.2 Resultado da aplicação da Inferência Bayesiana

Neste trabalho foram empregados os dados experimentais reais apresentados na figura 1, obtidos em [5]. Ressalta-se que o desvio padrão dos erros de medição, ${\sigma _{C^{exp}}}$, não foram reportados, de modo que neste trabalho considerou-se ${\sigma _{{\bf {C}}^{exp}}}=0,1mg/ml$, o que equivale a cerca de 0,8% do pico de concentração de glicose, e cerca de 1,25% do pico de concentração de frututose, sendo um ruído experimental factível.

Os gráficos das figuras subfig_A e 3a mostram a evolução da cadeia de Markov com os parâmetros da Glicose e Frutose respectivamente. A cadeia, para ambas as substâncias, parece ter convergido rapidamente e com isso foram descartadas apenas 5% das iterações iniciais (1250 iterações) para cada parâmetro, considerando que estas amostras iniciais são influenciadas pelos estados iniciais das cadeias. Este percentual é superior ao constatado pelo Critério de Raftery-Lewis, que indica que sejam descartados apenas 2,1% dos estados para a Glicose e 1,9% para Frutose.

(a) Glicose

(b) Frutose

Figura 3. Evolução das cadeias de Markov dos parâmetros de transferência de massa. (a) Glicose. (b) Frutose juntamente a média a posteriori e o valor de referência [5,7].

O número mínimo estimado de iterações necessárias para se obter a convergência, ainda segundo o critério de Raftery-Lewis, é de ${\textstyle N_{MCMC}=3746}$ para ambas as substâncias e um máximo de ${\textstyle N_{MCMC}=12842}$ para Glicose e ${\textstyle N_{MCMC}=13076}$ para Frutose. Isso mostra que para o problema proposto não é necessário um elevado número de iterações para a convergência da cadeia e para a inferência da densidade à posteriori. Na representação gráfica, figuras subfig_A e 3a, pode-se observar que não houve discrepância durante o processo de geração de amostras ou tendências de instabilidades no processo iterativo.

Com base nas probabilidades amostrais e descartando os estados iniciais da cadeia, foram obtidos os histogramas de densidade a posteriori, a normal ajustada e a densidade de convergência correspondente a cada parâmetro das substâncias. Na figura 4 estão apresentados os histogramas juntamente com o valor de referência [5,7] dos parâmetros e a média a posteriori.

(a) Glicose	(b) Glicose
(c) Frutose	(d) Frutose
Figura 4. Histograma das observações simuladas, normal ajustada, gráfico da densidade de convergência para (a) e (b) Glicose e (c) e (d) Frutose e o valor de referência [5,7].

Em todos os casos as distribuições a posteriori apresentaram um comportamento simétrico e os gráficos da função de densidade de probabilidade estimada tendem para uma distribuição normal. Ao mesmo tempo verifica-se uma pequena diferença das estimativas obtidas para o parâmetro de dessorção com relação ao valor de referência [5,7]. Este resultado está associado tanto ao modelo de coluna empregado quanto ao conjunto de dados experimentais utilizados. É importante ressaltar que o valor de referência claramente está dentro do intervalo de confiança de ${\textstyle 95\%}$ da distribuição a posteriori estimada.

Os gráficos de dispersão na figura 5 apresentam a região de confiança de percentil 95 dos parâmetros obtidos. As regiões foram calculadas usando o método Minimum Volume Ellipsoid [21] utilizado para filtrar dados que desviam de forma significativa do padrão normal do conjunto de parâmetros (outliers) e se baseia em um elipsoide que cobre grande parte das observações com menor dispersão.

Verifica-se no gráfico a tendência de correlação entre ${\textstyle k_{2}}$ e ${\textstyle k_{1}}$ devido a indicação de dependência dessas condicionais, ou seja, elas têm as mesmas variações, mas ${\textstyle k_{2}}$ varia com ${\textstyle k_{1}}$ no sentido em que ambos tendem a aumentar em conjunto. Este comportamento pode ser observado nos gráficos de sensibilidade na figura 2 se, ao aumentar ou diminuir a variação dos parâmetros, ambos tendem a compensar o efeito do outro em sentido oposto.

(a) Glicose	(b) Frutose
Figura 5. Gráfico de dispersão com a região de intensidade das estimativas obtidas para as substâncias. (a) Glicose. (b) Frutose.

5.3 Intervalos de Confiança

Finalmente, após o uso das abordagens foi possível calcular os intervalos de 95% de confiabilidade para cada parâmetro do modelo. Para a distribuição a posteriori calculada com o MCMC, os intervalos foram obtidos desprezando 2,5% das amostras iniciais e finais, estando estas ordenadas.

Em relação aos métodos estocásticos (PSO e PCA), para estimar os intervalos de confiança de 95% foram realizadas 15 execuções dos algoritmos com diferentes conjuntos de dados experimentais. O propósito de simular conjuntos distintos de dados consiste em tentar capturar a incerteza tanto devido à estocasticidade do método quanto devido ao ruído experimental. Para isso, ruídos foram simulados nos dados experimentais existentes usando o seguinte esquema:

${\hat {C}}^{exp}=C^{exp}+\epsilon _{C^{exp}},\epsilon _{C^{exp}}\sim N(0,\sigma _{C^{exp}})$

(13)

Dessa forma, diferentes conjuntos de amostras aleatórias são obtidos e em seguida os algoritmos estocásticos são executados para cada conjunto, simulando a realização de experimentos independentes em cada execução. Essa simulação difere no entanto, do caso apresentado na avaliação isolada do efeito da estocasticidade do método, tabela 1. Os desvios padrão contidos nesta simulação possuem duas componentes, uma devido à natureza estocástica dos métodos e outra devido às flutuações nos dados experimentais.

Dois esquemas foram usados na obtenção do intervalo de confiança para os algoritmos estocásticos (PSO e PCA). Inicialmente considerou-se o uso da média e o desvio padrão das amostras (${\textstyle \mu \pm {1},96\sigma }$), resultante das 15 execuções com os diferentes conjuntos de dados experimentais. No outro esquema, o desvio padrão para as estimativas é calculado através da seguinte equação [17]

$\sigma _{K_{i}}=\sigma _{C^{exp}}{\sqrt {\left[J^{T}J\right]_{ii}^{-1}}},i=1,2,\dots ,N$

(14)

onde ${\textstyle J_{ji}=\partial C_{j}/\partial K_{i}}$ são os elementos da matriz de sensibilidade calculados pela equação (11):

${\frac {\partial C(K)}{\partial K}}=\left|{\begin{array}{cccc}{\dfrac {\partial C_{1}}{\partial k_{1}}}&{\dfrac {\partial C_{1}}{\partial k_{2}}}&\cdots &{\dfrac {\partial C_{1}}{\partial k_{N}}}\\{\dfrac {\partial C_{2}}{\partial k_{1}}}&{\dfrac {\partial C_{2}}{\partial k_{2}}}&\cdots &{\dfrac {\partial C_{2}}{\partial k_{N}}}\\\vdots &\vdots &\ddots &\vdots \\{\dfrac {\partial C_{M}}{\partial k_{1}}}&{\dfrac {\partial C_{M}}{\partial k_{2}}}&\cdots &{\dfrac {\partial C_{M}}{\partial k_{N}}}\end{array}}\right|$

(15)

e a confiabilidade de 95% do intervalo para cada parâmetro é obtido por meio da seguinte expressão:

$K_{i}-1,96\sigma _{K_{i}}\leq K_{i}\leq K_{i}+1,96\sigma _{K_{i}}i=1,\dots ,N$

(16)

A média e o desvio padrão obtidos com os algoritmos PSO e PCA usando os 15 conjuntos de dados distintos estão apresentados na tabela 2. Além disso, são apresentadas também as médias e os desvios calculados após a remoção de outliers dos conjuntos de parâmetros obtidos com o algoritmo Metropolis-Hastings. Verifica-se que a presença destes outliers não influenciaram significativamente na variação da média. Isto confirma o que foi mencionado na seção 5.2 sobre a baixa discrepância no processo de geração de amostras.

Em cada execução dos algoritmos PCA e PSO foram consideradas 2700 avaliações da função objetivo e como já mencionado, os desvios padrão calculados para as constantes cinéticas neste caso, possuem duas componentes e refletem, portanto, a aleatoriedade dos resultados dos métodos estocásticos e o efeito da variação dos dados experimentais usados nas estimativas.

A tabela 2 também apresenta os coeficientes de variação, ${\textstyle \sigma /\mu \times 100\%}$. Para os resultados obtidos com os algoritmos PSO e PCA, os desvios relativamente à média apresentam menor variação comparado ao MCMC, onde os desvios são mais expressivos. Isso se justifica pelo baixo número de execuções dos algoritmos estocásticos de otimização, implicando, portanto, menor confiabilidade no valor do desvio padrão calculado.

Tabela 2. Média e desvio padrão para $k_{1}$ e $k_{2}$ obtidas com 2700 avaliações da função objetivo nos intervalos $k_{1}=(0,01;0,05)$ e $k_{2}=(0,001;0,06)$ com diferentes conjuntos de dados experimentais.

Alg.	Subst.	$k_{1}$($ml/mg.min$)			$k_{2}$($ml/mg.min$)
		$\mu$	$\sigma$	$\sigma /\mu \times 100\%$	$\mu$	$\sigma$	$\sigma /\mu \times 100\%$
PSO	Glicose	$0,01346$	$0,00119$	$8,85$	$0,02452$	$0,00204$	$8,35$
	Frutose	$0,01454$	$0,00164$	$11,27$	$0,01171$	$0,00123$	$10,5$
PCA	Glicose	$0,01391$	$0,00115$	$8,33$	$0,02540$	$0,00212$	$8,36$
	Frutose	$0,01452$	$0,00183$	$12,6$	$0,01167$	$0,00146$	$12,5$
MCMC	Glicose	$0,01498$	$0,00145$	$9,67$	$0,02583$	$0,00256$	$9,91$
	Frutose	$0,01445$	$0,00214$	$14,8$	$0,01162$	$0,00180$	$15,49$
MCMC$^{1}$	Glicose	$0,01492$	$0,00134$	$9,01$	$0,02570$	$0,00233$	$9,10$
	Frutose	$0,01427$	$0,00187$	$13,11$	$0,01147$	$0,00157$	$13,73$

${\textstyle ^{1}}$ Médias e desvios calculados após a remoção de outliers.

As estimativas dos parâmetros, apresentadas na tabela 2, foram usadas na solução do problema direto, assim como os valores de referência destes parâmetros [6,5]. Os valores calculados para a concentração são apresentados nos gráficos da figura 7.

(a) Glicose

(b) Frutose

Figura 6. Representação do perfil cromatográfico. (a) Glicose e (b) Frutose, calculados através do problema direto com os parâmetros apresentados na tabela (2).

Verifica-se nas duas curvas, a similaridade entre os dados experimentais e as concentrações calculadas utilizando os parâmetros estimados. Uma maior discrepância é observada na figura 7(b), correspondente à Frutose. Este fato, muito provavelmente, pode ser explicado pela abordagem simplificada adotada neste trabalho para descrever a cinética da coluna, na qual são desconsiderados o coeficiente de dispersão axial e a resistência à transferência de massa.

Para identificar a representatividade dos dados da tabela 2 utilizou-se o teste de Kruskal-Wallis [31], com o objetivo de verificar a existência de diferenças significativas entre conjuntos de parâmetros dos algoritmos, sendo o resultado apresentado na tabela 3, onde ${\textstyle H}$ denota a estatística de Kruskal-Wallis e ${\textstyle {\mathcal {X}}^{2}}$ é o valor crítico da distribuição qui-quadrado.

Tabela 3. Comparação das amostras dos algoritmos PSO, PCA e MCMC usando o teste de Kruskal-Wallis

Substância	$[ml/mg.min]$	$H<{\mathcal {X}}^{2}$	p-valor
Glicose	$k_{1}$	$1,3881<5,9914$	$0,4995$
	$k_{2}$	$1,7456<5,9914$	$0,4177$
Frutose	$k_{1}$	$0,0964<5,9914$	$0,9529$
	$k_{2}$	$0,2031<5,9914$	$0,9034$

De acordo com as regras do teste, verifica-se na tabela 3 que não existe diferença significativa das estimativas entre os algoritmos, indicando portanto que as amostras provêm de populações com a mesma mediana. Além do valor estatístico analisado, ${\textstyle H<{\mathcal {X}}^{2}}$, em todos os casos obteve-se p-valor maior que a significância ${\textstyle \alpha =0,05}$ considerada.

Na tabela 4 estão representados os intervalos de confiança de 95% para os parâmetros ${\textstyle k_{1}}$ e ${\textstyle k_{2}}$ calculados da distribuição a posteriori com o MCMC e os intervalos calculados a partir das execuções dos algoritmos PSO e PCA com os conjuntos de dados simulados com a equação (13). Os intervalos sinalizados com (*) foram calculados com a equação (14) e (#) representa o intervalo de confiança para o total de iterações definido pelo critério de Raftery-Lewis. Também são apresentados os intervalos de confiança para o conjunto de dados após o filtro de outliers, onde se verifica que estes não influenciaram na representação das amostras.

Tabela 4. Intervalos de confiança de $95$% para os parâmetros $k_{1}$ e $k_{2}$ calculados com a média e desvio padrão das 15 simulações estocásticas e com a abordagem Bayesiana.

Alg.	Glicose		Frutose
	$k_{1}[ml/mg.min]$	$k_{2}[ml/mg.min]$	$k_{1}[ml/mg.min]$	$k_{2}[ml/mg.min]$
PSO	$(0,01112;0,01580)$	$(0,02050;0,02853)$	$(0,01133;0,01776)$	$(0,00928;0,01414)$
	$*(0,01059;0,01634)$	$*(0,02077;0,02826)$	$*(0,01198;0,01710)$	$*(0,00862;0,01480)$
PCA	$(0,01164;0,01618)$	$(0,02124;0,02956)$	$(0,01093;0,01811)$	$(0,00880;0,01454)$
	$*(0,01103;0,01678)$	$*(0,02165;0,02915)$	$*(0,01196;0,01708)$	$*(0,00858;0,01476)$
MCMC	$(0,01228;0,01755)$	$(0,02101;0,03028)$	$(0,01060;0,01794)$	$(0,00838;0,01456)$
	${}^{1}(0,01233;0,01798)$	${}^{1}{0},0211;0,03095)$	${}^{1}(0,01011;0,01850)$	${}^{1}(0,00813;0,01508)$
	$\#(0,01122;0,01567)$	$\#(0,02095;0,02933)$	$\#(0,01088;0,01972)$	$\#(0,00834;0,01572)$

${\textstyle {}^{1}}$ intervalos calculados sem a remoção de outliers.

Como no problema tratado não foi considerado o uso de informação a priori informativa, pode-se observar na tabela 4 que, com as execuções realizadas com o PSO e PCA usando diferentes amostras, foi possível obter intervalos similares aos obtidos com a inferência Bayesiana. Por outro lado, o uso da inferência fornece mais informações estatísticas a respeito dos parâmetros.

Quando se aumenta o conjunto de amostras para inferir o desvio padrão para as amostras obtidas via PSO e PCA, destaca-se o aumento do custo computacional dos algoritmos estocásticos frente à simulação com a abordagem Bayesiana, visto que, na execução do MCMC com 25000 avaliações da função objetivo, obtém-se um conjunto de amostras das estimativas e isso é equivalente a uma execução do método estocástico com igual número de avaliações da função, porém com a obtenção de apenas uma solução. Ou seja, o custo computacional para se obter um conjunto de amostras confiável usando o PCA ou PSO seria muito elevado. Para a obtenção dos dados apresentados na tabela 2 por exemplo, cada algoritmo (PSO e PCA) realizou ${\textstyle 15\times {2700}}$ avaliações da função objetivo.

Ao atentar-se às iterações sugeridas pelo critério de convergência usado no MCMC, nota-se que os intervalos de confiança estão, numericamente, de acordo com aqueles encontrados com os métodos estocásticos e ainda, com um menor número de avaliação da função objetivo.

6. Conclusões

Neste trabalho foi resolvido um problema inverso visando estimar parâmetros de transferência de massa em um modelo de cromatografia empregando a abordagem Bayesiana. Nesta metodologia, foi utilizado um conjunto de dados experimentais reais e foi realizada uma comparação crítica com os resultados obtidos usando o procedimento da Máxima Verossimilhança e métodos estocásticos de otimização. Observou-se que os intervalos de confiança não apresentaram variação significativa de amplitude ao adotar as especificações do critério de convergência de Raftery-Lewis nas cadeias de Markov, havendo assim uma redução expressiva do custo computacional. Com base nas análises estatísticas estudadas e considerando que os algoritmos estocásticos utilizados necessitam de um grande número de execuções e de avaliação da função objetivo para se obter o mesmo um conjunto de amostras, a abordagem Bayesiana se apresentou como uma alternativa vantajosa para a estimação de parâmetros em problema de transferência de massa. Além disso, o uso desta estratégia fornece mais confiabilidade na determinação dos parâmetros do modelo devido à representatividade estatística das amostras.

Agradecimentos

Os autores agradecem o apoio financeiro da FAPERJ, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, da CAPES, Fundação Coordenação de Aperfeiçoamento de Pessoal de Nível Superior e do CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico.

Referências

[1] M. Schwaab e J. C. Pinto, Análise de dados experimentais I - Fundamentos de estatística e estimação de parâmetros, E-papers, Rio de Janeiro, 2007.

[2] M. Schwaab, E. C. Biscaia Jr, J. L. Monteiro and J. C. Pinto, Nonlinear parameter estimation through particle swarm optimization, Chemical Engineering Science, v.63, n.6, (2008) 1542-1552.

[3] D. C. S. Azevedo e A. Rodrigues. SMB chromatography applied to the separation/purification of fructose from cashew apple juice. Brazilian Journal of Chemical Engineering [online]. ISSN: 0104-6632, 17, 4-7, (2000)507-516.

[4] J.V. Beck, B. Blackwell, and C.R. St. Clair Jr, Inverse Heat Conduction - Ill-Posed Problems, John Wiley & Sons, New York, 1985

[5] A. L. J. Bihain, L. D. T. Câmara, A. J. Silva Neto. Avaliação da rotina inversa R2W na estimação de parâmetros de transferência de massa no processo de adsorção de Glicose e Frutose. TEMA, São Carlos, v.13, n.3, (2012).

[6] A. L. J. Bihain, “Desenvolvimento e avaliação de novas abordagens de modelagem de processos de separação em Leito Móvel Simulado”, Tese de Doutorado, IPRJ, UERJ, Nova Friburgo, RJ, 2014.

[7] L. D. T. Câmara, J. Lugon Junior, A. J. Silva Neto. Modelagem Computacional de Fenômenos de Separação por Adsorção. ISBN: 978-85-63749-25-3, Ed. Alternativa, Niterói, 2014

[8] L. D. T. Câmara, Stepwise Model Evaluation in Simulated Moving-Bed Separation of Ketamine, Chemical Engineering & Technology, v.37, n.2, (2014) 1521-4125.

[9] A. P. C. Cuco, A. J. Silva Neto, H. F. Campos Velho & F. L. de Sousa. Solution of an inverse adsorption problem with an epidemic genetic algorithm and the generalized extremal optimization algorithm. Inverse Problems in Science and Engineering, 17 (2009)289-302.

[10] S. Jackman. “Estimation and Inference via Bayesian Simulation: An Introduction to Markov Chain Monte Carlo.” American Journal of Political Science. 44 (2000)375-404.

[11] J. Kaipio, E. Somersalo. Statistical and Computational Inverse Problems, New York: Springer, 2005, pp.339.

[12] J. Kalyanaraman, et al. Bayesian estimation of parametric uncertainties, quantification and reduction using optimal design of experiments for ${\textstyle CO_{2}}$ adsorption on amine sorbents. Computers and Chemical Engineering (2015)1-13.

[13] J. Kennedy, R. Eberhart. Particle Swarm Optimization: Proc of IEEE International Conference on Neural Network, Perth, Austrália, IEEE Service Center Piscataway NJ. (1995)1942-1948.

[14] J. Kennedy. The particle swarm: Social Adaptation of knowledge, IEEE International Conference on Evolutionary Computation. (1997)303-308.

[15] E. F. P. Luz, J. C. Becceneri, H. F. de Campos Velho. A new multi-particle collision algorithm for optimization in a high performance environment. Journal of Computational Interdisciplinary Sciences, DOI: 10.6062/jcis.2008.01.01.0001, 1 (2008)3-10.

[16] E. Medova. Bayesian Analysis and Markov Chain Monte Carlo Simulation. Wiley StatsRef: Statistics Reference Online. Nova Iorque: John Wiley e Sons. (2007)1-20.

[17] M. N. ${\textstyle {\ddot {O}}}$zisik e H. R. B. Orlande, Inverse Heat Transfer: Fundamentals and Applications, Taylor & Francis, New York, 2000

[18] M. Strand. Metropolis-Hastings Markov Chain Monte Carlo, Chapman University. (2009)1-9.

[19] A. Prieto-moreno, O. L. Santiago, L. D. T. Câmara, A. J. Silva Neto e C. Oliveira. Uncertainty analysis in mass transfer parameter estimations for chromatographic separation of glucose and fructose. Anais do Congresso Nacional de Matemática Aplicada e Computacional - CNMAC, 2014

[20] A. L. Raftery, S. Lewis. How Many Iterations in the Gibbs Sampler? In Bayesian Statistics 4 (J.M. Bernardo et al., editors). Oxford University Press, (1992)763-773.

[21] P. J. Rousseeuw & B. C. V. Zomeren. Unmasking multivariate outliers and leverage points. Journal of the American Statistical Association, 85 (411), (1990)633-639.

[22] W. F. Sacco , C. R. E. de Oliveira. A new stochastic optimization algorithm based on particle collisions. Proceedings of the 2005 ANS Annual Meeting. Transactions of the American Nuclear Society. 92 (2005)657-659.

[23] W. F. Sacco, C. R. E. de Oliveira. A new stochastic optimization algorithm based on particle collisions metaheuristic. In ${\textstyle 6^{th}}$ World Congresses of Structural and Multidisciplinary Optimization, Rio de Janeiro, Brazil, (2005)1-6.

[24] W. F. Sacco, C. M. F. Lapa, C. M.  N. A. Pereira, H. Alves Filho. A Metropolis Algorithm applied to a Nuclear Power Plant Auxiliary Feedwater System surveillance tests policy optimization. ISSN: 0149-1970, Progress in Nuclear Energy. 50 (2008)15-21.

[25] M. A. G. Santos, V. Veredas, I. J. Silva Jr, C. R. D. Correia, L. T. Furlan, C. C. Santana. Simulated Moving-Bed Adsorption for Separation of Racemic Mixtures. Brazilian Journal of Chemical Engineering. 21 (2004)127-136.

[26] L. G. Silva, D. C. Knupp, L. Bevilacqua, A. C. N. R. Galeão, A. J. Silva Neto. Formulação e Solução de um Problema Inverso de difusão Anômala com Técnicas Estocásticas, Ciência e Natura. 36 (2014)82-96.

[27] A. J. Silva Neto e J. C. Becceneri. Técnicas de Inteligência Computacional Inspiradas na Natureza - Aplicação em Problemas Inversos em Transferência Radiativa, 2ª Edição. Sociedade Brasileira de Matemática Aplicada e Computacional - SBMAC (Notas em Matemática Aplicada, v. 41), São Carlos (2012).

[28] G. J. Simões, N. F. F. Ebecken, Algoritmo genético e enxame de partículas para a otimização de suportes laterais de fornos, Rev. int. métodos numér cálc.diseño ing. (2015)1-6.${\textstyle <}$http://dx.doi.org/10.1016/j.rimni.2014.07.001${\textstyle >}$

[29] A. M. Stuart. Inverse problems: a Bayesian perspective. Acta Numerica 19 (2010)451-559.

[30] A. K. Singh, E. R. M. Kedor-Hackmann, M. I. R. M. Santoro. Cromatografia líquida com fase quiral aplicada na separação enantiomérica de fármacos cardiovasculares. Brazilian Journal of Pharmaceutical Sciences. 42 (2006)553-566.

[31] N. Thomson. Understanding ANOVA the apl way. Em APL ’93 International Conference on Array-Programming Languages, New York, NY, USA, (1993)295-303.

[32] J. Xu et. al, Determination of competitive adsorption isotherm of enantiomers on preparative chromatographic columns using inverse method Journal of Chromatography A, 1273 (2013) 49-56.

[33] J. F. V. Vasconcellos, A. J. Silva Neto & C. C. Santana. An Inverse Mass Transfer Problem in Solid-Liquid Adsorption Systems: Inverse Problems in Engineering. 11 (2003)391-408.

[34] M. I. Mesa, L. D. T. Câmara, D. C. Knupp & A. J. Silva Neto. Uncertainty Quantification in Chromatography Process Identification Based on Markov Chain Monte Carlo. In: A. J. Silva Neto, O. Llanes-Santiago, G. N. Silva. Mathematical Modeling and Computational Intelligence in Engineering Applications (2016) pp. 77-88. Springer International Publishing.