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 <span id='citeF-30'></span>[[#cite-30|[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 <span id='citeF-25'></span>[[#cite-25|[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 <span id='citeF-12'></span>[[#cite-12|[12]]].

Na estimativa de parâmetros em processos cromatográficos, a abordagem da máxima verossimilhança geralmente é empregada <span id='citeF-1'></span>[[#cite-1|[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 <span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-19'></span>[[#cite-5|[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 <span id='citeF-27'></span><span id='citeF-32'></span>[[#cite-27|[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 <span id='citeF-2'></span>[[#cite-2|[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 <span id='citeF-11'></span><span id='citeF-18'></span>[[#cite-11|[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 <span id='citeF-10'></span>[[#cite-10|[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'') <span id='citeF-6'></span><span id='citeF-34'></span>[[#cite-6|[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 <span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-19'></span>[[#cite-5|[5,6,8,19]]]:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dC_j}{dt}=-k_{1j}C_j+k_{2j}q_j   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dq_j}{dt}=+k_{1j}C_j-k_{2j}q_j   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

onde <math display="inline">C_j</math> e <math display="inline">q_j</math> são, respectivamente, as concentrações nas fases líquida e sólida, o índice <math display="inline">j</math> corresponde às espécies das misturas a serem separadas (Glicose e Frutose), <math display="inline">k_{1j}</math> e <math display="inline">k_{2j}</math> representam os parâmetros cinéticos de adsorção e dessorção respectivamente, aqui representados por <math display="inline">{\bf K} = \{k_1, k_2\}</math>.

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 <span id='citeF-6'></span>[[#cite-6|[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, <math display="inline">v = Q/\varepsilon A_T</math>, onde <math display="inline">Q</math> é a vazão volumétrica, <math display="inline">A_T</math> é a área total da coluna,  e <math display="inline">\varepsilon </math> é 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 <math>C(t)</math> e <math>q(t)</math>. 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 <span id='citeF-7'></span>[[#cite-7|[7]]].

Na configuração experimental usada neste trabalho, tem-se a vazão volumétrica <math display="inline">Q = 30 ml/min</math>, área total da coluna <math display="inline">A_T = 2,12 cm^2</math>, volume <math display="inline">V = 300ml</math>,  altura da coluna <math display="inline">h=319cm</math>, porosidade <math display="inline">\varepsilon = 0,4</math>,  concentração inicial <math display="inline">C(0)=C_0=15mg/ml</math> e diâmetro <math display="inline">D=2,6cm</math> conforme descrito em <span id='citeF-19'></span>[[#cite-19|[19]]]. Para o cálculo dos perfis individuais das concentrações, as equações diferenciais, equações ([[#eq-1|1]])-([[#eq-2|2]]), foram resolvidas através de um esquema numérico usando o Método de Euler Melhorado.

Na figura [[#img-1|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 <span id='citeF-3'></span>[[#cite-3|[3]]].  <div id='img-1a'></div>
<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-Glucose_Frutose.png|500px|Dados experimentais das substâncias Glicose e Frutose obtidos em <span id='citeF-3'></span>[[#cite-3|[3]]].]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;" | '''Figura 1.''' Dados experimentais das substâncias Glicose e Frutose obtidos em <span id='citeF-3'></span>[[#cite-3|[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 <span id='citeF-26'></span>[[#cite-26|[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 ([[#eq-1|1]]-[[#eq-2|2]]), considera-se que alguns dados experimentais de concentração de cada substância estão disponíveis, <math display="inline">{\bf C}^{exp}</math>, como indicado na figura [[#img-1|1]], e que <math display="inline">{\bf K}</math> é 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 <math display="inline">y</math>, deseja-se encontrar a variável <math display="inline">x</math>, relacionada com <math display="inline">y</math> por <math display="inline">y = f(x, \epsilon )</math>, onde <math display="inline">y\in \Re ^M</math> representa a quantidade medida, <math display="inline">x\in \Re ^n</math> é a quantidade para a qual se busca obter informações, <math display="inline">\epsilon </math> representa o erro associado às medições experimentais e imprecisões contidas no modelo adotado e <math display="inline">f:\Re ^n\times \Re ^m\longrightarrow \Re ^M</math> é o modelo.

O critério adotado na abordagem consiste em estimar <math display="inline">x</math> de modo a minimizar os desvios entre as predições do modelo e os dados experimentais, <math display="inline">(y-f(x))</math>.  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:

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_{obj}({\bf K}) = \sum _{i=0}^{N_c} \left[C^{exp}_{i} - C^{calc}_{i}({\bf K})\right]^2  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

onde <math display="inline">C^{exp}_{i}</math> e <math display="inline">C^{calc}_{i}</math> são as concentrações experimentais e calculadas, respectivamente, <math display="inline">{\bf K} = (k_1,..., k_{N_p})</math> é o vetor de interesse com <math display="inline">N_p</math> parâmetros desconhecidos e <math display="inline">N_c</math> é o número de dados experimentais, ou seja, a dimensão do vetor de medidas experimentais, <math display="inline">\mathbf{C}^{exp}</math>. Para o presente trabalho tem-se <math display="inline">N_p = 2</math> e <math display="inline">N_c = 38</math> (24 para a Glicose e 14 para a Frutose). A minimização da função objetivo dada pela equação ([[#eq-3|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 <span id='citeF-13'></span>[[#cite-13|[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 <math display="inline">i</math> é inicializada com uma posição <math display="inline">x_i(t)</math> e uma velocidade <math display="inline">v_i(t)</math> no tempo <math display="inline">t</math>, que é modificada levando em conta a melhor posição da partícula (<math display="inline">p_{best}</math>) e a melhor posição do grupo (<math display="inline">g_{best}</math>) de acordo com a equação ([[#eq-4|4]])<span id='citeF-28'></span>[[#cite-28|[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:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

Na equação ([[#eq-4|4]]), <math display="inline">c_1</math> e <math display="inline">c_2</math> 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 <math display="inline">r_1</math> e <math display="inline">r_2</math> são números aleatórios entre <math display="inline">0</math> e <math display="inline">1</math> e <math display="inline">w</math> é o fator de inércia. O Algoritmo [[#algorithm-1|1]] apresenta um pseudocódigo simplificado do PSO. Mais detalhes deste algoritmo podem ser encontrados em <span id='citeF-14'></span>[[#cite-14|[14]]].


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|style="padding:10px;"| 
[[File:Draft_Samper_522746929_5375_Algoritmo1.png]]
|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-1'></span>'''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 <span id='citeF-23'></span>[[#cite-23|[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 (<math display="inline">Old\_Conf</math>) é selecionada e é definida como melhor solução (<math display="inline">Best\_Conf</math>). A configuração é modificada por uma perturbação estocástica e leva à construção de uma nova solução (<math display="inline">New\_Conf</math>). Em seguida essa nova solução é comparada por uma função de aptidão (<math display="inline">Fitness</math>), 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 <span id='citeF-15'></span>[[#cite-15|[15]]]. O pseudocódigo simplificado do PCA é apresentado pelo Algoritmo [[#algorithm-2|2]] e a descrição completa pode ser encontrada em <span id='citeF-15'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span>[[#cite-15|[15,22,23,24]]].


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|style="padding:10px;"|
[[File:Draft_Samper_522746929_5126_Algoritmo2.png]]
|-style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"|'''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 <math display="inline">p({\bf K})</math>, chamada de distribuição a priori. Quando fornecida com novas informações, tais como dados experimentais, que têm uma densidade de probabilidade correspondente <math display="inline">p({\bf C}^{exp}|{\bf K})</math>, os parâmetros do modelo podem ser atualizados, resultando numa distribuição a posteriori <math display="inline">p\left({\bf K}|{\bf C}^{exp}\right)</math>  dos parâmetros desconhecidos. O resultado fundamental é expresso pelo teorema de Bayes em problemas inversos, que pode ser expresso  por <span id='citeF-29'></span>[[#cite-29|[29]]]:

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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)}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

onde <math display="inline">p\left({\bf K}\right)</math> é a função de densidade a priori, <math display="inline">p\left({\bf K}|{\bf C}^{exp}\right)</math> é a função de densidade de probabilidade a posteriori, a partir da qual é possível obter informações sobre as incógnitas de <math display="inline">{\bf K}</math>, e <math display="inline">p\left({\bf C}^{exp}|{\bf K}\right)</math> é a probabilidade de que a  evidência <math display="inline">{\mathbf{\bf C}^{exp}}</math> seja observada para os valores dos parâmetros <math display="inline">{\bf K}</math>, e é conhecida como ''função de verossimilhança''. O denominador da  equação ([[#eq-5|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:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p\left({\bf K}|{\bf C}^{exp}\right)\propto p\left({\bf C}^{exp}|{\bf K}\right)p\left({\bf K}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <math>\sigma _{{\bf C}^{exp}}^2 {\bf I}</math>. Neste caso, a função de verossimilhança é dada por:

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p\left({\bf C}^{exp}|{\bf K}\right)= {\left(2\pi \sigma ^2_{{\bf C}^{exp}}\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 ^2_{{\bf C}^{exp}}}}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

onde <math display="inline">{\mathbf{C}^{exp}}</math> refere-se às concentrações experimentais e <math display="inline">{\sigma _{{\bf C}^{exp}}^2}</math> é a variância relacionada às incertezas das medidas experimentais.

Uma das formas de se determinar a densidade à posteriori, <math display="inline">p\left({K}|{C^{exp}}\right)</math>, é 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 <span id='citeF-18'></span>[[#cite-18|[18]]]. Neste algoritmo são gerados estados candidatos para a cadeia a partir de uma função densidade de probabilidade de transição <math display="inline">Q({K}|{K'})</math>, que neste trabalho é uma distribuição normal centrada em <math display="inline">{K}</math>, que fornece a probabilidade de se escolher um novo candidato <math display="inline">{K'}</math>, supondo que a cadeia esteja em um estado <math display="inline">{K}^{(i)}</math>, <math display="inline">i = 0, 1, \dots , N</math>. O novo valor de <math display="inline">{K'}</math> é aceito com a seguinte probabilidade:

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha = \min \left[1, \frac{p\left({K'}|{C^{exp}}\right)Q\left({K'}, {K}^{(i)}\right)}{p\left({K}^{(i)}|{C^{exp}}\right)Q\left({K}^{(i)},{K'}\right)}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

onde <math display="inline">p</math> é a distribuição a posteriori de interesse, calculada de acordo com a equação ([[#eq-5|5]]).

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha = \min \left[1, \frac{p\left({K'}|{C^{exp}}\right)}{p\left({K}^{(i)}|{C^{exp}}\right)}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

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

O resultado fornecido pelo algoritmo é a construção de uma cadeia de Markov com os estados <math display="inline">({K}^0, {K}^1,\dots ,{K}^{N_{MCMC}})</math>. 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 <span id='citeF-18'></span>[[#cite-18|[18]]].


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:InputDefinir  OutputSaída <math>N</math>, <math>{K}^{(0)}</math> N é Tamanho da cadeia de Markov <math>{K}^N</math>  '''for''' <math>i\gets </math> <math>N</math> '''do'''
::
:'''end'''
: <math>{K'} \gets  </math> amostra da distribuição <math>Q({K'}|{K}^{(i)})</math>
|-
| <math>u \gets  U_{[0,1]}</math>Dist. uniforme entre [0,1] <math>\alpha ({K}, {K'}) \gets </math> equação ([[#eq-8|8]])
|-
| <math>\alpha \geq u</math> <math>K^{(i+1)} \gets  {K'}</math>
: <math>{K}^{(i+1)} \gets  {K}^{(i)}</math>
: <math display="inline">{K}^N</math>;  


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-3'></span>'''Algorithm. 3''' Algoritmo de Metropolis-Hasting
|}

Na distribuição a priori deve-se representar o conhecimento que se tem sobre <math display="inline">{K}</math> antes da realização do experimento. Neste trabalho considera-se que a única informação a priori disponível sobre <math display="inline">\mathbf{K}</math> é 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 <span id='citeF-20'></span>[[#cite-20|[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 <span id='citeF-12'></span><span id='citeF-33'></span>[[#cite-12|[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:

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{X}_{k_i} = k_i\times \frac{\partial C_j({K})}{\partial k_i}, \quad i = 1, 2, \dots , N_p    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

onde <math display="inline">k_i</math> representa as incógnitas do problema inverso e o termo diferencial é resolvido usando diferenças finitas centradas conforme a equação ([[#eq-11|11]]), considerando <math display="inline">\varepsilon _i = 0,01 \times k_i</math>:

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{\partial C_j({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}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 ([[#eq-10|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  <span id='citeF-4'></span><span id='citeF-33'></span>[[#cite-4|[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 [[#img-2|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.    <div id='img-2a'></div>
<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-SensibilidadeR_Glucose.png|450px|]]
|[[Image:draft_Samper_522746929-SensibilidadeR_Frutose.png|400px|Resultado para a análise de sensibilidade com respeito aos parâmetros das componentes (a) Glicose e (b) Frutose respectivamente.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Resultado para a análise de sensibilidade com respeito aos parâmetros das componentes (a) Glicose e (b) Frutose respectivamente.
|}
Na tabela [[#table-1|1]] estão apresentadas a média (<math display="inline">\mu </math>), desvio padrão (<math display="inline">\sigma </math>) e a dispersão (<math display="inline">d</math>) 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 [[#img-1|1]]. Os valores de referência obtidos em <span id='citeF-3'></span>[[#cite-3|[3]]] também estão apresentados nesta tabela.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 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 ([[#img-1|1]]).
|- style="border-top: 2px solid;"
| style="text-align: left;" |   
|  
| colspan='2' | <math>k_{1}</math>(<math>ml/mg.min</math>)
| 
|-
| colspan='2' | <math>k_{2}</math>(<math>ml/mg.min</math>)
| 
| <math>\mu </math>
| <math>\sigma </math>
| <math>\sigma /\mu \times 100%</math> 
| <math>\mu </math>
| <math>\sigma </math>
| <math>\sigma /\mu \times 100%</math>
|- style="border-top: 2px solid;"
| style="text-align: left;" | <math>d</math>
| Glic. 
| <math>0,014300</math>
| <math>1,8970\times{10}^{-6}</math>
| <math>0,013</math>
| <math>0,023581</math>
| <math>3,016\times{10}^{-6}</math>
| <math>0,012</math>
| <math>0,00013</math>
| Frut. 
| <math>0,014934</math>
| <math>7,931\times{10}^{-4} </math>
| <math>5,41</math>
| <math>0,011299</math>
| <math>1,349\times{10}^{-3}</math>
| <math>5,55 </math>
|- style="border-top: 2px solid;"
| style="text-align: left;" | <math>0,05487</math>
| Glic. 
| <math>0,014658</math>
| <math>1,1\times{10}^{-4}</math>
| <math>0,73</math>
| <math>0,024265</math>
| <math>8,161\times{10}^{-9}</math>
| <math>0,00</math>
| <math>0,00769</math>
| Frut. 
| <math>0,014619</math>
| <math>6,4576\times{10}^{-4}</math>
| <math>4,41</math>
| <math>0,010981</math>
| <math>4,981\times{10}^{-4}</math>
| <math> 4,53</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="text-align: left;" | <math>0,04476</math>
| Glic. 
| colspan='2' | <math>k_{1}=0,01511</math> <math>ml/mg.min</math> 
| colspan='3' | <math>k_{2}=0,02810</math> <math>ml/mg.min</math>
|                       
| Frut. 
| colspan='2' | <math>k_{1}=0,01454</math> <math>ml/mg.min</math> 
| colspan='3' | <math>k_{2}=0,01332</math> <math>ml/mg.min</math>

|}

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 [[#table-1|1]] foi calculada empregando a equação ([[#eq-12|12]]) <span id='citeF-9'></span>[[#cite-9|[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.

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d = \sqrt{\frac{(1+(\sigma _{k_2}/\mu _{k_2}))^2 + (1+(\sigma _{k_1}/\mu _{k_1}))^2}{2}} - 1    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

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

Neste trabalho foram empregados os dados experimentais reais apresentados na figura [[#img-1|1]], obtidos em <span id='citeF-5'></span>[[#cite-5|[5]]]. Ressalta-se que o desvio padrão dos erros de medição, <math>{\sigma _{C^{exp}}}</math>, não foram reportados, de modo que neste trabalho considerou-se <math>{\sigma _{C^{exp}}}=0,1 mg/ml</math>, 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 [[#img-3a|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.    <div id='img-3a'></div>
<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-posteriore_GlucoseR.png|528px|]]
|[[Image:draft_Samper_522746929-Posteriore_FrutoseR.png|400px|Evolução das cadeias de Markov dos parâmetros de transferência de massa: (a) Glicose e (b) Frutose juntamente a média a posteriori e o valor de referência <span id='citeF-5'></span><span id='citeF-7'></span>[[#cite-5|[5,7]]].]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 3:''' Evolução das cadeias de Markov dos parâmetros de transferência de massa: (a) Glicose e (b) Frutose juntamente a média a posteriori e o valor de referência <span id='citeF-5'></span><span id='citeF-7'></span>[[#cite-5|[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 <math display="inline">N_{MCMC} = 3746</math> para ambas as substâncias e um máximo de <math display="inline">N_{MCMC} = 12842</math> para Glicose e <math display="inline">N_{MCMC} = 13076</math> 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 [[#img-3a|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 [[#img-4|4]] estão apresentados os histogramas juntamente com o valor de referência <span id='citeF-5'></span><span id='citeF-7'></span>[[#cite-5|[5,7]]] dos parâmetros e a média a posteriori.  <div id='img-4a'></div>
<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-hist_Glucose_k1_R.png|384px|]]
|[[Image:draft_Samper_522746929-hist_Glucose_k2_R.png|384px|]]
|-
|[[Image:draft_Samper_522746929-hist_Frutose_k1_R.png|384px|]]
|[[Image:draft_Samper_522746929-hist_Frutose_k2_R.png|384px|Histograma das observações simuladas, normal ajustada, gráfico da densidade de convergência para subfig_a1subfig_b1 Glicose e subfig_c1subfig_d1 Frutose e o valor de referência <span id='citeF-5'></span><span id='citeF-7'></span>[[#cite-5|[5,7]]].]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' Histograma das observações simuladas, normal ajustada, gráfico da densidade de convergência para subfig_a1subfig_b1 Glicose e subfig_c1subfig_d1 Frutose e o valor de referência <span id='citeF-5'></span><span id='citeF-7'></span>[[#cite-5|[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 <span id='citeF-5'></span><span id='citeF-7'></span>[[#cite-5|[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 <math display="inline">95%</math> da distribuição a posteriori estimada.

Os gráficos de dispersão na figura [[#img-5|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'' <span id='citeF-21'></span>[[#cite-21|[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 <math display="inline">k_2</math> e <math display="inline">k_1</math> devido a indicação de dependência dessas condicionais, ou seja, elas têm as mesmas variações, mas <math display="inline">k_2</math> varia com <math display="inline">k_1</math> no sentido em que ambos tendem a aumentar  em conjunto. Este comportamento pode ser observado nos gráficos de sensibilidade na figura [[#img-2|2]] se, ao aumentar ou diminuir a variação dos parâmetros, ambos tendem a compensar o efeito do outro em sentido oposto.  <div id='img-5a'></div>
<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-DispersaoGlR.png|420px|]]
|[[Image:draft_Samper_522746929-Dispersao_frR.png|400px|Gráfico de dispersão com a região de intensidade das estimativas obtidas para as substâncias (a) Glicose e (b) Frutose.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Gráfico de dispersão com a região de intensidade das estimativas obtidas para as substâncias (a) Glicose e (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:

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{C}^{exp} = C^{exp} + \epsilon _{C^{exp}},   \epsilon _{C^{exp}} \sim N(0, \sigma _{C^{exp}})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 [[#table-1|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 (<math display="inline">\mu \pm{1},96\sigma </math>), 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 <span id='citeF-17'></span>[[#cite-17|[17]]]

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{K_i} = \sigma _{C^{exp}} \sqrt{\left[J^TJ\right]^{-1}_{ii}},   i = 1, 2,\dots ,N   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

onde <math display="inline">J_{ji} = \partial C_j/ \partial K_i</math> são os elementos da matriz de sensibilidade calculados pela equação ([[#eq-11|11]]):

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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| </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K_i - 1,96\sigma _{K_i}\leq K_i \leq K_i + 1,96\sigma _{K_i}   i=1,\dots ,N  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 [[#table-2|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 Resultado da  aplicação da Inferência Bayesiana|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 [[#table-2|2]] também apresenta os coeficientes de variação, <math display="inline">\sigma /\mu \times 100%</math>. 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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Média e desvio padrão para <math>k_1</math> e <math>k_2</math> obtidas com 2700 avaliações da função objetivo nos intervalos <math>k_1 = (0,01; 0,05)</math> e <math>k_2 = (0,001; 0,06)</math> com diferentes conjuntos de dados experimentais.
|- style="border-top: 2px solid;"
| style="text-align: left;" |   
|  
| colspan='2' | <math>k_{1}</math>(<math>ml/mg.min</math>)
| 
|-
| colspan='2' | <math>k_{2}</math>(<math>ml/mg.min</math>)
| 
| <math>\mu </math>
| <math>\sigma </math>
| <math>\sigma /\mu \times 100%</math> 
| <math>\mu </math>
| <math>\sigma </math>
|- style="border-top: 2px solid;"
| style="text-align: left;" | <math>\sigma /\mu \times 100%</math>   
| Glicose 
| <math>0,01346</math>
| <math>0,00119</math>
| <math>8,85</math>
| <math>0,02452</math>
| <math>0,00204</math>
| <math>8,35</math>
| Frutose 
| <math>0,01454</math>
| <math>0,00164</math>
| <math>11,27</math>
| <math>0,01171</math>
| <math>0,00123</math>
| <math>10,5</math>
| Glicose 
| <math>0,01391</math>
| <math>0,00115</math>
| <math>8,33</math>
| <math>0,02540</math>
| <math>0,00212</math>
| <math>8,36</math>
| Frutose 
| <math>0,01452</math>
| <math>0,00183</math>
| <math>12,6</math>
| <math>0,01167</math>
| <math>0,00146</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="text-align: left;" | <math>12,5</math>
| Glicose 
| <math>0,01498</math>
| <math>0,00145</math>
| <math>9,67</math>
| <math>0,02583</math>
| <math>0,00256</math>
| <math>9,91</math>
| Frutose 
| <math>0,01445</math>
| <math>0,00214</math>
| <math>14,8</math>
| <math>0,01162</math>
| <math>0,00180</math>
| <math>15,49</math>
| Glicose 
| <math>0,01492</math>
| <math>0,00134</math>
| <math>9,01</math>
| <math>0,02570</math>
| <math>0,00233</math>
| <math>9,10</math>
| Frutose 
| <math>0,01427</math>
| <math>0,00187</math>
| <math>13,11</math>
| <math>0,01147</math>
| <math>0,00157</math>
|- style="border-top: 2px solid;"
| style="text-align: left;" | <math>13,73</math>
|}

<math display="inline">^1</math> Médias e desvios calculados após a remoção de ''outliers''.

As estimativas dos parâmetros, apresentadas na tabela [[#table-2|2]], foram usadas na solução do problema direto, assim como os valores de referência destes parâmetros <span id='citeF-6'></span><span id='citeF-5'></span>[[#cite-6|[6,5]]]. Os valores calculados para a concentração são apresentados nos gráficos da figura [[#img-7|7]].    <div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-Glicose.png|510px|Glicose]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Glicose
|}
<div id='img-7a'></div>
<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Samper_522746929-Frutose.png|510px|Frutose]]
|- style="text-align: center; font-size: 75%;"
| (a) Frutose
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Representação do perfil cromatográfico, (a) Glicose e (b) Frutose, calculados através do problema direto com os parâmetros apresentados na tabela ([[#table-2|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 [[#img-7|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 [[#table-2|2]] utilizou-se o teste de Kruskal-Wallis <span id='citeF-31'></span>[[#cite-31|[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 [[#table-3|3]], onde <math display="inline">H</math> denota a estatística de Kruskal-Wallis e <math display="inline">\mathcal{X}^2</math> é o valor crítico da distribuição ''qui-quadrado''.


==table 3==

<div class="center" style="font-size: 75%;">'''Tabela 3'''. Comparação das amostras dos algoritmos PSO, PCA e MCMC usando o teste de Kruskal-Wallis
</div>


{|  style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:100%;" 
|- style="border-top: 2px solid;"
|  Substância 
| <math>[ml/mg.min]</math>
| <math>H<\mathcal{X}^{2}</math>
| ''p''-valor 
|- style="border-top: 2px solid;border-bottom: 2px solid;"
|rowspan="2" | Glicose
| <math>k_{1}</math>
| <math>1,3881<5,9914</math>
| <math>0,4995</math>
|-
|rowspan="2" | Frutose
| <math>k_{2}</math>
| <math>1,7456<5,9914</math>
| <math>0,4177</math>
|-
| <math>k_{1}</math>
| <math>0,0964<5,9914</math>
| <math>0,9529</math>
|-
| style="border-top: 2px solid;" | <math>k_{2}</math>
| style="border-top: 2px solid;" <math>0,2031<5,9914</math>
| style="border-top: 2px solid;" <math>0,9034</math>
|}

De acordo com as regras do teste, verifica-se na tabela [[#table-3|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, <math display="inline">H<\mathcal{X}^{2}</math>, em todos os casos  obteve-se ''p''-valor maior que a significância <math display="inline">\alpha=0,05</math> considerada.

Na tabela [[#table-4|4]] estão representados os intervalos de confiança de 95% para os parâmetros <math display="inline">k_1</math> e <math display="inline">k_2</math> 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 ([[#eq-13|13]]). Os intervalos sinalizados com (*) foram calculados com a equação ([[#eq-14|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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 Intervalos de confiança de <math>95</math>% para os parâmetros <math>k_1</math> e <math>k_2</math> calculados com a média e desvio padrão das 15 simulações estocásticas e com a abordagem Bayesiana.
|- style="border-top: 2px solid;"
| style="text-align: left;" |   
| colspan='2' | Glicose
|-
| colspan='2' | Frutose
| <math>k_{1} [ml/mg.min]</math>
| <math>k_{2} [ml/mg.min]</math>
| <math>k_{1} [ml/mg.min]</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="text-align: left;" | <math>k_{2} [ml/mg.min]</math>
| <math>(0,01112; 0,01580)</math>
| <math>(0,02050;  0,02853)</math>
| <math>(0,01133; 0,01776)</math>
| <math>(0,00928; 0,01414)</math>
| <math>*(0,01059; 0,01634)</math>
| <math>*(0,02077; 0,02826)</math>
| <math>*(0,01198; 0,01710)</math>
| <math>*(0,00862; 0,01480)</math>
| <math>(0,01164; 0,01618)</math>
| <math>(0,02124;  0,02956)</math>
| <math>(0,01093; 0,01811)</math>
| <math> (0,00880; 0,01454)</math>
| <math>*(0,01103; 0,01678)</math>
| <math>*(0,02165; 0,02915)</math>
| <math>*(0,01196; 0,01708)</math>
| <math>*(0,00858; 0,01476)</math>
| <math> (0,01228; 0,01755)</math>
| <math>(0,02101;  0,03028)</math>
| <math>(0,01060; 0,01794)</math>
| <math> (0,00838; 0,01456)</math>
| <math>\dag (0,01233; 0,01798)</math>
| <math>\dag{0},0211;0,03095)</math>
| <math>\dag (0,01011; 0,01850)</math>
| <math>\dag (0,00813; 0,01508)</math>
| <math>\# (0,01122; 0,01567)</math>
| <math>\# (0,02095; 0,02933)</math>
| <math>\# (0,01088; 0,01972)</math>
|- style="border-top: 2px solid;"
| style="text-align: left;" | <math>\# (0,00834; 0,01572)</math>
|}

<math display="inline">\dag </math>  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 [[#table-4|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 [[#table-2|2]] por exemplo, cada algoritmo (PSO e PCA) realizou <math display="inline">15\times{2700}</math> 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==


<div id="cite-1"></div>
'''[[#citeF-1|[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.

<div id="cite-2"></div>
'''[[#citeF-2|[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.

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''   D.&nbsp;C.&nbsp;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.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''   J.V. Beck, B. Blackwell, and C.R. St. Clair Jr, Inverse Heat Conduction - Ill-Posed Problems, John Wiley & Sons, New York, 1985

<div id="cite-5"></div>
'''[[#citeF-5|[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).

<div id="cite-6"></div>
'''[[#citeF-6|[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.

<div id="cite-7"></div>
'''[[#citeF-7|[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

<div id="cite-8"></div>
'''[[#citeF-8|[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.

<div id="cite-9"></div>
'''[[#citeF-9|[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.

<div id="cite-10"></div>
'''[[#citeF-10|[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.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''   J. Kaipio, E. Somersalo. Statistical and Computational Inverse Problems, New York: Springer,  2005, pp.339.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''   J. Kalyanaraman, et al. Bayesian estimation of parametric uncertainties, quantification and reduction using optimal design of experiments for <math display="inline">CO_2</math> adsorption on amine sorbents. Computers and Chemical Engineering (2015)1-13.

<div id="cite-13"></div>
'''[[#citeF-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.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''   J. Kennedy. The particle swarm: Social Adaptation of knowledge, IEEE International Conference on Evolutionary Computation. (1997)303-308.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''   E.&nbsp;F.&nbsp;P. Luz, J.&nbsp;C.&nbsp;Becceneri, H.&nbsp;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.

<div id="cite-16"></div>
'''[[#citeF-16|[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.

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''   M.&nbsp;N. <math display="inline">\ddot{O}</math>zisik e H.&nbsp;R.&nbsp;B. Orlande, Inverse Heat Transfer: Fundamentals and Applications, Taylor & Francis, New York, 2000

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''   M. Strand. Metropolis-Hastings Markov Chain Monte Carlo, Chapman University.  (2009)1-9.

<div id="cite-19"></div>
'''[[#citeF-19|[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

<div id="cite-20"></div>
'''[[#citeF-20|[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.

<div id="cite-21"></div>
'''[[#citeF-21|[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.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''   W.&nbsp;F. Sacco , C.&nbsp;R.&nbsp;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.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''   W.&nbsp;F. Sacco, C.&nbsp;R.&nbsp;E. de Oliveira. A new stochastic optimization algorithm based on particle collisions metaheuristic. In <math display="inline">6^{th}</math> World Congresses of Structural and Multidisciplinary Optimization, Rio de Janeiro, Brazil, (2005)1-6.

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''   W. F.&nbsp;Sacco, C.&nbsp;M.&nbsp;F. Lapa, C.&nbsp;M.&nbsp; N.&nbsp;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.

<div id="cite-25"></div>
'''[[#citeF-25|[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.

<div id="cite-26"></div>
'''[[#citeF-26|[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.

<div id="cite-27"></div>
'''[[#citeF-27|[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).

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''   G.&nbsp;J.&nbsp;Simões, N.&nbsp;F.&nbsp;F.&nbsp;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.<math display="inline"><</math>''http://dx.doi.org/10.1016/j.rimni.2014.07.001''<math display="inline">></math>

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''   A.&nbsp;M. Stuart. Inverse problems: a Bayesian perspective. Acta Numerica 19 (2010)451-559.

<div id="cite-30"></div>
'''[[#citeF-30|[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.

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''   N. Thomson. Understanding ANOVA the apl way. Em APL ’93 International Conference on Array-Programming Languages, New York, NY, USA, (1993)295-303.

<div id="cite-32"></div>
'''[[#citeF-32|[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.

<div id="cite-33"></div>
'''[[#citeF-33|[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.

<div id="cite-34"></div>
'''[[#citeF-34|[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.