Resolução de um Problema Inverso Diferencial Fracionário no Processo de Fermentação Batelada Usando o Método da Colocação Ortogonal e o Algoritmo de Busca Fractal Estocástica

RESUMO

A presente contribuição tem como objetivo desenvolver uma ferramenta numérica para a resolução de problemas inversos modelados por equações diferenciais ordinárias com ordem fracionária. Esta consiste da associação entre o Método da Colocação Ortogonal no contexto fracionário com o algoritmo de Busca Fractal Estocástica. Os resultados obtidos com a extensão do Método da Colocação Ortogonal em funções matemáticas demonstraram a habilidade desta estratégia em comparação com outras abordagens numéricas. Para fins de ilustração, um problema inverso que consiste na determinação dos parâmetros de um modelo e da ordem fracionária do processo de fermentação da enzima lacase é proposto e resolvido. Em relação a este estudo pode-se concluir que o aumento do número de graus de liberdade (ordem fracionária é mais uma variável de projeto) aumenta a chance de um melhor ajuste do modelo aos pontos experimentais.

Palavras-chave:
problemas inversos; equações diferenciais ordinárias fracionárias; busca fractal estocástica; fermentação batelada

ABSTRACT

This contribution aims to develop a numerical tool for solving inverse problems modeled by ordinary differential equations with fractional order. This consists of association between the Orthogonal Collocation Method in fractional context and the Stochastic Fractal Search algorithm. The results obtained with the extension of the Orthogonal Collocation Method in mathematical functions demonstrated the ability of this strategy in comparison with other numerical approaches. For the purposes of illustration, an inverse problem for the determination of model parameters and fractional order in laccase enzyme fermentation process considering real experimental data is proposed and solved. In relation to this study, it can be concluded that the increase in number of freedom degrees (fractional order is considered as a new design variable) increases the chance of a better fit between the model and experimental data.

Keywords:
inverse problems; fractional ordinary differential equation; stochastic fractal search; batch fermentation

1 INTRODUÇÃO

Nas últimas décadas, o aprimoramento/desenvolvimento de métodos analíticos e numéricos empregados para a análise do cálculo fracionário têm atraído cada vez mais adeptos. O cálculo fracionário é uma generalização da diferenciação e integração ordinárias para uma ordem não inteira arbitrária. O interesse por seu estudo se deve ao fato de que vários fenômenos observados na natureza não podem ser explicados pelo tradicional cálculo inteiro. Por exemplo, a relação entre a força e o deslocamento em uma mola ideal que, sob pequenas deformações em materiais elásticos, é linear. Por outro lado, em um amortecedor ideal a força é proporcional à velocidade de extensão ou compressão da mesma, i.e., a força está relacionada à primeira derivada da deformação. Todavia, qual lei física governa os materiais com propriedades mecânicas intermediárias, i.e., entre a mola ideal e o amortecedor ideal? ⁶ . Objetivando responder esta e outras questões, inúmeras pesquisas têm sido desenvolvidas com o intuito de modelar e analisar o comportamento de sistemas não lineares por meio de cálculo fracionário ¹⁴ .

Historicamente, a primeira aplicação do cálculo fracionário é atribuída a Niels Henrik Abel em 1826 ¹⁶ . Este pesquisador apresentou uma abordagem dedutiva empregada para encontrar a solução de uma equação integral resultante da formulação de um problema de engenharia mecânica (denominado de problema da isócrona). Desde então, várias aplicações no contexto fracionário podem ser encontradas, dentre as quais pode-se citar: materiais viscoelásticos ² , processos eletroquímicos ⁹ , polarização dielétrica ²² , processamento de sinais ¹⁷ , teoria de controle ¹⁹ , advecção e dispersão de solutos em meios porosos ³ e caos ¹⁵ .

Desde o seu surgimento, inúmeras sugestões para interpretação do cálculo fracionário podem ser encontradas na literatura ^{8), (19} . Todavia, a maioria destas são abstratas e sem ligação com aspectos físicos ⁸ . Uma interpretação bastante útil para a interpretação física para o cálculo fracionário é baseada no conceito de memória. Para entender este conceito considere um sistema em que a saída em cada instante de tempo depende apenas da entrada no mesmo intervalo de tempo. Neste caso, esse sistema é considerado sem memória. Por outro lado, quando o sistema precisa lembrar os valores anteriores da entrada para determinar o valor atual da saída, esse é denominado de sistemas com memória ⁸ . Do ponto de vista físico, o que é a memória e como esta é definida em um sistema depende de um profunda compreensão dos fenômenos envolvidos. Além disso, não existem regras e métodos para selecionar o tipo de provável contribuição fracionária durante a etapa da modelagem de um sistema. Entretanto, o uso de pontos experimentais pode ser útil para a caracterização destes modelos ¹³ . Como destacado por Li e colaboradores ¹¹ , devido à propriedade não local de uma derivada fracionária, esta pode ser usada para descrever sistemas de engenharia com dinâmicas complexas por envolverem uma longa memória no tempo.

Na literatura especializada, inúmeros métodos analíticos e numéricos podem ser encontrados para a resolução de uma Equação Diferencial Ordinária Fracionária (EDOF). Dentre estes, Lin e Liu ¹³ introduziram um método numérico baseado no conceito de múltiplos passos para resolver EDOFs lineares. Li e colaboradores ¹¹ propuseram novas abordagens numéricas baseadas no método Simpson para a resolução de EDOF gerais. Li e Zeng ¹² estudaram a estabilidade e a convergência de diferentes metodologias para a resolução de EDOFs. Para essa finalidade, os Métodos de Euler, Adams e de alta ordem com base na fórmula da quadratura foram analisados. Oldham e Spanier ¹⁸ introduziram abordagens numéricas baseadas nas definições de Grünwald-Letnikov e Riemann-Liouville para aproximar as derivadas fracionárias. Podlubny ¹⁹ introduziu propôs uma abordagem numérica para aproximar derivadas com ordens arbitrárias usando a definição de Riemann-Liouville. Deng ⁷ propôs um procedimento para o uso do princípio de memória usando um método preditor-corretor para a resolução de EDOFs. Yuan e Agrawal ²⁶ desenvolveram uma nova abordagem numérica para minimizar o efeito prejudicial de propriedades não locais no contexto fracionário.

No contexto dos problemas inversos diferenciais fracionários, o número de trabalhos tem aumentado significativamente nos últimos anos. Zheng e Wei 2011 ²⁷ propuseram um problema inverso para uma equação de difusão tempo-fracionária em um domínio semi-infinito unidimensional. Para essa finalidade, o perfil de temperatura e o fluxo de calor formam determinados considerando pontos experimentais sintéticos e uma nova metodologia para lidar com o problema do mal condicionamento numérico. Xiong e colaboradores ²⁵ propuseram um problema inverso para uma equação de difusão fracionária mal condicionada numericamente. Tal problema foi obtido a partir da equação clássica de difusão substituindo a derivada de tempo de primeira ordem pela derivada fracionária de Caputo, bem como o uso de pontos experimentais sintéticos. Chi e colaboradores ⁵ determinaram o termo fonte a partir da formulação e resolução de um problema inverso para a equação de advecção-dispersão fracionária considerando pontos experimentais sintéticos. Para essa finalidade, um esquema baseado em diferenças finitas associado ao conceito de derivadas fracionárias de Grünwald foi proposto. Leitoles ¹⁰ propôs e resolveu um problema inverso para fins da determinação dos parâmetros cinéticos em um sistema reacional, bem como da ordem fracionária considerando pontos experimentais reais.

Diante do que foi apresentado, a presente contribuição tem por objetivo desenvolver uma ferramenta numérica para a resolução de problemas inversos sujeitos à equações diferenciais ordinárias que apresentam ordem fracionária. Em linhas gerais, esta consiste da associação entre o Método da Colocação Ortogonal (MCO), estendido para o contexto fracionário, com o algoritmo de Busca Fractal Estocástica (BFE). A metodologia proposta é empregada para determinar os parâmetros cinéticos e a ordem da derivada fracionária no processo de fermentação da enzima lacase. Além disso, também é analisada a influência da ordem fracionária nos perfis simulados para o processo em questão. Este trabalho está estruturado como segue. Na seção 2 é apresentada uma breve descrição matemática do processo de fermentação da lacase. A descrição do MCO para ambos os contextos inteiro e fracionários são apresentados na seção 3. Uma breve descrição do algoritmo de BFE é apresentado na seção 4. A metodologia e os resultados e discussão são apresentados nas seções 5 e 6, respectivamente. A última seção apresenta as conclusões deste trabalho.

2 MODELAGEM MATEMÁTICA DO PROBLEMA DE INTERESSE

Nesta contribuição considera-se o processo de fermentação da enzima lacase. Esta apresenta uma série de aplicações industriais, dentre as quais pode-se citar o seu uso na indústria de alimentos e no tratamento de efluentes industriais. Especificamente na indústria de alimentos, esta enzima é adicionada para ressaltar as cores, o aroma, o gosto, a estabilidade nos sucos das frutas e nas bebidas alcoólicas fermentadas ¹ . Matematicamente, este processo pode ser representado pelos seguintes balanços de massa ¹:

$\frac{dX}{dt}=\frac{{\mu }_{\max }SX}{\left(K_s+S\right)}\left(X\right(0)=X_0)$ (2.1)

$\frac{dS}{dt}=-\frac{1}{Y_{XS}}\frac{dX}{dt}\left(S\right(0)=S_0)$ (2.2)

$\frac{dP}{dt}=\left\{\begin{array}{l}{K}_1\left(1-exp\left(-(t-0,1)/0,1\right)\right)S-K_{p}Pset\le t_s\\ \left(K_1-K_2\right)\left(1-exp\left(-(t-0,1)/0,1\right)\right)Sset>t_s\end{array}\right.\left(P\right(0)=P_0)$ (2.3)

em que t é o tempo (dias), X é a concentração de biomassa (g/L), P é a concentração da enzima lacase (U/L), S é a concentração de substrato (g/L), µmax e Ks são as constantes relacionadas ao modelo de inibição considerado, Y _X/S é o rendimento de biomassa em relação ao substrato, K ₁, K ₂ e K _p são parâmetros que ponderam os modelos no que tange o produto e t _s é o tempo de atuação de cada um destes modelos. O subscrito 0 representa a condição inicial associada ao vetor de variáveis dependentes ( X , S e P ).

Se todos os parâmetros ( µmax , K _s , Y _X/S , K ₁, K ₂, K _p , t _s ) e as condições iniciais ( X ₀, S ₀ e P ₀) do modelo forem conhecidos, tem-se um problema de simulação ou problema direto (em oposição a problema inverso). Por outro lado, se todos ou qualquer combinação destes não forem conhecidos, deve-se formular e resolver um problema inverso considerando pontos experimentais para a determinação do conjunto de incógnitas (variáveis de projeto). Mais detalhes sobre a formulação do problema inverso diferencial fracionário proposto neste trabalho serão apresentados na seção 5.

3 MÉTODO DA COLOCAÇÃO ORTOGONAL

Em linhas gerais, o MCO consiste na definição de uma função para aproximar a solução do problema em análise em determinados pontos do domínio de interesse, denominados de pontos de colocação. Neste caso, obriga-se que a equação original coincida com a solução aproximada nestes pontos. Apesar da escolha da função de aproximação e dos pontos de colocação poderem ser realizadas de forma aleatória, tem utilizado como função de aproximação o Polinômio de Lagrange (PL) e como pontos de colocação as raízes de um polinômio ortogonal. Cabe ressaltar que a escolha pelo PL se deve ao fato de que este possibilita a redução do custo computacional associado com a necessidade de obtenção de derivadas em relação a outros tipos de aproximações ²³ .

Conforme destacado anteriormente, o objetivo desta contribuição é a resolução de um problema inverso diferencial fracionário considerando como ferramenta para a simulação do modelo fracionário o MCO estendido para o contexto fracionário. Antes de apresentar a metodologia proposta, faz-se necessário discutir aspectos gerais do MCO para o contexto inteiro, conforme a próxima subseção.

3.1 MCO para o Contexto Inteiro

O MCO é baseado na definição de uma função de aproximação (geralmente uma função polinomial), na qual sua solução numérica é avaliada considerando um determinado número de pontos dentro do domínio de interesse, i.e., o Número de Pontos de Colocação ou NPC²⁴ . A equação original deve satisfazer a função de aproximação nos pontos considerados, bem como nas condição inicial e de contorno (caso o problema em questão seja de valor no contorno). Uma função ortogonal é usada para determinar a posição ótima desses pontos no domínio X . Para esse objetivo, considere a seguinte relação recursiva ^{23), (24}:

${\Pi }_i^{(\psi ,\eta )}\left(X\right)=(X+{\psi }_i){\Pi }_{i-1}^{}\left(X\right)+{\eta }_i{\Pi }_{i-2}^{}\left(X\right)$ (3.1)

em que Π é uma função polinomial e ψ e η são coeficientes definidos a partir da aproximação considerada.

As raízes desta equação são obtidas considerando um polinômio ortogonal de grau NPC , o peso W ( X ) e a seguinte condição de Galerkin ²⁴:

$\underset{0}{\overset{1}{\int }}W\left(X\right)(\psi X+\eta ){\Pi }_i^{(\psi ,\eta )}\left(X\right)dX=0,i=0,...,NPC-2$ (3.2)

Multiplicando a Eq.(3.1) por Π_i−2 e integrando, a seguinte relação é obtida:

${\eta }_i=-\frac{{\int }_0^{1}XW\left(X\right){\Pi }_{i-1}^{}\left(X\right){\Pi }_{i-2}^{}\left(X\right)dX}{{\int }_0^{1}W\left(X\right){\Pi }_{i-2}^2\left(X\right)dX}$ (3.3)

Analogamente, multiplicando a Eq.(3.1) por Π_i−1 e integrando, o parâmetro ψ pode ser obtido, i.e.:

${\psi }_i=-\frac{{\int }_0^{1}XW\left(X\right){\Pi }_{i-1}^2\left(X\right)dX}{{\int }_0^{1}W\left(X\right){\Pi }_{i-1}^2\left(X\right)dX}$ (3.4)

Considerando, por exemplo, ${\Pi }_{i-2}=0,{\Pi }_{i-1}=1,W\left(X\right)=1$ e ${\eta }_1=0$ , o procedimento acima pode ser usado para calcular ψ e η e, consequentemente, determinar o polinômio ortogonal em NPC , i.e., ${\Pi }_{NPC}^{(\psi ,\eta )}\left(X\right)$ ^{23), (24} . As raízes deste polinômio são tomadas como sendo os pontos de colocação.

Sabe-se que a função de aproximação e os pontos de colocação podem ser escolhidos usando abordagens diferentes. A metodologia do polinômio de Lagrange (PL) tem sido tradicionalmente usada como uma função de aproximação. Essa escolha se deve à redução do custo computacional associado à aproximação numérica dos derivativos em comparação com outras aproximações ^{23), (24} . Considerando o conjunto de pontos de dados $(X_1,Y_1),(X_2,Y_2),...,(X_{NPC+1},Y_{NPC+1})$ , uma fórmula de interpolação que passa por esses pontos (um polinômio de interpolação de NPC -ésimo grau) é dada por:

$Y_{NPC}\left(X\right)=\sum _{i=1}^{NPC+1}{Y}_i{l}_i\left(X\right)$ (3.5)

em que $l_i\left(X\right)$ é o polinômio de interpolação de Lagrange definido como:

$l_i\left(X\right)=\prod _{j=1}^{NPC+1}\frac{X-X_j}{X_i-X_j}$ (3.6)

Se o subscrito i for igual a $j,l_i\left(X\right)$ é igual a 1. Caso contrário, $l_i\left(X\right)$ é igual a 0. Como essa aproximação é uma função contínua, ela pode ser diferenciada e integrada. Assim, a primeira derivada para uma raiz específica X _j pode ser expressa como:

$\frac{d{Y}_{NPC}\left(X_j\right)}{dX}=\sum _{i=1}^{NPC+1}{Y}_i\frac{d{l}_i\left(X_j\right)}{dX},j=1,2,...,NPC+1$ (3.7)

Substituindo essas aproximações no modelo original, é possível obter expressões (residuais) para cada ponto de colocação. Estes resíduos devem ser minimizados para cada raiz do polinômio ortogonal, i.e., devem ser zerados em cada ponto de colocação, bem como devem satisfazer as condições iniciais e de contorno (caso o problema seja de valor no contorno). O sistema algébrico resultante, geralmente não linear, deve ser resolvido considerando uma técnica numérica apropriada, como por exemplo, o Método Newton.

A seguir é apresentado um consolidado dos passos necessários para a aplicação do MCO para a resolução de equações diferenciais ordinárias com ordem inteira:

1. Definir os parâmetros de entrada: domínio do problema, função peso para a determinação do polinômio ortogonal, grau NPC da aproximação (o polinômio tem $NPC+1$ coeficientes);

2. Calcular o polinômio ortogonal de grau $NPC+1-NC$ , onde NC é o número de condições que precisam ser satisfeitas;

3. Calcular as NPC raízes (pontos de colocação) do polinômio ortogonal;

4. Determinar as equações (resíduo nos pontos de colocação);

5. Resolver o sistema de equações resultantes;

6. Verificar se a solução obtida não é modificada com o aumento do grau da aproximação.

É importante ressaltar que, como o MCO apresentado é fundamentado no uso de um polinômio ortogonal, a variável independente no modelo diferencial ordinário deve estar definida no intervalo [0 1] para que o MCO possa ser empregado. Caso isso não aconteça, inicialmente deve-se realizar uma mudança de variável de modo que o modelo seja sempre integrado no intervalo [0 1].

Finalmente, cabe ressaltar que a qualidade do resultado obtido é função da aproximação considerada, sendo que um aumento excessivo do grau da mesma não, necessariamente, implica na melhora da qualidade da solução obtida. Além disso, esse aumento pode ocasionar um comportamento oscilatório nas proximidades de regiões onde a solução não experimenta variações pronunciadas ²³ .

3.2 Extensão do MCO para o Contexto Fracionário

Na subseção anterior foi apresentado o algoritmo para a resolução de uma equação diferencial ordinária com ordem inteira. Para a sua aplicação em um problema em que o termo diferencial apresenta ordem fracionária, é necessária uma pequena adaptação no algoritmo apresentado. Neste caso, a partir da definição da derivada fracionária a ser considerada, durante a caracterização das derivadas considerando o PL no algoritmo apresentado, deve-se empregar uma das definições de derivadas fracionárias para a representação deste termo. Com esta mudança, o algoritmo que antes tratava apenas de derivadas inteiras, agora pode lidar com derivadas fracionárias. Para validar a metodologia proposta, a mesma será aplicada em problemas puramente matemáticos e que apresentam solução analítica, conforme apresentado e discutido na seção 5.

4 BUSCA FRACTAL ESTOCÁSTICA

Para resolver o problema inverso de interesse neste trabalho, considera-se o algoritmo de Busca Fractal Estocástica (BFE) como ferramenta de otimização. Esta recente técnica de otimização inspirada na natureza foi proposta por Salimi ²⁰ . Em linhas gerais, esta estratégia baseia-se no uso de fractais para imitar o fenômeno de deslocamento de partículas (candidatos a solução do problema de otimização) no processo de difusão. Neste movimento, as partículas tendem a explorar o espaço de busca considerando uma propriedade de difusão empregada em fractais aleatórios ⁴ . Segundo Salimi ²⁰ , o algoritmo de BFE faz uso de fractais aleatórios via emprego do método de agregação limitada por difusão como metodologia de busca.

No algoritmo de BFE, dois processos são aplicados para gerar novos candidatos, a saber, o processo de difusão e o processo de atualização. No primeiro, cada partícula se difunde em torno de sua posição atual para garantir a propriedade de exploração. Este procedimento evita que o algoritmo fique preso em mínimos locais, além de aumentar a chance de encontrar o mínimo global. No segundo processo, uma determinada partícula tem a sua posição atualizada em função das localizações atuais de outras partículas vizinhas. Além da exploração proposto no algoritmo de BFE, o mesmo ainda emprega estratégias puramente aleatórias como metodologias para a atualização dos candidatos em potencial. Uma breve descrição do algoritmo de BFE é apresentada a seguir ²⁰ .

O processo de difusão consiste na geração de caminhos aleatórios gaussianos para gerar pontos em torno de cada partícula até que um número máximo de difusões, pré-determinado pelo usuário, seja alcançado. Existem dois tipos de caminhos que podem ser seguidos, a saber, GW ₁ e GW ₂, definidos como:

$G{W}_1=Gaussian\left({\mu }_{BP},\delta \right)+\left(rand(0,1)\times P_{best}-rand\times P_i\right)$ (4.1)

$G{W}_2=Gaussian\left({\mu }_P,\delta \right)$ (4.2)

em que P _i e P _best são as posições referentes a i -ésima e a melhor partículas, respectivamente ( $i=1,2,...,NP$ , onde NP é o tamanho da população), rand (0,1) é um número aleatório gerado no intervalo $[0,1].{\mu }_{BP}$ e µP são parâmetros definidos pelo usuário. O desvio padrão δ é ajustado dinamicamente com base no número da geração G:

$\delta =\left|\frac{\log \left(G\right)}{G}\times \left(P_i-P_{best}\right)\right|$ (4.3)

No algoritmo de BFE o processo para a atualização de um candidato em potencial emprega dois procedimentos estatísticos. No primeiro processo de atualização cada posição de partícula é atualizada como segue:

$P_i^{new}\left(j\right)=\left\{\begin{array}{cc}{P}_{1r}\left(j\right)-rand(0,1)\times \left(P_{2r}\left(j\right)-P_i\left(j\right)\right),& \mathrm{se}{\gamma }_i\le rand(0,1)\\ P_j\left(j\right),& \mathrm{se}{\gamma }_i>rand(0,1)\end{array}\right.$ (4.4)

em que j pertence ao conjunto $\{1,2,...,D\}$ que define o número de variáveis de projeto $D.P_i^{new}$ é a nova posição da i -ésima partícula P _i , computada em função das posições das partículas P _1r e P _2r , selecionadas aleatoriamente, γ _i é a probabilidade de seleção da partícula P _i , calculada como:

${\gamma }_i=1-\frac{ran{k}_i}{NP}$ (4.5)

em que o operador rank _i representa o ordenamento da população de partículas em relação ao valor da função objetivo (o ordenamento é realizado da melhor para a pior solução, i.e., a partícula com pior valor de função objetivo tem um menor valor para o rank _i ). No segundo processo de atualização, a probabilidade γ _i é recalculada e usada para determinar se uma partícula P _i deve ser atualizada, conforme a seguinte relação:

$P_i^{new}\left(j\right)=\left\{\begin{array}{ll}{P}_i\left(j\right)-rand(0,1)\times \left(Q_{1r}\left(j\right)-P_{best}\left(j\right)\right),& \mathrm{se}rand(0,1)\le 0.5\\ P_i\left(j\right)+rand(0,1)\times \left(Q_{1r}\left(j\right)-Q_{2r}\left(j\right)\right),& \mathrm{se}rand(0,1)>0.5\end{array}\right.$ (4.6)

em que Q _1r e Q _2r são duas posições das partículas, selecionadas aleatoriamente. Este procedimento deve ser executado até que um determinado critério de parada seja satisfeito. Neste algoritmo considera-se o número máximo de gerações como critério de parada. Mais detalhes sobre o desenvolvimento teórico do algoritmo BFE podem ser encontrados no trabalho de Salimi ²⁰ .

5 METODOLOGIA

5.1 Resolução dos Problemas Diretos

Para validar a metodologia proposta para resolver uma equação diferencial ordinária fracionária e detalhada anteriormente, dois estudos de caso puramente matemáticos e que apresentam solução analítica conhecida são considerados. Neste cenário, a influência da ordem fracionária e do número de pontos de colocação ( NPC ) são analisados. Para avaliar o termo fracionário será utilizado a derivada do tipo Riemann-Liouville ${(}_{RL}{D}^{\alpha }f\left(t\right))$ com ordem fracionária $\alpha (\alpha >0)$ aplicada a uma função genérica f no domínio $(0,\infty )$ . Esta é definida como:

${}_{RL}{D}^{\alpha }f\left(t\right)=\frac{d^n}{d{t}^n}\left(\frac{1}{\Gamma (n-\alpha )}\underset{0}{\overset{t}{\int }}{\left(t-\tau \right)}^{n-\alpha -1}f\left(\tau \right)d\tau \right)$ (5.1)

em que t é uma variável independente e $n=\left[\alpha \right]+1$ e [ α ] é um operador que representa a parte inteira de α .

5.2 Resolução dos Problemas Inversos

Após a validação da abordagem proposta para a integração do modelo diferencial fracionário, a mesma será associada ao algoritmo de BFE para a resolução de dois problemas inversos, a saber, um puramente matemático em que serão gerados pontos experimentais sintéticos e um outro cujo objetivo é a resolução de um problema inverso diferencial fracionário aplicado ao processo de fermentação batelada para a produção da enzima lacase. Para essa finalidade, alguns pontos devem ser mencionados:

• Para a formulação do problema inverso matemático considera-se minimizar o funcional FO ₁ a partir da determinação da ordem fracionária α no MCO, conforme a seguinte equação:

$F{O}_1\equiv \sum _{i=1}^N(Y_i^{\exp }-Y_i^{cal}{)}^2$ (5.2)

• em que N é o número de pontos experimentais, $Y_i^{\exp }$ e $Y_i^{cal}$ representam os pontos experimentais e os calculados numericamente pelo MCO, respectivamente.

• • Como não são conhecidos pontos experimentais para este problema em particular, pontos experimentais sintéticos são gerados conforme a seguinte relação:

$Y_i^{\exp }=Y_i+\xi {\rm Y}$ (5.3)

• em que ξ é o desvio padrão associado aos erros de medida e ϒ é um número pseudo-aleatório gerado com uma distribuição gaussiana de média zero e o desvio padrão igual a um. Para determinar os pontos experimentais sintéticos, a solução analítica ( Y _i ) conhecida será utilizada para essa finalidade.

• Para a caracterização do problema inverso diferencial fracionário referente ao processo fermentativo, as derivadas de ordem inteira (Eq.(2.1)-Eq.(2.3)) são substituídas por derivadas fracionárias com ordem igual a α , i.e., $dX/dt\to d^{\alpha }X/d{t}^{\alpha },dS/dt\to d^{\alpha }S/d{t}^{\alpha }$ e $dP/dt\to d^{\alpha }P/d{t}^{\alpha }$ . Neste caso, para que o modelo modificado continue tendo consistência dimensional, visto que a generalização da ordem fracionária interfere nas unidades, é necessário corrigir as dimensões de alguns dos parâmetros, conforme apresentado a seguir no item que define o espaço de projeto no problema de otimização;

• Para aplicar o MCO no contexto fracionário será considerada definição da derivada do tipo Riemann-Liouville, conforme a Eq.(5.1);

• Para resolver o sistema de equações resultantes da aplicação do MCO, o Método de Newton é considerado como metodologia numérica para a resolução do sistema não linear;

• O problema inverso diferencial fracionário consiste da determinação dos parâmetros do modelo ( µ _max , K _s , Y _X/S , K ₁, K ₂, K _p e t _s ) e da ordem fracionária ( α ) que minimize o somatório do erro absoluto quadrático entre os pontos experimentais e simulados ponderado pelo maior valor experimental observado, i.e.:

$F{O}_2\equiv \sum _{i=1}^M\frac{{\left(X_i^{\exp }-X_i^{cal}\right)}^2}{{\left(X_{\max }\right)}^2}+\sum _{i=1}^M\frac{{\left(S_i^{\exp }-S_i^{cal}\right)}^2}{{\left(S_{\max }\right)}^2}+\sum _{i=1}^M\frac{{\left(P_i^{\exp }-P_i^{cal}\right)}^2}{{\left(P_{\max }\right)}^2}$ (5.4)

• em que FO ₂ é o funcional que deve ser minimizado, M é o número de pontos experimentais, $X_i^{\exp }$ e $X_i^{cal}$ representam os i -ésimos valores experimentais e simulados pelo modelo para a concentração de células, $S_i^{\exp }$ e $S_i^{cal}$ representam os i -ésimos valores experimentais e simulados pelo modelo para a concentração de substrato e $P_i^{\exp }$ e $P_i^{cal}$ representam os i -ésimos valores experimentais e simulados pelo modelo para a concentração de produto. X _max, S _max e P _max representam os maiores valores experimentais observados para X , S e P , respectivamente.

• Espaço de projeto considerado: $0,001\le {\mu }_{max}\le 5\left({\mathrm{dia}}^{-\mathrm{\alpha }}\right)$ ; $0,001\le K_s\le 10(\mathrm{g}/\mathrm{L})$ ; $0,002\le Y_{X/S}\le 3(\mathrm{g}/\mathrm{g})$ ; $30\le K_1\le 30000(\mathrm{U}{\mathrm{g}}^{-1}{\mathrm{dia}}^{-\alpha })$ ; $50\le K_2\le 50000(\mathrm{U}{\mathrm{g}}^{-1}{\mathrm{dia}}^{-\alpha })$ ; $0,001\le K_p\le 2\left({\mathrm{dia}}^{-\alpha }\right)$ ; $6\le t_s\le 11\left(\mathrm{dia}\right)$ e $0,9\le \alpha \le 1$ ;

• Os pontos experimentais considerados para a formulação do problema inverso diferencial fracionário proposto são apresentados na tabela 1 ¹: Estes representam os pontos experimentais obtidos com a fermentação para a produção da lacase por Tramares versicolor em cultura submersa descontínua. Neste processo estudado por Alves ¹ foram analisados o efeito da concentração inicial de extrato de malte e o efeito do pH do meio na atividade da lacase.

• O vetor de condições iniciais, i.e., para o tempo inicial ( $t=1$ Dia), corresponde a $[X_0=0,1S_0=2P_0=0]$ , conforme a tabela 1 ;

• Para avaliar o algoritmo de BFE, duas configurações são consideradas, a saber, BFE¹ e BFE², respectivamente. BFE¹ apresenta os seguintes parâmetros: tamanho da população (50), número de gerações (500), número máximo de difusões permitidas em torno de cada partícula (2), primeiro caminho gaussiano - Eq.(4.4). BFE² apresenta os mesmos parâmetros considerados em BFE1, todavia, faz uso do segundo caminho gaussiano - Eq.(4.6);

• Para fins de comparação considera-se o algoritmo de Evolução Diferencial ²¹ com os seguintes parâmetros: tamanho da população (50), número de gerações (500), taxa de perturbação (0,8) e probabilidade de cruzamento (0,8);

Tabela 1:
Dados experimentais para o processo de fermentação da enzima lacase.

t (Dia)	X (g/L)	S (g/L)	P (U/L)
1	0,1	2,0	0
2	0,1	2,0	125,3
3	0,3	1,9	1233,8
4	0,4	1,6	2699,1
5	0,7	1,3	5656,6
6	0,8	1,0	7089,1
7	1,1	0,8	5984,7
8	1,1	0,5	2897,9
9	1,2	0,4	1175,5
10	1,2	0,2	112,2
11	1,2	0,1	32,5
12	1,1	0,1	1,1

• Cada algoritmo foi executado 20 vezes e os melhores resultados são apresentados nas tabelas a seguir. Além disso, com os parâmetros considerados, são necessárias $50+50\times 500$ avaliações da função objetivo em cada execução de cada algoritmo. O tempo de processamento médio será avaliado considerando um microcomputador Intel Core i7 com 8 GB de memória.

6 RESULTADOS E DISCUSSÃO

Na presente seção são apresentados os resultados obtidos com a metodologia proposta. Inicialmente são apresentados e discutidos os resultados referentes a simulação de duas equações diferenciais ordinárias fracionárias para fins de comparação entre o MCO para o contexto fracionário e outras abordagens numéricas. Em seguida, são apresentados e discutidos os resultados obtidos com a resolução de dois problemas inversos usando o algoritmo de BEF. Finalmente, na última subseção é apresentada a influência da ordem fracionária nos perfis obtidos considerando a simulação do processo de fermentação da lacase.

6.1 Problema Direto

Para avaliar o desempenho da metodologia proposta para a resolução de equações diferenciais ordinárias fracionárias, dois modelos puramente matemáticos são estudados. O primeiro, proposto por Lin e Liu ¹³ e estudado por Li e Zeng ¹² , é dado pelo seguinte equacionamento:

$D^{\alpha }y\left(t\right)+y^2\left(t\right)=f(t,\alpha ),0\le \alpha \le 2,0<t\le 2$ (6.1)

em que:

$f(t,\alpha )=\frac{\Gamma \left(6\right)}{\Gamma (6-\alpha )}{t}^{5-\alpha }-\frac{3\Gamma \left(5\right)}{\Gamma (5-\alpha )}{t}^{4-\alpha }+\frac{2\Gamma \left(4\right)}{\Gamma (4-\alpha )}{t}^{3-\alpha }+(t^5-3t^4+2t^3{)}^2$ (6.2)

Neste caso, se α pertencer ao intervalo [0 1], a condição inicial para este modelo é igual a $y\left(0\right)=0$ . Todavia, se α for maior do que a unidade e menor do que dois, o problema passa a ter duas condições iniciais, i.e., $y\left(0\right)=y\text{'}\left(0\right)=0$ . Conforme reportado por Lin e Liu ¹³ , a solução analítica para este estudo considerando α igual a 3/2 é igual a $y\left(t\right)=t^5-3t^4+2t^3$ .

A tabela 2 apresenta os erros absolutos obtidos com a aplicação do MCO estendido para o contexto fracionário considerando diferentes valores para o parâmetro NPC . Nesta tabela também são apresentados os resultados obtidos por duas metodologias numéricas. A primeira, denominada de Método p -HOFLMSM (M p -H), foi proposta por Lin e Liu ¹³ e consiste do uso de uma fórmula de quadratura associada a derivada do tipo de Riemman-Liouville para discretizar a EDOF. Já a segunda estratégia, denominada de Método de Alta Ordem (MAO), foi proposta por Li e Zeng ¹² e é baseada no uso fórmula de Gronwall discreta generalizada para aproximar as derivadas fracionárias. Para essa finalidade, em ambas estratégias (M p -H e MAO) foi adotado um passo de integração da ordem de 0,2 e ordem igual a 3 $(p=3)$ . Para este passo de integração, um sistema com 10 equações algébricas deve ser resolvido.

Tabela 2:
Erros absolutos obtidos para a primeira EDOF considerando $\alpha =3/2$ e diferentes abordagens numéricas.

t	M p -H 13	MAO 12	$NPC=5$	$NPC=10$	$NPC=15$	$NPC=20$
0,2	3,22×10-6	3,22×10-6	1,03×10-5	3,19×10-8	7,28×10-10	9,08×10-14
0,4	2,91×10-5	2,91×10-5	1,11×10-5	3,20×10-8	6,24×10-10	6,98×10-14
0,6	3,97×10-5	3,97×10-5	1,15×10-5	3,11×10-8	6,76×10-10	6,11×10-14
0,8	3,48×10-5	3,49×10-5	1,10×10-5	2,55×10-8	5,34×10-10	5,63×10-14
1,0	1,54×10-5	1,54×10-5	1,31×10-5	3,45×10-8	4,98×10-10	4,65×10-14
1,2	1,81×10-5	1,82×10-5	1,47×10-5	3,25×10-8	3,34×10-10	3,25×10-14
1,4	6,83×10-5	6,87×10-5	1,55×10-5	2,99×10-7	2,24×10-10	6,22×10-13
1,6	1,43×10-4	1,44×10-4	1,67×10-5	3,05×10-7	4,41×10-9	6,98×10-13
1,8	2,57×10-4	2,58×10-4	1,23×10-4	3,28×10-7	7,89×10-9	7,23×10-13
2,0	4,05×10-4	4,07×10-4	3,38×10-4	7,99×10-7	9,33×10-9	8,66×10-13

Nesta tabela, como esperado, na medida em que se aumenta o valor do NPC , o erro absoluto em cada ponto de discretização t diminui. Em comparação com os resultados reportados por Lin e Liu ¹³ e por Li e Zeng ¹² , observa-se que, para NPC igual a 5, os erros tem a mesma ordem de grandeza e para NPC igual a 20, os erros obtidos são bem menores. Todavia, como o NPC significa, na prática, a dimensão do sistema algébrico a ser resolvido, NPC igual a 5 resulta em um sistema com 5 equações algébricas, o que implica em um sistema com menor dimensão do que os métodos propostos por Lin e Liu ¹³ e por Li e Zeng ¹² . Ao se considerar NPC igual a 10, a dimensão do sistema a ser resolvido pelo MCO é idêntico as abordagens M p -H e MAO, todavia, com um erro absoluto bem menor, i.e., da ordem de 10⁻⁸. É importante ressaltar que os resultados obtidos considerando as abordagens propostas por Lin e Liu ¹³ e por Li e Zeng ¹² podem ser melhorados na medida em que o número de pontos de discretização e/ou a ordem dos métodos propostos por este autores são aumentados. Em termos do tempo de processamento, em média, a metodologia proposta para a integração do modelo diferencial fracionário encontrou a solução em, aproximadamente, 0,12 segundos para o referido estudo de caso.

Já o último estudo de caso matemático foi estudado por Li e Zeng ¹² , sendo o mesmo descrito como segue:

$D^{\alpha }y\left(t\right)=-y\left(t\right)+\frac{t^{4-\alpha }}{\Gamma (5-\alpha )},0\le \alpha \le 1,0<t\le 1$ (6.3)

Para este estudo de caso, a condição inicial associada a esta equação diferencial fracionária é dada por $y\left(0\right)=0$ . Conforme apresentado por Li e Zeng ¹² , a solução analítica para este estudo de caso é dado por $y\left(t\right)=t^4{E}_{\alpha ,5}(-t^{\alpha })$ , onde E _α,β é a função de Mittag-Leffler, definida como sendo:

$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )},\alpha ,\beta >0$ (6.4)

em que z é a variável independente, α é a ordem fracionária e β é um parâmetro necessário para a avaliação desta função.

A tabela 3 apresenta os erros absolutos obtidos com a aplicação do MCO estendido para o contexto fracionário e aqueles obtidos pelo Método de Alta Ordem (MAO) ¹² no tempo t igual a unidade. Em ambas as estratégias, diferentes valores para o tamanho do passo de integração h (e ordem igual a 3) e diferentes valores para o NPC são considerados. Neste cenário, quanto menor o valor do parâmetro h , maior é a dimensão do modelo algébrico a ser resolvido.

Tabela 3:
Erros absolutos obtidos para a segunda EDOF considerando $t=1$ e diferentes valores para h e NPC .

Método	$\alpha =0,1$	$\alpha =0,3$	$\alpha =0,5$	$\alpha =0,7$	$\alpha =0,9$
$h=1/10$ 12	2,43×10-5	1,91×10-5	3,49×10-5	5,38×10-5	7,72×10-5
$h=1/20$ 12	3,19×10-6	2,50×10-6	4,52×10-6	6,88×10-6	9,76×10-6
$h=1/40$ 12	4,05×10-7	3,20×10-7	5,75×10-7	8,69×10-7	1,22×10-6
$h=1/80$ 12	5,06×10-8	4,05×10-8	7,24×10-8	1,09×10-7	1,53×10-7
$h=1/160$ 12	6,25×10-9	5,09×10-9	9,09×10-9	1,36×10-8	1,91×10-8
$h=1/320$ 12	7,71×10-10	6,38×10-10	1,13×10-9	1,71×10-9	2,39×10-9
$h=1/640$ 12	9,50×10-11	7,99×10-11	1,42×10-10	2,14×10-10	3,00×10-10
$NPC=5$	1,33×10-5	1,45×10-5	2,65×10-5	1,45×10-5	1,02×10-5
$NPC=10$	2,01×10-8	1,15×10-8	5,66×10-8	2,88×10-8	1,34×10-8
$NPC=15$	2,11×10-10	1,77×10-10	1,09×10-9	3,44×10-10	1,76×10-10
$NPC=20$	5,66×10-14	3,14×10-14	2,15×10-14	3,33×10-14	4,45×10-14

Conforme observado para o primeiro estudo de caso, na medida em que se aumenta o valor do NPC no MCO ou se reduz o valor de h no MAO ¹² , o erro absoluto no ponto t igual a unidade é reduzido. Em comparação com os resultados reportados por Li e Zeng ¹² observa-se que, para NPC igual a 5, os erros tem a mesma ordem de grandeza daqueles obtidos para h igual a 10. Todavia, como mencionado anteriormente, a dimensão do sistema considerando o MCO é menor do que aquele considerado pelo MAO. Ao se considerar a mesma dimensão do sistema a ser resolvido como parâmetro de comparação observa-se que, para NPC iguais a 10 e 20, os erros obtidos pelo MCO (da ordem de 10⁻⁸ e 10⁻¹⁴, respectivamente) são bem menores do que os reportados por Li e Zeng ¹² (da ordem de 10⁻⁵ e 10⁻⁶, para $h=1/10$ e $h=1/20$ , respectivamente). Assim, considera-se que os resultados obtidos pela metodologia proposta são melhores do que aqueles reportados por Li e Zeng ¹² . Finalmente, em termos do tempo de processamento, em média, o MCO encontrou a solução em, aproximadamente, 0,15 segundos.

6.2 Problema Inverso

6.2.1 Problema Inverso Matemático

Para avaliar a metodologia proposta para a resolução de problemas inversos formulados por equações diferenciais ordinárias fracionárias, esta seção apresenta os resultados obtidos para um estudo de caso puramente matemático onde foram gerados pontos experimentais sintéticos considerando desvio padrão ξ igual a zero (sem ruído) e ξ igual a 5% (com ruído), conforme a Eq. (5.3). Conforme destacado anteriormente, o objetivo é determinar a ordem fracionária α de modo a minimizar o funcional FO ₁, dado pela Eq. (5.2). A tabela 4 apresenta os valores médios para a ordem fracionária e para a função objetivo considerando o problema matemático estudado por Li e Zeng ¹² e descrito pela Eq. (6.3). De forma geral, observa-se nesta tabela que ambos os algoritmos (ED e BEF) foram capazes de estimar o valor de α para o problema sem ruído e obtiveram boas estimativas para o problema com ruído, para um mesmo número de avaliações da função objetivo, i.e., $50+50\times 500$ avaliações. Além disso, como esperado, o aumento no valor de ξ implica no aumento da dispersão dos pontos experimentais sintéticos e, consequentemente, no incremento do valor da função objetivo. Em termos do tempo de processamento, em média, a metodologia proposta encontrou a solução do problema inverso matemático em, aproximadamente, 85 segundos, contra 79 segundos do algoritmo de ED. Finalmente, cabe ressaltar que os resultados apresentados nesta tabela apresentaram um desvio padrão, em relação as 20 execuções dos algoritmos de otimização, da ordem de 10⁻⁶ para ambos os algoritmos.

Tabela 4:
Função objetivo estimada em um problema puramente matemático considerando $NPC=10$ .

		$\alpha =0,1$		$\alpha =0,5$		$\alpha =0,9$
		$\xi =0\%$	$\xi =5\%$	$\xi =0\%$	$\xi =5\%$	$\xi =0\%$	$\xi =5\%$
ED	α	0,0999	0,0801	0,5000	0,4781	0,9000	0,8245
	FO 1	9,2344x10-7	0,24434	2,9442x10-7	0,42323	2,4656x10-7	0,86543
BFE1	a	0,1000	0,0812	0,4999	0,4744	0,9000	0,8222
	FO 1	3,9834x10-8	0,23444	6,8877x10-8	0,20984	2,8887x10-8	0,71232
BFE2	a	0,1000	0,0812	0,5000	0,4786	0,9000	0,8244
	FO 1	1,4098x10-8	0,17431	3,9878x10-8	0,41554	1,9332x10-8	0,61112

6.2.2 Produção da Lacase

A tabela 5 apresenta os melhores valores obtidos com a execução de cada estratégia apresentada na seção 5. Nesta tabela também considera-se, para cada abordagem numérica, a ordem α constante e inteira.

Tabela 5:
Parâmetros estimados considerando diferentes metodologias numéricas e os modelos inteiro e fracionário.

	ED		BFE1		BFE2
α	1	0,9225	1	0,9227	1	0,9226
µmax (dia-α )	2,1860	2,3065	2,1857	2,1877	2,1865	2,1966
Ks (g/L)	7,6441	7,6156	7,6265	7,6245	7,6243	7,5960
YX/S (g/g)	0,6321	0,6410	0,6333	0,6408	0,6324	0,6304
K1(Ug − 1dia-α )	3743,4397	3915,1765	3743,4933	3914,3694	3743,4744	3915,5343
K2(Ug − 1dia-α )	7046,0178	7292,5605	7046,0123	7293,7587	7046,0822	7291,5644
Kp (dia-α )	0,9999	0,9341	0,9999	0,9259	0,9999	0,9349
ts (dia)	7,8050	7,7149	7,8065	7,7944	7,8054	7,7832
FO 2	0,6172	0,5896	0,6173	0,5863	0,6173	0,5898

De forma geral, observa-se que, para um mesmo valor de avaliações da função objetivo $(50+50\times 500)$ , os resultados obtidos pelo algoritmo de BFE estão em concordância com aqueles obtidos pelo algoritmo de ED, independentemente da ordem constante ou como variável de projeto. Em relação ao desvio padrão obtido com as execuções realizadas para cada um dos algoritmos, foi observado, no mínimo, um valor da ordem de 10⁻⁵ para o vetor de variáveis de projeto e para a função objetivo. Neste caso, constata-se que ambos os algoritmos convergiram para a mesma solução em cada uma das execuções realizadas. Em termos do tempo de processamento, em média, o algoritmo de BEF convergiu em, aproximadamente, 112 segundos, contra 97 segundos do algoritmo de ED. Ao se considerar a ordem como sendo uma variável de projeto, observa-se um valor de função objetivo melhor do que quando essa é mantida constante e igual a unidade. Isto se deve ao fato do aumento do graus de liberdade do sistema, i.e., ao se aumentar o número de variáveis de projeto (com a inclusão da ordem), flexibiliza-se a acomodação do modelos aos pontos experimentais, mesmo que essa não seja tão discrepante. Em termos da configuração do algoritmo de BFE, observa-se, para esta aplicação, um melhor desempenho da estratégia BFE¹, a qual faz uso do caminho gaussiano descrito pela Eq.(4.4). É importante ressaltar que, em termos da robustez de cada algoritmo, ambos são considerados satisfatórios visto o valor do desvio padrão obtido. As figuras 1(a)-(c) apresentam os perfis de concentração de biomassa, substrato e produto experimentais e preditas pelo modelo considerando as variáveis de projeto encontradas pela configuração BFE¹. Matematicamente, observa-se nestas figuras uma boa concordância entre os pontos experimentais e computados pelo modelo ajustado, principalmente em relação aos perfis de X e de S . Fisicamente, na figura 1(b) , a concentração de substrato diminui ao longo do processo, indicando que esta é consumida pela biomassa e, consequentemente, esta aumenta com o decorrer do processo, conforme a figura 1(a) . Já no perfil de produto P aumenta devido ao crescimento da concentração de biomassa até chegar ao tempo t _s , a partir do qual observa-se a sua redução ao longo do processo, conforme observado na figura 1(c) . É importante ressaltar que, como o modelo cinético adotado para o balanço de massa de produto (Eq. (2.3)) é definido por expressões distintas e que dependem do parâmetro t _s , observa-se na figura 1(c) a descontinuidade na derivada para o perfil de produto. Fisicamente, não existe razão aparente para o aparecimento desta descontinuidade no processo analisado. Por outro lado, do ponto de vista matemático, esta pode ser eliminada a partir do uso de um único modelo cinético para todo o intervalo de tempo ou a partir do atendimento da condição de diferenciabilidade no tempo ts para ambas as expressões definidas na Eq. (2.3).

Figura 1:
Perfis de concentração de biomassa, substrato e produto experimentais e preditos.

6.3 Influência da Ordem Fracionária nos Perfis Simulados

Para avaliar a influência da ordem fracionária nos perfis simulados para esse processo, nas figuras 2(a)-(c) são apresentadas as concentrações de biomassa, substrato e produto considerando os seguintes parâmetros: ${\mu }_{max}=2,1860{\mathrm{dia}}^{-\alpha }$ ; $K_s=7,6441\mathrm{g}/\mathrm{L}$ ; $Y_{X/S}=0,6321\mathrm{g}/\mathrm{g}$ ; $K_1=3743,4397{\mathrm{Ug}}^{-1}{\mathrm{dia}}^{-\alpha }$ ; $K_2=7046,0178{\mathrm{Ug}}^{-1}{\mathrm{dia}}^{-\alpha }$ ; $K_p=0,9999{\mathrm{dia}}^{-\alpha }$ ; $t_s=7,8050$ dia e diferentes valores para a ordem fracionária $\alpha \left(\right[0.40.5...0.911.1...1.92\left]\right)$ .

Figura 2:
Perfis de concentração de biomassa, substrato e produto em função da ordem fracionária α.

Nestas figuras é importante observar que o aumento no valor do parâmetro α implica em um comportamento oscilatório, i.e., se o valor de α de aproxima de 2, o modelo diferencial se torna um com ordem 2 e, naturalmente, este tende a apresentar um comportamento oscilatório, como por exemplo, de um oscilador harmônico clássico. Além disso, também é possível observar, para algumas variáveis dependentes e para alguns valores de α , perfis negativos, o que implica em um comportamento fisicamente inviável para este processo. Neste caso, ao se propor um problema inverso onde a ordem é uma variável de projeto deve-se observar os perfis obtidos para avaliar, do ponto de vista físico, se os mesmo são coerentes ou não, visto que podem ser obtidos perfis inviáveis do ponto de vista físico.

7 CONCLUSÕES

O presente trabalho teve por objetivo associar Método da Colocação Ortogonal com o algoritmo de Busca Fractal Estocástico para a resolução de problemas inversos modelados por equações diferenciais fracionárias. Para essa finalidade, o Método da Colocação Ortogonal, proposto originalmente para a resolução de equações diferenciais com ordem inteira, foi estendido para resolver de uma equação diferencial ordinária com ordem fracionária no tempo. Para avaliar a metodologia proposta para a resolução de equações diferenciais ordinárias fracionárias, dois estudos de caso puramente matemáticos, cujas soluções analíticas são conhecidas, foram apresentados. Em linhas gerais, foi possível observar que a metodologia proposta foi capaz de obter resultados satisfatórios em relação aos reportados considerando outros métodos numéricos, conforme observado pelo valor do erro encontrado. Além disso, como esperado, o aumento do número de pontos de colocação implica no aumento da dimensão do sistema a ser resolvido, mas também implica no aumento da precisão do método. Assim, pode-se avaliar o custo benefício para a escolha entre o número de pontos de colocação e a dimensão do problema a ser resolvido. Ressalta-se que a metodologia proposta pode ser estendida para problemas formulados como equações diferenciais parciais fracionárias no tempo e/ou no espaço.

Com relação aos problemas inversos propostos observa-se que o algoritmo de Busca Fractal Estocástico foi capaz de obter resultados satisfatórios em relação ao já consolidado algoritmo de Evolução Diferencial, i.e., foram observados bons resultados em termos do valor da função objetivo para um mesmo número de avaliações da função objetivo em ambos os algoritmos. No que tange o problema inverso usando pontos experimentais reais, embora não tenha sido observada uma grande discrepância entre os valores referentes a função objetivo para α igual a unidade e para α como variável de projeto, os resultados obtidos foram promissores. Do ponto de vista físico, o valor de α computado está próximo ao que a literatura tem utilizado durante a etapa de modelagem deste tipo de processo. Do ponto de vista matemático, o aumento do número de variáveis de projeto implica no aumento do número de graus de liberdade do processo e, a priori, no aumento da chance de uma maior acomodação dos pontos experimentais ao modelo.

Cabe ressaltar que o estudo de modelos fracionários configura uma linha de pesquisa que visa generalizar a descrição dada pela respectiva de uma equação diferencial de ordem inteira. Do ponto de vista físico, a generalização destes modelos pode implicar na obtenção de representações matemáticas mais realísticas. Neste cenário, pode-se: i ) avaliar a influência da ordem fracionária com relação aos perfis físicos; e ii ) analisar aspectos relacionados com a análise dimensional destes modelos.

Como sugestões para trabalhos futuros pode-se citar a formulação e a resolução de problemas inversos relacionados com a identificação de modelos difusivos anômalos. Além disso, avaliar a influência de outros tipos de aproximação para as derivadas fracionárias no que tange a simulação dos processos e a resolução de problemas inversos.

Agradecimentos

O autor agradece o apoio financeiro da FAPEMIG, CAPES e do CNPq para a realização deste trabalho.

REFERÊNCIAS

• ¹
  F. Alves. “Modelagem e Simulação de Biorreator Operando com Fungos Trametes versicolor para Produção de Enzima Lacase”. Ph.D. thesis, Escola de Engenharia Mauá do Centro Universitário do Instituto Mauá de Tecnologia, São Caetano do Sul, Dissertação de Mestrado (2010).
• ²
  R.L. Bagley & R.A. Calico. Fractional order state equations for the control of viscoelastic strucures. J. Guid. Control Dyn., 2(14) (1999), 1-15.
• ³
  D.A. Benson, S.W. Wheatcraft & M.M. Meerschaert. Application of a fractional advection-despersion equation. Water Resour. Res., 6(36) (2000), 1403-1412.
• ⁴
  X. Chen, H. Yue & K. Yu. Perturbed stochastic fractal search for solar PV parameter estimation. Energy, 3(https://doi.org/10.1016/j.energy.2019.116247) (2020), 1-14.
  » https://doi.org/https://doi.org/10.1016/j.energy.2019.116247
• ⁵
  G. Chi, G. Li & X. Jia. Numerical inversions of a source term in the FADE with a Dirichlet boundary condition using final observations. Computers and Mathematics with Applications, 6(62) (2011), 1619-1626.
• ⁶
  M. Dalir & M. Bashour. Applications of fractional calculus. Applied Mathematical Sciences, 21(4) (2010), 1021-1032.
• ⁷
  W. Deng. Short memory principle and a predictor-corrector approach for fractional differential equations. J. Comput. Appl. Math, 3(206) (2007), 174-188.
• ⁸
  R. Herrmann. “Fractional calculus: An introduction for physicists”. World Scientific (2011).
• ⁹
  M. Ichise, Y. Nagayanagi & T. Kojima. An analog simulation of non-integer order transfer functions for analysis of electrode processes. J. Electroanal. Chem. Interfacial Electrochem, 1(33) (1971), 253- 265.
• ¹⁰
  D.P. Leitoles. “Análise de Sistemas Reacionais e de Separação usando Cálculo Fracionário”. Ph.D. thesis, EESC, Universidade de São Paulo, Universidade Federal do Paraná (2015).
• ¹¹
  C. Li, A. Chen & J. Ye. Numerical approaches to fractional calculus and fractional ordinary differential equation. Journal of Computational Physics, 2(230) (2011), 3352-3368.
• ¹²
  C. Li & F. Zeng. The Finite Difference Methods for Fractional Ordinary Differential Equations. Numerical Functional Analysis and Optimization, 2(34) (2013), 149-179.
• ¹³
  R. Lin & F. Liu. Fractional high order methods for the nonlinear fractional ordinary differential equation. Nonlinear Analysis, 1(66) (2007), 856-869.
• ¹⁴
  R.L. Magin. “Fractional calculus in bioengineering”. Begell House Redding (2006).
• ¹⁵
  B. Mandelbrot. Some noises with 1/f spectrum, a bridgre between direct current and white noise. IEEE Trans. Inform. Theory, 2(13) (1967), 289-298.
• ¹⁶
  B.B. Mandelbrot . “The fractal geometry of nature”. W. H. Freeman and Company, New York (2003).
• ¹⁷
  R.J. Mark & M.W. Hall. Differentegral interpolation from a bandlimited signal’s samples. IEEE Trans.Acoust. Speech Signal Process, 1(29) (1981), 872-877.
• ¹⁸
  K.B. Oldham & J. Spanier. “The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order”. Academic Press, New York (1974).
• ¹⁹
  I. Podlubny. “Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications”. Academic Press, San Diego, CA, USA (1999).
• ²⁰
  H. Salimi. Stochastic Fractal Search: A powerful metaheuristic algorithm. Knowledge-Based Systems, 3(75) (2015), 1-18.
• ²¹
  R. Storn & K. Price. Differential Evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical Report TR-95-012, International Computer Science Institute, Berkeley, 2(2) (1995), 1-17.
• ²²
  H.H. Sun, N. Onaral & Y. Tsao. Application of position reality principle to metal electrode linear polarization phenomena. IEEE Trans. Ciomed. Eng., 10(31) (1984), 664-674.
• ²³
  J. Villadsen & M.L. Michelsen. “Solution of differential equation models by polynomial approximation”. Prentice-Hall, Englewood Cliffs (1978).
• ²⁴
  J. Villadsen & W.E. Stewart. Solution of boundary-value problems by orthogonal collocation. Chemical Engineering Science, 3(22) (1967), 1483-1501.
• ²⁵
  X. Xiong, H. Guo & X. Liu. An inverse problem for a fractional diffusion equation. Journal of Computational and Applied Mathematics, 4(236) (2012), 4474-4484.
• ²⁶
  L. Yuan & O.P. Agrawal. A numerical scheme for dynamic systems containing fractional derivatives. ASME J. Vibr. Acoust., 2(124) (2002), 321-324.
• ²⁷
  G.H. Zheng & T. Wei. A new regularization method for solving a time-fractional inverse diffusion problem. J. Math. Anal. Appl., 2(378) (2011), 418-431.

Datas de Publicação

• Publicação nesta coleção
  30 Nov 2020
• Data do Fascículo
  Sep-Dec 2020

Histórico

• Recebido
  14 Abr 2020
• Aceito
  01 Set 2020