UNIVERSIDADE ESTADUAL PAULISTA "JULIO DE MESQUITA FILHO" FACULDADE DE CIÊNCIAS FARMACÊUTICAS CÂMPUS DE ARARAQUARA PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA DE BIOMATERIAIS E BIOPROCESSOS MESTRADO PROFISSIONAL OTIMIZAÇÃO DE BIOPROCESSOS BASEADA EM MODELOS MATEMÁTICOS E CÁLCULO VARIACIONAL ANGEL GUSTAVO TOLABA ORIENTADOR: Prof. Dr. Samuel Conceição de Oliveira ARARAQUARA - SP 2019 UNIVERSIDADE ESTADUAL PAULISTA "JULIO DE MESQUITA FILHO" FACULDADE DE CIÊNCIAS FARMACÊUTICAS CÂMPUS DE ARARAQUARA OTIMIZAÇÃO DE BIOPROCESSOS BASEADA EM MODELOS MATEMÁTICOS E CÁLCULO VARIACIONAL ANGEL GUSTAVO TOLABA Dissertação apresentada ao Programa de Pós-Graduação em Engenharia de Biomateriais e Bioprocessos (Mestrado Profissional), Área de Biomateriais, Bioprocessos, Bioprodutos da Faculdade de Ciências Farmacêuticas, UNESP, como parte dos requisitos para obtenção do Título de Mestre em Engenharia de Biomateriais e Bioprocessos. ORIENTADOR: Prof. Dr. Samuel Conceição de Oliveira ARARAQUARA - SP 2019 Tolaba, Angel Gustavo. T647o Otimização de bioprocessos baseada em modelos matemáticos e cálculo variacional / Angel Gustavo Tolaba. – Araraquara: [S.n.], 2019. 66 f. : il. Dissertação (Mestrado Profissional) – Universidade Estadual Paulista. “Júlio de Mesquita Filho”. Faculdade de Ciências Farmacêuticas. Programa de Pós Graduação em Engenharia de Biomateriais e Bioprocessos. Área de Biomateriais, Bioprocessos, Bioprodutos. Orientador: Samuel Conceição de Oliveira. 1. Modelagem matemática. 2. Simulação. 3. Otimização. 4. Controle. 5. Bioprocessos. 6. Cálculo variacional. I. Oliveira, Samuel Conceição de, orient. II. Título. Diretoria do Serviço Técnico de Biblioteca e Documentação - Faculdade de Ciências Farmacêuticas UNESP - Campus de Araraquara CAPES: 33004030170P0 Esta ficha não pode ser modificada RESUMO A determinação da estratégia de controle adequada para bioprocessos batelada e batelada alimentada é uma questão prática importante devido ao alto valor agregado de alguns bioprodutos. Desde que é altamente desejável otimizar a produção de bioprodutos, vários métodos têm sido propostos para esse objetivo. Uma vez dispondo de um modelo matemático adequado para o bioprocesso, o problema de otimização pode ser formulado no âmbito do princípio do máximo de Pontryagin e da teoria de controle ótimo para determinar a melhor trajetória de controle para certas variáveis manipuladas, como temperatura, pH e taxa de alimentação do substrato. Neste estudo, duas aplicações dessas técnicas baseadas em modelos matemáticos para otimizar e controlar bioprocessos de produção de antibióticos são revisadas e novos aspectos são enfatizados. Os casos analisados incluem a otimização da taxa de alimentação de substrato em um reator batelada alimentada e da temperatura em um reator batelada durante fermentações penicilínicas. Os principais resultados obtidos neste estudo foram: (i) a constatação de que métodos numéricos simples (Runge-Kutta, Newton-Raphson) podem ser aplicados para resolver satisfatoriamente os problemas de valor no contorno propostos; (ii) a demonstração de que a operação não isotérmica é mais produtiva em antibiótico do que a operação sob temperatura constante; (iii) a necessidade de adoção de um modelo matemático apropriado para o bioprocesso visando à resolução do problema de controle ótimo. O principal objetivo deste estudo é fornecer uma base teórica para a aplicação de uma metodologia baseada em modelos matemáticos que possa ser usada para a otimização e controle de outros bioprocessos com diferentes características. Palavras-chave: Modelagem matemática. Simulação. Otimização. Controle. Bioprocessos. Cálculo variacional. 1 ABSTRACT The control strategy determination for batch and fed-batch bioprocesses is an important practical issue due to the high added value of some bioproducts. Since it is highly desirable to optimize the bioproduct production, several methods have been proposed aimed at this objective. Once having a mathematical model for the bioprocess, the optimization problem can be formulated within the framework of Pontryagin's maximum principle and of the optimal control theory to determinate the best control trajectory for certain key manipulated variables, such as temperature, pH, and substrate feed rate. In this study, two applications of these mathematical- model-based techniques to optimize and control antibiotics production bioprocesses are reviewed and new aspects are emphasized. The cases analyzed included the optimization of the substrate feed rate in a fed-batch reactor and of the temperature in a batch reactor during penicillin fermentations. The main results obtained in this study were: (i) the realization that simple numerical methods (Runge-Kutta, Newton- Raphson) can be applied to satisfactorily solve the proposed boundary value problems, (ii) the demonstration that the non-isothermal operation is more productive in antibiotic than the operation under constant temperature, and (iii) the need to adopt an appropriate mathematical model for bioprocess to solve the optimal control problem. The main aim of this study is to provide a theoretical basis for the application of a mathematical-model-based methodology that can be used for the optimization and control of other bioprocesses with different characteristics. Keywords: Mathematical modeling. Simulation. Optimization. Control. Bioprocesses. Variational calculus. 2 LISTA DE FIGURAS Figura 1. Perfis das variáves de estado S, X e P durante o bioprocesso de produção de penicilina em modo batelada alimentada ............................................................. 28 Figura 2. Perfil da variável de estado V e da variável de controle u durante o bioprocesso de produção de penicilina em modo batelada alimentada .................... 29 Figura 3. Perfil da variável adjunta 1 durante uma fermentação penicilínica não isotérmica .................................................................................................................. 32 Figura 4. Perfil de concentração adimensional de células durante uma fermentação penicilínica não isotérmica ........................................................................................ 33 Figura 5. Perfil de concentração adimensional de produto durante uma fermentação penicilínica não isotérmica ........................................................................................ 33 Figura 6. Perfil de temperatura durante uma fermentação penicilínica não isotérmica .................................................................................................................................. 34 Figura 7. Perfis de concentração de células e de produto durante fermentações penicilínicas isotérmicas a diferentes temperaturas .................................................. 36 Figura 8. Perfil de concentração adimensional de células durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular ........................................................................................................................ 38 Figura 9. Perfil de concentração adimensional de produto durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular ........................................................................................................................ 39 Figura 10. Perfil de temperatura durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular ............... 39 Figura 11. Perfil da variável adjunta 1 durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular ............... 40 3 LISTA DE TABELAS Tabela 1. Parâmetros cinéticos e variáveis operacionais utilizadas no estudo de caso da produção de penicilina em modo batelada alimentada (Fonte: Costa, 1996) ....... 27 Tabela 2. Dados da integração numérica das equações de balanço de massa utilizados para determinação do 10 tempo de troca durante a produção de penicilina em modo batelada alimentada .................................................................................. 28 Tabela 3. Valores finais da concentração de células (y1) e de produto (y2) em fermentações penicilínicas isotérmicas e não isotérmicas ........................................ 35 4 SUMÁRIO 1 INTRODUÇÃO ..................................................................................................... 6 2 OBJETIVOS ......................................................................................................... 8 3 REVISÃO BIBLIOGRÁFICA ................................................................................. 9 3.1 Bioprocesso ................................................................................................ 9 3.2 Modelagem de bioprocessos ................................................................... 12 3.3 Otimização de Processos ......................................................................... 13 3.4 Controle Ótimo e Aplicação do Princípio do Máximo ............................ 14 4 METODOLOGIA ................................................................................................ 17 4.1 Estudo dos fundamentos teóricos do Princípio do Máximo de Pontryagin ............................................................................................................ 17 4.2 Análise dos modelos matemáticos referentes aos bioprocessos escolhidos para estudo ....................................................................................... 17 4.3 Estudos de casos e simulação dos bioprocessos nas condições otimizadas ............................................................................................................ 17 5 RESULTADOS E DISCUSSÃO ......................................................................... 18 5.1 Estudo dos fundamentos teóricos do Princípio do Máximo de Pontryagin ............................................................................................................ 18 5.2 Análise dos modelos matemáticos referentes aos bioprocessos escolhidos para estudo ....................................................................................... 18 5.2.1 Modelo de Constantinides e Mostoufi (1999) ............................................ 18 5.2.2 Modelo de Costa et al. (1996) ................................................................... 19 5.3 Otimização de bioprocessos usando o Princípio do Máximo de Pontryagin ............................................................................................................ 20 5.3.1 Determinação do perfil ótimo de vazão específica de alimentação de substrato ............................................................................................................ 20 5.3.1.1 Simulação do bioprocesso de produção de penicilina nas condições otimizadas – Caso 1 ....................................................................................... 28 5.3.2 Determinação do perfil ótimo de temperatura durante a produção em batelada de penicilina ........................................................................................ 29 5.3.2.1 Simulação do bioprocesso de produção de penicilina nas condições otimizadas – Caso 2 ....................................................................................... 31 5 5.3.2.2 Simulação do bioprocesso de produção de penicilina em outras condições diferentes da ótima – Caso 2 ......................................................... 35 5.3.2.3 Modelo matemático com restrição – correção da curva de crescimento ....................................................................................................................... 36 6 CONCLUSÕES .................................................................................................. 41 7 REFERÊNCIAS .................................................................................................. 42 APÊNDICE ................................................................................................................ 47 6 1 INTRODUÇÃO A melhoria na produtividade de muitos bioprocessos é realizada pela manipulação de parâmetros nutricionais e físicos, como composição do meio de cultivo, velocidade de agitação, taxa de aeração, pH e temperatura (SIRCAR et al., 1998; CHEN, 1998; OLIVEIRA, 2018). Embora a determinação das condições ótimas para um bioprocesso multivariável seja muitas vezes tediosa, é possível aplicar métodos racionais para este fim tais como os delineamentos experimentais (CHEN, 1998; OLIVEIRA, 2018). Delineamentos experimentais podem ser divididos em dois grupos (BERKHOLZ; GUTHKE, 2002; KENNEDY; KROUSE, 1999; OLIVEIRA, 2018): (i) delineamentos experimentais baseados em modelos e (ii) delineamentos experimentais estatísticos. Em delineamentos experimentais baseados em modelo, as previsões de um modelo matemático são usadas para determinar como um experimento ou processo deve ser realizado ou conduzido, enquanto que em delineamentos experimentais estatísticos, essas previsões de modelo não são explicitamente requeridas (OLIVEIRA, 2018). A otimização de bioprocessos desempenha um papel fundamental no setor de biotecnologia devido à forte concorrência entre as empresas (OLIVEIRA, 2018). Metabólitos secundários, como antibióticos e outros produtos farmacêuticos, possuem alto valor agregado e, assim, melhorias na produção desses bioprodutos são de grande interesse industrial (OLIVEIRA, 2018). Para alcançar operações de alto desempenho, a otimização das variáveis manipuladas que afetam o bioprocesso torna-se uma tarefa indispensável (OLIVEIRA, 2018). Em geral, os problemas de otimização podem ser classificados em duas categorias: otimização do ponto de operação e otimização de perfis (LIM; LEE, 1982; OLIVEIRA, 2018). O primeiro caso consiste em determinar o melhor conjunto de valores das variáveis manipuladas que levam à maximização do índice de desempenho. O segundo caso consiste em determinar funções temporais ou espaciais (perfis), em vez de um ponto no espaço n-dimensional, que fornecem um valor ótimo para o índice de desempenho (LIM; LEE, 1982; OLIVEIRA, 2018). Como exemplo do segundo tipo de otimização, sabe-se que nos bioprocessos de produção de antibióticos, as condições de temperatura e pH para a máxima produção de produto são diferentes daquelas para o máximo crescimento celular, 7 requerendo, deste modo, a implementação de perfis de temperatura e pH para um melhor desempenho do bioprocesso (OLIVEIRA, 2018). Considerando que o principal objetivo de um bioprocesso é a produção rentável de bioprodutos, é importante selecionar o modo de operação mais apropriado que permita a obtenção de uma alta concentração, produtividade e rendimento no produto desejado (LEE et al.,1999; OLIVEIRA, 2018). Nesse contexto, os processos em batelada alimentada (fed batch) têm sido amplamente utilizados para a produção de vários bioprodutos, incluindo metabólitos primários e secundários. No caso particular de metabólitos secundários, tais como os antibióticos, a interação entre o metabolismo de crescimento e a biossíntese de produto é fortemente influenciada pelas concentrações dos nutrientes limitantes do crescimento, isto é, nesses processos, uma alimentação insuficiente ou excessiva de nutrientes prejudica o crescimento celular e a formação de produto devido à possibilidade de ocorrência de fenômenos tais como inanição celular e repressão catabólica, o que requer uma estratégia de alimentação adequada para esses bioprocessos. Em sendo requerida uma particular sequência de valores no tempo para as variáveis de controle visando conduzir o bioprocesso ao longo de uma trajetória que proporcione a maior produtividade, é possível que, em alguns casos, esses perfis temporais sejam muito complexos a ponto de inviabilizar a determinação dos mesmos apenas experimentalmente. Assim, métodos matemáticos e numéricos podem ser aplicados para a determinação desses perfis, reduzindo o esforço experimental e o tempo necessários para a otimização (OLIVEIRA, 2018). A determinação dos perfis ótimos de pH, temperatura e taxa de alimentação de substrato em bioprocessos batelada e batelada alimentada é um típico problema de otimização cuja solução requer o uso de modelos cinéticos e técnicas matemáticas poderosas. Várias abordagens para a determinação de perfis temporais ótimos para variáveis de controle têm sido relatados na literatura. Nestes relatos, o problema de otimização é geralmente formulado com base no Princípio do Máximo de Pontryagin (PMP), tomando como ponto de partida um modelo matemático fenomenológico do bioprocesso. Para modelos matemáticos simples, o problema pode ser resolvido analiticamente, a partir do Hamiltoniano do sistema, aplicando um procedimento iterativo na variável de controle para determinar o seu perfil ótimo (OLIVEIRA, 2018). 8 2 OBJETIVOS O objetivo geral desta dissertação é revisitar aplicações do Princípio do Máximo de Pontryagin à otimização de bioprocessos visando à compreensão e ao domínio completo do método bem como buscar e discutir aspectos não explorados originalmente nos trabalhos revisitados. Para isso, são propostos os seguintes objetivos específicos: ● estudar os fundamentos teóricos do Princípio do Máximo de Pontryagin; ● analisar os modelos matemáticos referentes aos bioprocessos escolhidos para estudo; ● realizar estudos de casos e simular os bioprocessos nas condições otimizadas. 9 3 REVISÃO BIBLIOGRÁFICA 3.1 Bioprocesso Um bioprocesso é um processo no qual são usadas células vivas ou alguns de seus componentes como agentes biológicos para transformação de matérias- primas em produtos de interesse econômico e/ou social tais como enzimas, ácidos e solventes orgânicos, antibióticos, vacinas e hormônios, fermento, bebidas, biopolímeros, bioplásticos, aminoácidos, biocombustíveis e outros (OLIVEIRA, 2018). Os bioprocessos podem ocorrer na presença ou ausência de oxigênio, originando as denominações clássicas dadas a estes processos de aeróbios e anaeróbios, respectivamente. Os bioprocessos anaeróbios são ainda denominados de fermentação ou processo fermentativo (OLIVEIRA, 2018). Os bioprocessos são conduzidos em equipamentos denominados de biorreatores, os quais são geralmente tanques cilíndricos agitados, aerados ou não e dotados de dispositivos de monitoramento e controle de pH, temperatura, oxigênio dissolvido, etc. (OLIVEIRA, 2018). Dependendo do modo como o biorreator é operado, quanto à adição e retirada de materiais, os bioprocessos são classificados em três tipos: batelada, contínuo e batelada alimentada. No primeiro caso, todos os materiais usados na biorreação são alimentados ao biorreator no início da operação e nenhuma adição ou retirada de materiais ocorre durante o bioprocesso. No segundo caso, um sistema aberto é configurado. Meio de cultivo não convertido é adicionado continuamente ao biorreator enquanto que uma igual quantidade de meio convertido, contendo microrganismos, é simultaneamente removida do sistema. No bioprocesso descontínuo alimentado, o substrato é adicionado segundo um perfil de alimentação pré-estabelecido (constante, linear, exponencial, etc.) até a capacidade volumétrica máxima do biorreator ser atingida (OLIVEIRA, 2018; BAYER et al., 2015; SBARCIOG et al., 2012; ZHANG et al., 2010; CRISTALDI et al., 2009). Os bioprocessos existem desde o início da civilização humana (PALANKI et al., 1993). Fazer pão e vinho são alguns dos processos fermentativos que os seres humanos sempre utilizaram para a sobrevivência e o prazer. Embora ligados fortemente à vida cotidiana do homem, os bioprocessos não receberam muita 10 atenção nas atividades de pesquisa em biotecnologia e bioengenharia até a segunda metade do século XX (CHEN et al.,2006). Uma aplicação importante e bem sucedida de bioprocessos na história é a produção de penicilina. Em 1941, apenas baixas concentrações de penicilina (1 mg/L) eram possíveis de ser obtidas por técnicas de cultura de superfície em escala de laboratório, mesmo quando cepas de alto rendimento eram usadas. A demanda por penicilina naquele momento excedia a capacidade de produção. Assim, durante a segunda guerra mundial, um enorme esforço foi feito pelos aliados para desenvolver tecnologia com vistas ao aumento de escala do bioprocesso de produção de penicilina a nível industrial, culminando com o advento da Engenharia de Bioprocessos e o estabelecimento do protocolo do desenvolvimento de processos industriais de fermentação. Posteriormente, na década de 70, utilizando- se técnicas de mutação induzida das linhagens produtoras, a produtividade em antibiótico foi significativamente aumentada para mais de 50 g/L (CHEN et al., 2006). Com o avanço das técnicas biotecnológicas, incluindo a manipulação genética de microrganismos, vários produtos inovadores, tais como aqueles de “química fina” (alcaloides, glicosídeos cardíacos, antitumorais, antimicrobianos, esteroides, corantes, adoçantes, aromatizantes, etc.) e biológicos (proteínas recombinantes, tecidos, vacinas, soros, etc.) foram obtidos a partir de bioprocessos, contribuindo significativamente para a melhoria da saúde e da qualidade de vida humana (BAYEN et al., 2015; SBARCIOG et al., 2012 ; ZHANG et al., 2010; CRISTALDI et al., 2009; CHEN et al., 2006). Contudo, os altos custos associados a muitos bioprocessos representam um gargalo para o desenvolvimento, produção e aplicação de certos produtos. Desenvolver um método de cultivo econômico e ambientalmente saudável é o principal objetivo atualmente pretendido no desenvolvimento de bioprocessos, isto é, o objetivo é controlar o bioprocesso em sua trajetória ótima de modo a atingir produtividade máxima com mínimo custo de produção e garantia da qualidade do produto. Um bioprocesso deve ser, em tese, sempre operado de forma otimizada por várias razões. Um bioprocesso otimizado e controlado proporciona a obtenção de produto com alta pureza, segurança operacional, atendimento às regulamentações ambientais e redução de custos. A otimização e controle de bioprocessos ainda constituem tarefas desafiadoras, principalmente devido às seguintes razões (CHEN et al., 2006): 11  a natureza inerente não linear e variável no tempo (dinâmica) torna o processo extremamente complexo;  modelos de processos precisos raramente estão disponíveis devido à complexidade dos processos bioquímicos subjacentes;  as respostas das variáveis do bioprocesso, em particular das concentrações celulares e metabólicas, são muito lentas e os parâmetros do modelo variam de maneira complexa;  sensores on-line confiáveis que possam medir com precisão as variáveis de estado importantes raramente estão disponíveis. É muito importante monitorar o bioprocesso para estimar seu estado visando a implementação de estratégias de controle. Por exemplo, oxigênio dissolvido (OD) e pH são as variáveis mais comumente medidas com sensores eletroquímicos. No entanto, algumas variáveis de estado chave, como a concentração de biomassa, não podem ser medidas diretamente devido à falta de sensores adequados ou ao custo elevado de tais sensores. Nos últimos anos, muitos esforços foram feitos no sentido de desenvolver sensores de software online (softsensor). O conceito chave de técnicas de softsensor é estimar estados não medidos a partir de estados medidos. Os estados não medidos são normalmente inacessíveis ou difíceis de medir por meio de biossensores ou sensores de hardware, enquanto os estados medidos são relativamente fáceis de monitorar on-line usando instrumentos confiáveis. Com base nesta filosofia, várias técnicas softsensor têm sido propostas na literatura, a saber (CHEN et al.,2006):  estimação usando balanços elementares;  observadores de estado adaptativos;  estimação usando filtros (filtro de Kalman e filtro de Kalman estendido). Os dois primeiros métodos apresentam dificuldades relacionadas às imprecisões dos instrumentos e modelos disponíveis. O terceiro método requer muito esforço computacional, estimativas anteriores do ruído da medição e características da incerteza do modelo. Também apresenta alguns problemas numéricos e dificuldades de convergência devido à aproximação associada à linearização de modelos (CHEN et al.,2006). 12 3.2 Modelagem de bioprocessos Por muitos anos, a dinâmica dos bioprocessos em geral foi modelada por um conjunto de equações diferenciais não lineares de primeira ou de maior ordem. Esses modelos matemáticos podem ser divididos em duas categorias diferentes: modelos estruturados e modelos não estruturados. Os modelos estruturados representam os processos ao nível subcelular, enquanto os modelos não estruturados representam os processos ao nível da população (extracelular), considerando as células como mais um soluto presente no meio (CHEN et al.,2006). Lei et al. (2001) propuseram um modelo bioquimicamente estruturado para leveduras cultivadas em um quimiostato, baseado em equações cinéticas do tipo de Monod. Um conjunto de dados experimentais obtidos em estado estacionário foi bem descrito pelo modelo. No entanto, quando o modelo foi aplicado a um cultivo em batelada alimentada, erros relativamente grandes foram observados entre os valores experimentais e aqueles preditos pelo modelo (CHEN et al.,2006). Outro modelo estruturado para simular o crescimento de levedura de panificação em biorreatores industriais foi apresentado por Di Serio et al. (2001). A modelagem detalhada dos processos de regulação foi substituída por uma estrutura de modelagem cibernética, baseada na hipótese de que os microrganismos otimizam a utilização dos substratos disponíveis para maximizar sua taxa de crescimento em todos os instantes. As previsões do modelo concordaram razoavelmente bem com dados experimentais obtidos em fermentações batelada em escala de laboratório e industrial. A limitação do modelo, como apontado pelos autores, era que o modelo e seus parâmetros precisavam ser melhorados para uma aplicação mais geral (CHEN et al., 2006). Um modelo não estruturado para fermentações industriais utilizando leveduras foi apresentado por Pertev et al. (1997). As cinéticas do metabolismo das leveduras consideradas para construção do modelo foram baseadas na hipótese de capacidade respiratória limitada, proposta por Sonnleitner e Kepelli (1986). O modelo foi testado para dois diferentes modos de fermentação industrial: batelada e batelada alimentada, sendo capaz de prever, com precisão suficiente, o comportamento do processo em ambas modalidades de fermentação. Posteriormente, um estudo de Berber et al. (1999) utilizando esse modelo, mostrou que um melhor perfil da taxa de alimentação do substrato poderia ser implementado 13 visando aumentar a produção de biomassa e diminuir a formação de etanol, uma vez que o objetivo do bioprocesso é produzir biomassa (CHEN et al., 2006). Como as variáveis de estado intracelulares (isto é, enzimas e outros componentes) não estão envolvidas em modelos não estruturados, é relativamente fácil validar esses tipos de modelos por meio de experimentos. É por isso que modelos não estruturados são mais preferíveis do que modelos estruturados para otimização e controle de bioprocessos (CHEN et al., 2006). Entretanto, modelos não estruturados também apresentam problemas de identificação de parâmetros e grandes erros de previsão em alguns casos. Os parâmetros de modelos matemáticos variam de uma cultura para outra. Os métodos convencionais de identificação de parâmetros do sistema, como Mínimos Quadrados, Mínimos Quadrados Recursivos, Máxima Verossimilhança ou Variável de Instrumento, funcionam bem para sistemas lineares. Esses esquemas, no entanto, são, em essência, técnicas de busca local e muitas vezes falham na busca pelo ótimo global se o espaço de busca não for diferenciável ou não for linear nos parâmetros. Embora um esforço considerável tenha sido feito no desenvolvimento de modelos matemáticos detalhados, os bioprocessos são muito complexos para serem completamente descritos matematicamente (CHEN et al., 2006). Do ponto de vista de aplicação, as limitações dos modelos matemáticos são: ● necessidade de um conhecimento preliminar sobre a fisiologia do microrganismo e a influência de variáveis físicas; ● número reduzido de metabólitos que podem ser incluídos na descrição; ● capacidade reduzida de descrever variações entre bateladas; ● aplicação restrita a condições ideais de fermentação; ● elevado número de equações diferenciais e parâmetros mesmo em descrições matemáticas moderadamente complexas. 3.3 Otimização de Processos O desenvolvimento sistemático de estratégias de controle ótimas para bioprocessos, especialmente aqueles em batelada alimentada, é de particular interesse não só acadêmico, mas também, tecnológico. 14 Muitos produtos obtidos por processos biotecnológicos, tais como produtos farmacêuticos, agrícolas, químicos e biológicos, são produzidos comercialmente utilizando-se fermentações em batelada alimentada. Geralmente, o rendimento final da fermentação em batelada alimentada é superior ao da fermentação em batelada porque um melhor controle do processo pode ser realizado. Entretanto, manter o equilíbrio entre a taxa de alimentação de substrato e a capacidade respiratória do microrganismo é uma tarefa crítica, pois, como já mencionado, tanto a superalimentação quanto a subnutrição podem ser prejudiciais ao crescimento e à formação do produto desejado. Do ponto de vista de engenharia de bioprocessos, tal questão configura-se como um problema de otimização no qual deseja-se maximizar a produtividade de bioproduto pela determinação do perfil ideal da variável de controle, no caso a taxa de alimentação de substrato (OLIVEIRA, 2018). Controlar uma fermentação batelada alimentada em sua trajetória de estados ótimos não é tarefa simples devido ao caráter dinâmico e não linear do biossistema e à sua representação matemática por meio de um conjunto de equações diferenciais não lineares de primeira ordem ou de ordem superior (RODRIGUEZ FERNANDEZ, 2006). Diversas técnicas de otimização têm sido propostas na literatura para resolver tais tipos de problemas. Entretanto, os métodos convencionais de otimização são frequentemente incapazes de funcionar bem para tais biossistemas, requerendo a utilização de métodos matemáticos e numéricos mais avançados para a resolução do problema (OLIVEIRA, 2018). 3.4 Controle Ótimo e Aplicação do Princípio do Máximo Um sistema sob controle é um sistema dinâmico, que evolui no tempo. A teoria de controle analisa as propriedades de tais sistemas, com o intuito de conduzi- los de um estado inicial até um estado final, respeitando eventualmente certas restrições. A origem de tais sistemas pode ser muito diversa: mecânica, elétrica, biológica, química, econômica, dentre outras. O objetivo pode ser o de estabilizar o sistema tornando-o insensível a certas perturbações (problema de estabilização) ou ainda determinar as soluções ótimas relativamente a um determinado critério de otimização (problema de controle ótimo). Para o desenvolvimento da teoria de controle pode-se recorrer a equações diferenciais, integrais, funcionais, de 15 diferenças finitas, determinísticas ou estocásticas, dentre outras. Portanto, a teoria de controle envolve vários domínios da matemática (CORON et al., 2004; BRYSON, 1996; SOTANG, 1990; MACKI et al., 1982), sendo uma área bastante complexa. Um sistema sob controle é dito controlável se ele puder ser conduzido, em tempo finito, de um determinado estado inicial até um estado final desejado. Em relação ao problema da controlabilidade, Kalman apresentou em 1949 um critério para caracterizar a controlabilidade de sistemas lineares de dimensão finita (Condição de Kalman). Para sistemas não lineares, o problema matemático da controlabilidade é muito mais difícil e constitui um domínio de investigação ainda ativo nos dias de hoje. Assegurada a propriedade de controlabilidade, pode-se desejar conduzir o sistema de um estado inicial a um estado final minimizando ou maximizando um determinado critério, o que configura então um problema de controle ótimo. A teoria de controle ótimo muitas vezes está incorporada a métodos matemáticos de otimização tais como aqueles baseados no Princípio do Máximo de Pontryagin (PMP). O método de otimização de bioprocessos empregando o PMP tem sido realizada desde os anos 60 (PALANKI et al., 1993). Embora essa técnica tenha sido aplicada na área de Biotecnologia, existem relatos de aplicações em outras áreas tais como química, farmacêutica e alimentícia. Por exemplo, no campo da química, Palanki et al. (1993) relatam várias aplicações do PMP na otimização de processos de polimerização de radicais livres em reatores batelada, determinando os perfis ótimos de temperatura e de adição de precursores. Por outro lado, Ciannella et al. (2012) realizou a simulação de um processo de craqueamento catalítico de diesel, aplicando cálculo variacional para a determinação do perfil de temperatura ótimo que proporcionasse máxima produção de gasolina no processo. Na área de Biotecnologia, o PMP tem sido amplamente utilizado para otimização de diversos bioprocessos. Palanki et al.(1993) relatam a aplicação do PMP para determinar a estratégia ótima de alimentação de glicose na produção de penicilina em reator fed-batch enquanto que em outra aplicação, a estratégia de alimentação ótima foi determinada para uma fermentação descontínua alimentada repetida visando maximizar a produção de massa celular. Guthke et al. (1981) usaram o PMP para determinar o perfil ótimo de concentração de substrato durante a produção do antibiótico turimicina em reator batelada com alimentação 16 intermitente, empregando para isso um modelo matemático simplificado para o crescimento microbiano e formação de produtos. Vanichsriratana (2000) aplicou o PMP para determinar o perfil ótimo de alimentação de substrato em fermentações batelada de metabólitos primários, exemplificando o estudo com o caso da maximização da produção de biomassa ao final do processo. Os autores incluíram o fator custo por unidade de tempo de operação justificando como necessária para a otimização de processos de produção de biomassa em reatores fed-batch. Skolpap et al. (2005) utilizaram o PMP para realizar estudos de otimização dos perfis de velocidade de alimentação de substrato em processos batelada alimentada visando a maximização da produtividade de penicilina no bioprocessos (BAYER et al., 2015; SBARCIOG et al., 2012; ZHANG et al., 2010; CRISTALDI et al., 2009). Lobato (2011) determinou o perfil ótimo de alimentação de substrato em um processo fed- batch de fermentação alcoólica, aplicando o PMP e procedimentos de redução do índice diferencial do sistema algébrico-diferencial para estabelecer diferentes estratégias de alimentação de substratos em função da concentração inicial de substrato, gerando novas possibilidades para a otimização do bioprocesso em escala industrial. A partir dos exemplos relatados, pode-se constatar que a aplicação do PMP é bastante atual, configurando-se como uma ferramenta matemática eficaz para a otimização de processos químicos e bioprocessos. Desse modo, torna-se necessário e justificável todo e qualquer aprofundamento nos fundamentos teóricos desse método tendo em vista o enorme potencial de aplicação que ele oferece. 17 4 METODOLOGIA 4.1 Estudo dos fundamentos teóricos do Princípio do Máximo de Pontryagin Nesta etapa, os fundamentos teóricos nos quais se baseia o PMP foram estudados visando à compreensão e ao domínio completo do método. Após tais estudos, elaborou-se uma formulação matemática simples visando facilitar a aplicação generalizada do PMP à otimização de bioprocessos baseada em modelos matemáticos e cálculo variacional. 4.2 Análise dos modelos matemáticos referentes aos bioprocessos escolhidos para estudo Nesta etapa realizou-se uma análise dos modelos matemáticos referentes aos bioprocessos escolhidos para estudo. A escolha dos bioprocessos foi realizada com base em referências que relatam aplicações do PMP à otimização desses processos, buscando abranger diferentes casos quanto aos objetivos de otimização e controle, variável manipulada, tipo de processo (batch, fed-batch ou contínuo) e bioproduto obtido (metabólito primário ou secundário). A análise dos modelos referentes a esses bioprocessos compreendeu o estudo detalhado da concepção, estrutura, variáveis, parâmetros, natureza e métodos de resolução das equações constitutivas de cada modelo. Os modelos estudados foram aqueles apresentados por Constantinides e Mostoufi (1999) para a produção em batelada de penicilina e por Costa et al., (1996) para a produção deste antibiótico em modo batelada alimentada. 4.3 Estudos de casos e simulação dos bioprocessos nas condições otimizadas Após a etapa de análise dos modelos matemáticos, estes foram usados para realizar estudos de casos a partir da simulação computacional dos bioprocessos nas condições otimizadas. 18 5 RESULTADOS E DISCUSSÃO 5.1 Estudo dos fundamentos teóricos do Princípio do Máximo de Pontryagin Alguns conceitos importantes do cálculo variacional bem como os fundamentos teóricos do PMP são apresentados no apêndice desta dissertação, o qual deve ser lido atenciosamente para uma melhor compreensão dos métodos matemáticos aplicados neste estudo. Também no apêndice, é apresentado um exemplo de cálculo da vazão de alimentação singular para um bioprocesso batelada alimentada de produção de metabólito primário, introduzindo-se na análise um fator custo por unidade de tempo operacional. 5.2 Análise dos modelos matemáticos referentes aos bioprocessos escolhidos para estudo 5.2.1 Modelo de Constantinides e Mostoufi (1999) O modelo matemático apresentado por Constantinides e Mostoufi (1999) para o crescimento de Penicillium chrysogenum, um microrganismo produtor de penicilina, é baseado na lei logística. Neste modelo, a limitação do crescimento microbiano não é descrita explicitamente em função da dinâmica de variação da concentração do substrato limitante, mas sim, implicitamente pela introdução do termo (1-X/Xm) na equação da taxa de crescimento (rx), a qual diminui ao longo do tempo e se torna nula quando X atinge o valor máximo Xm, condição em que o crescimento celular é encerrado. A produção de penicilina é também modelada, ao considerar que a cinética de formação de antibiótico é dada pelo modelo misto de Luedeking e Piret (rp=rx+βX), com =0 e β≠0 (β>0), uma vez que o antibiótico é um metabólito secundário produzido majoritariamente após a fase de crescimento Sinclair e Kristiansen (1987). Considera-se também o fato comumente observado em fermentações penicilínicas de que o produto é degradado por hidrólise a uma taxa proporcional à sua própria concentração, sendo este processo descrito por uma cinética de primeira ordem. Como a dinâmica de variação da concentração de O2 dissolvido não é descrita no modelo, conclui-se que o processo não é limitado por oxigênio. Da mesma forma, condições isotérmicas são implicitamente consideradas, uma vez que a temperatura também não é modelada. O modelo completo, constituído por duas equações diferenciais ordinárias correspondentes aos balanços de massa de células e de produto em um biorreator batelada, é dado a seguir: 19 𝑑𝑋 𝑑𝑡 = 𝑟𝑋 = 𝜇𝑚𝑋 (1 − 𝑋 𝑋𝑚 ) 𝑑𝑃 𝑑𝑡 = 𝑟𝑃 − 𝑟ℎ = 𝛽𝑋 − 𝑘ℎ𝑃 (1) onde 𝑡 é o tempo, 𝑋 é a concentração celular (% ms, ms = massa seca), 𝑃 é a concentração de antibiótico (U mL-1), 𝑟𝑋 é a taxa de crescimento celular (% ms h-1), 𝑟𝑃 é a taxa de produção de antibiótico (U mL-1 h-1), 𝑟ℎ é a taxa de hidrólise do produto (U mL-1 h-1) e 𝜇𝑚, 𝑋𝑚, 𝛽 e 𝑘ℎ são os parâmetros do modelo, os quais possuem os seguintes significados: 𝜇𝑚 é a velocidade específica máxima de crescimento (h-1), 𝑋𝑚 é a máxima concentração celular possível de ser atingida (% ms), 𝛽 é a constante de formação de produto não associada ao crescimento (U mL-1 h-1 % ms-1) e 𝑘ℎ é a constante cinética da reação de hidrólise do antibiótico (h-1). As variáveis do modelo foram então adimensionalizadas e expressões descrevendo os parâmetros cinéticos (bi) em função da temperatura () foram incorporadas ao modelo com o objetivo de ampliar sua faixa de validade para condições não isotérmicas. Na versão adimensionalizada do modelo, a hidrólise do antibiótico foi desconsiderada, sendo essa versão representada pelas seguintes equações: 13 2 2 1 2 1 11 1 yb d yd y b b yb d yd     (2) onde:             C20711940 C300050113 2501 01 2501 01 2501 01 o 654 o 321 2 62 2 62 532 32 2 32 422 32 2 32 11                            w ; .w ; .w w ; .w ; .w ww. ww. wb ; ww. ww. wb ; ww. ww. wb  (3) 5.2.2 Modelo de Costa et al. (1996) O modelo apresentado por Costa et al. (1996) é baseado no modelo clássico proposto por Bajpai e Reuss para o bioprocesso de produção de penicilina por Penicillium chrysogenum. Neste modelo, a velocidade específica de crescimento (µ) leva em conta limitações difusionais que ocorrem na biomassa fúngica filamentosa, 20 sendo descrita pelo modelo de Contois. A velocidade específica de formação de produto (π) considera a repressão da produção de penicilina pelo substrato, sendo representada pela equação de Andrews. A degradação da penicilina por hidrólise é também considerada, admitindo-se cinética de primeira ordem para esta reação, isto é, rh=kP. A velocidade específica de consumo de substrato (σ) é representada pelo modelo generalizado de Herbert-Pirt, segundo o qual o consumo de substrato é dado pelas parcelas consumidas para crescimento celular, manutenção e formação de produto. As equações do modelo considerando são as seguintes:     ; ; ; 2 / i m m S SPX/S f m k S Sk S PukPX dt dP m YY uSSX dt dS SBX S XuX dt dX            (4) Sendo: u=D=F/V (u é a variável de controle, D é a taxa de diluição, F é a vazão volumétrica de alimentação e V é o volume de cultura no biorreator). 5.3 Otimização de bioprocessos usando o Princípio do Máximo de Pontryagin Após a fase de análise dos modelos matemáticos referentes aos dois primeiros bioprocessos escolhidos para estudo, estes foram usados para ilustrar a aplicação do Princípio do Máximo de Pontryagin para a determinação do perfil ótimo da variável de controle (temperatura ou vazão específica de alimentação de substrato). Estes casos são apresentados a seguir. 5.3.1 Determinação do perfil ótimo de vazão específica de alimentação de substrato Para este estudo de caso, o processo de produção de penicilina foi representado pelo modelo matemático de Costa et al. (1996), dado em forma matricial por: 21    uXgXf dt Xd  (5) Sendo:       VFDu P SS X Xg kPX X X Xf P S X X f / ; ; ;                                          (6) As restrições para a variável de controle u são: maxmin uuu  (7) onde: VFu u /e0 maxmaxmin  As condições iniciais são dadas por:         0000 VVPP SSX X  0 ;0;0 ;0 . O objetivo do problema de otimização/controle é determinar o perfil temporal da variável de controle que maximiza a concentração de antibiótico ao final do processo. De acordo com o Princípio do Máximo de Pontryagin, o perfil ótimo deve maximizar o Hamiltoniano (H) , dado por:       321;  uXgXfH T T  (8) ou:                   Xgt XftH uttHH T T    0 0 . (9) Sendo              3 2 1 f f f Xf e              3 2 1 g g g Xg , chega-se a: 22    kPXXXffftH   3213322110     PSSXgggt f 321332211   (10) Para o controle ótimo tem-se que:  se   0t , então máxuu  ;  se   0t , então minuu  ;  se   ,0t então sinuu  . No intervalo singular:       0 ;0 ;0 0  tHttH  . Como  t  , tem-se que:              T TT dt d Xg X uXf X uxgxf Xdt d                            (11)         Xg dt d Xg dt d dt d Xgt T T T           (12)           uXgXfXg tdt Xd Xg X Xg dt d        (13)              uXgXg X XfXg X Xg dt d       (14) Substituindo (11) e (14) em (12):                                      uXgXg X XfXg X XgXg X uXf Xdt d TT   (15)                    XgXf X XfXg Xdt d T   (16) Em forma matricial: 23                                                                                                    3 2 1 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 3 2 1 3 3 2 3 1 3 3 2 2 2 1 2 3 1 2 1 1 1 321 - g g g X f X f X f X f X f X f X f X f X f f f f X g X g X g X g X g X g X g X g X g dt d   (17)   PXSXXXPgSSgXg f  321321 , , , , , (18) 0,0,1 3 1 2 1 1 1          X g X g X g (19) 0,1,0 3 2 2 2 1 2          X g X g X g (20) 1,0,0 3 3 2 3 1 3          X g X g X g (21)      S,X,S,S,X;kPXf,Xf,Xf   321 (22) 0,',' 3 1 2 1 1 1          X f X X f X X f SX  (23)   0,',' 3 2 2 2 1 2          X f X X f X X f SX  (24) k X f X X f X f S          3 3 2 3 1 3 ,',  (25) Então:                                                 3 2 1 3 2 1 321 ' 0 ' ' 0 ' ' - 1- 0 0 0 1- 0 0 0 1- g g g kX XX XXμμ f f f dt d S SX SX      (26)       ' '' ' - 321 21 2'1 3 2 1 321                            kgXgg XgXg XgXμg f f f dt d S SX SX      (27)        kgXggf XgXgfXgXμgf dt d S SXSX 32133 21222'111 ' '''      (28)         Sf SfXSfX XSS XSSXXSSX dt d ' '''.' 3 2 2 2 1      (29) No intervalo singular: 0 dt d (30)   03322110  ffftH  (31)   0332211  gggt  (32) 24 De (31):   1 3322 1332211 0 f ff fff     (33) Substituindo-se a Equação (33) na Equação (32), obtém-se:   3 1221 3113 2 1 1 2 2 331 1 3 2331 1 3 21 1 2 2 33221 1 33 1 1 22 33221 1 322 00                                                    gfgf gfgf g f f g gg f f gg f f g f f g ggg f f g f f ggg f ff (34) Introduzindo-se a expressão de 𝜆2 na de 𝜆1, chega-se a: 3 1 32 1221 3113 1                              f ff gfgf gfgf (35) Retornando à expressão de dt d :          0' '''' 3 2 3 2 3   Sf SfXSfX XSS XSSXXSSX dt d    (36)                XQ XSSXSSXXSSX dt d SfSfXSfX ' '''' 22 3    (37)     XQ XQ dt d 003    (Expressão do arco singular independente das variáveis adjuntas 1, 2, 3). Na expressão de  XQ , tem-se que:      2 ' 2 ' ; SBX SSBX SBX BS SBX S mm S m X m              (38) 25 2 2 2 2 2 1                            i m i mm i m ' S i m m k S Sk S k S k S Sk k S Sk S     (39) SP S SX S S SX X X SPX/S YY Y m YY / ' / ' ' / ' ' / ;        (40) Para a determinação da taxa de diluição singular (usin), inicialmente parte-se da definição de colchete de Lie para escrever a equação da derivada primeira de  , isto é:              Xg X Xf Xf X Xg XgfXgf dt d T       , ; ,  (41) A segunda derivada de  fica então dada por:       Xgf dt d Xgf dt d dt d T T ,, 2 2          (42) A derivada do colchete de Lie com relação ao tempo é dada por:                 uXgXfXgf Xd d dt Xd Xgf Xd d Xgf dt d                    ,,, (43) Substituindo T dt d        e    Xgf dt d , na expressão de 2 2 dt d  , obtém-se:                                                       ugf Xd gd ggf Xd d gf Xd fd fgf Xd d dt d TT ,,,, 2 2   (44) No arco singular, 2 2 dt d  =0 e u=usin. Assim:       0,,,, sin2 2  ugfggff dt d TT   (45) Utilizando um procedimento análogo ao utilizado na obtenção da expressão de  XQ , chega-se a: 26         Xb Xa uuXbXa dt d 3  sinsin2 2 0  (46) A expressão obtida para usin é bastante complexa e extensa, a qual foi implementada em linguagem de programação FORTRAN para a avaliação do valor desta variável ao longo do intervalo singular. A condição de parada da integração das equações de balanço de massa durante o período que sucede o intervalo singular, é determinada a partir do fato de que para problemas de tempo final livre, como é o caso,   0ftH . Assim:        dt dP dt dS dt dX HfffH XfHuXgXfH T u T 321332211 0      (47) Como as condições finais das variáveis adjuntas são 1e0,0 321   , a condição de parada em t=tf é dP/dt=0. Assim, a resolução do problema consistiu no seguinte algoritmo: 1) Integrar as equações de balanço de massa com 0u (operação em batelada) até os valores das variáveis de estado satisfazerem   0XQ , sendo o instante em que isto ocorre o 10 tempo de troca 1t ; 2) A partir do instante 1t , fazer sinuu  até o volume de meio no reator chegar ao máximo, o que corresponde ao 20 tempo de troca t2; 3) A partir do instante t2 , retomar a operação do reator com 0u até a condição de parada ser atingida em tf. Considerando os dados apresentados no trabalho de Costa et al. (1996) e sumarizados na Tabela 1, as equações de balanço de massa foram integradas numericamente para determinação de t1, isto é, o instante em que   0XQ . De acordo com os dados da Tabela 2, este instante de tempo pode ser considerado como sendo t1=28.74 h, uma vez que  XQ muda de sinal em torno deste valor. 27 Utilizando a equação deduzida para usin e por meio do programa desenvolvido em linguagem FORTRAN para o cálculo desta variável verificou-se que o valor de usin foi praticamente constante durante o intervalo singular sendo o valor de 0.00213h-1 bastante representativo. Com relação ao 20 tempo de troca, este foi determinado a partir de usin e dos volumes inicial e final (máximo) como segue: h71.9700213.0 121.8 10 lnln tttu V V udt V dV uVF dt dV i f                (48) Desta forma, a duração do intervalo singular é de 97.71 h e o 20 tempo de troca (t2) é t2= t1+∆t=28.74+97.71=126.45h. As equações de balanço de massa foram então integradas até este tempo e verificou-se que ao final da integração, a condição de parada tinha sido satisfeita, dispensando o período complementar de operação em batelada e atingindo-se uma concentração final de produto (P) de 6.35 g/L. Tabela 1. Parâmetros cinéticos e variáveis operacionais utilizadas no estudo de caso da produção de penicilina em modo batelada alimentada (Fonte: Costa, 1996) Parâmetros Cinéticos/Variáveis Operacionais Valores μm (h-1); B (g/L) 0.11; 6x10-3 πm (g.g-1.h-1); km (g/L); ki (g/L) 4x10-3; 1x10-4; 1x10-1 YX/S (g/g) ; YP/S (g/g) ; m(g.g-1.h-1) 0.47; 1.2; 2.9x10-2 X0 (g/L) ; S0 (g/L); P0 (g/L); 1.3; 69.0; 0.0 V0 (L) ; Vf (L) 8.121; 10.0 28 Tabela 2. Dados da integração numérica das equações de balanço de massa utilizados para determinação do 10 tempo de troca durante a produção de penicilina em modo batelada alimentada t (h) X (g/L) S (g/L) P (g/L)  XQ 28.742 30.09 4.16x10-3 1.39x10-2 125799.24 28.743 30.09 3.06 x10-3 1.40x10-2 -86336.84 5.3.1.1 Simulação do bioprocesso de produção de penicilina nas condições otimizadas – Caso 1 A simulação do bioprocesso nas condições otimizadas foi realizada em linguagem de programação FORTRAN. Para a integração numérica das equações diferenciais ordinárias correspondentes aos balanços de massa foi utilizado o método de Runge-Kutta-Gill de 4a ordem (CONSTANTINIDES; MOSTOUFI, 1999). Os perfis completos das variáveis de estado durante o bioprocesso estão apresentados nas Figuras 1 e 2. Figura 1. Perfis das variáves de estado S, X e P durante o bioprocesso de produção de penicilina em modo batelada alimentada 0 1 2 3 4 5 6 7 0 10 20 30 40 50 60 70 80 0 20 40 60 80 100 120 P ( g /L ) X , S ( g /L ) t (h) S X P 29 Figura 2.Perfil da variável de estado V e da variável de controle u durante o bioprocesso de produção de penicilina em modo batelada alimentada Perfis característicos de fermentações penicilínicas foram obtidos para as variáveis de estado, isto é, observa-se uma fase inicial de acumulação de células na qual o substrato é quase completamente consumido para esta finalidade. Após esta fase, o substrato alimentado é praticamente utilizado para a produção de penicilina uma vez que não há mais repressão catabólica da síntese do antibiótico pelos baixos níveis de concentração de substrato que se estabelecem no reator. Além disso, o padrão cinético observado está de acordo com aquele esperado para um metabólito secundário, isto é, a produção ocorre majoritariamente após o crescimento celular. 5.3.2 Determinação do perfil ótimo de temperatura durante a produção em batelada de penicilina Para este estudo de caso, o processo de produção de penicilina foi representado pelo modelo matemático proposto por Constantinides e Mostoufi (1999) em sua versão adimensionalizada e com parâmetros cinéticos (bi) dados em função da temperatura (), o qual foi apresentado no item 5.2. 30 Como neste caso   0Xg e u=0 devido ao reator ser operado em modo descontínuo (batelada), tem-se que:  Xf dt Xd  (49) Na Equação (49):                           13 2 1 2 1 11 2 1 ; yb y b b yb f f Xf y y X 2 1 (50) Como já visto no estudo de caso anterior, o Hamiltoniano é dado por:            132 2 1 2 1 11121 2 1 1 0 ; yby b b ybff f f H XfH u.XgXfH 21 1 2 2 TTT                           (51) Por outro lado, as taxas de variação temporal das variáveis adjuntas 1 e 2 são dadas por:                      0 2 321 2 1 11 2 1 by b b b d d X H d d      (52) A partir da equação anterior, as seguintes equações podem ser obtidas: d d b-y b b b d d 23 0 2 2 11 2 1 11 1        (53) A condição necessária para a otimização do bioprocesso é:   0 / 0 3 12 212 1 1 11                                    b y bb y b y HH (54) A partir das expressões dos bi em função da temperatura apresentadas anteriormente, chega-se a: 31                                        2 62 625321 2 32 3211 250,1 2 ;0 / ; 250,1 2 ww wwwb bb ww wwwb     (55) Substituindo os resultados anteriores na expressão de H/=0, obtém-se a expressão do perfil ótimo de temperatura  :         ww, wwy ww, wwy ww, wwwy ww, wwwy                    2 62 251 2 32 2111 2 62 62512 2 32 32111 2501 2 2501 2 2501 2 2501 2   (56) Como já demonstrado anteriormente, quando o objetivo é maximizar a concentração de antibiótico ao final do processo é necessário que   011   e   0.112  Sendo d2/dt=0, a segunda condição impõe que 2 seja constante e igual a 1.0 em todo o domínio do tempo, isto é, 2=1.0 para 0    1. Assim, o algoritmo para a resolução do problema consistiu dos seguintes passos: 1) Atribuição de um valor inicial para 1(0); 2) Integração do sistema de EDOs desde =0 até =1 e verificação se 1(1)=0. Caso não, atribuição de novo valor a 1(0) até a condição final ser satisfeita. 5.3.2.1 Simulação do bioprocesso de produção de penicilina nas condições otimizadas – Caso 2 O algoritmo proposto foi implementado em linguagem de programação FORTRAN e os perfis da variável adjunta 1 e das variáveis de estado (y1, y2 e  ) estão mostrados nas Figuras 3 a 6. Os valores iniciais das variáveis de estado foram aqueles utilizados por Constantinides e Mostoufi (1999): y1(0)=0.03 e y2(0)=0. Encontrou-se grande dificuldade em determinar, por tentativa e erro, o valor de 1 em =0, o qual foi determinado como sendo igual a 3.61210035. Em razão disso, o método de Newton-Raphson para solução de equações algébricas não lineares (CONSTANTINIDES; MOSTOUFI, 1999) foi acoplado ao método de integração, 32 tornando o algoritmo computacional autônomo para a determinação de 1(0), o qual mostrou-se eficaz para resolver o problema formulado uma vez que a condição de contorno 1=0 foi satisfeita ao final do bioprocesso (=1), conforme pode-se observar na Figura 3. Figura 3. Perfil da variável adjunta 1 durante uma fermentação penicilínica não isotérmica O perfil de concentração celular mostrado na Figura 4 incorpora as principais fases envolvidas em uma curva de crescimento microbiano típica, isto é, as fases exponencial, estacionária e de declínio. A fase de declínio pode ser atribuída aos efeitos negativos das baixas temperaturas sobre os parâmetros cinéticos de crescimento. Na Figura 5, referente à dinâmica de produção de penicilina, pode-se observar uma fase lag inicial curta, seguida por uma fase de transição na qual a produção de penicilina é iniciada até que uma fase final de produção aproximadamente linear é alcançada. De acordo com a formulação apresentada (modelo do processo e Princípio do Máximo de Pontryagin), o perfil ótimo de temperatura é declinante, variando entre 30 e 20oC segundo o comportamento mostrado na Figura 6 (a). Esse perfil sugere uma temperatura operacional variável durante as fases de crescimento e produção de penicilina, contradizendo a prática industrial de manter constante a temperatura 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0   ( -)  (-) 33 durante todo o bioprocesso. Particularmente durante a fase de produção da penicilina, o perfil prescreve uma diminuição na temperatura de operação para que uma alta concentração de antibiótico seja atingida ao final do bioprocesso. Figura 4. Perfil de concentração adimensional de células durante uma fermentação penicilínica não isotérmica Figura 5. Perfil de concentração adimensional de produto durante uma fermentação penicilínica não isotérmica 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 y 1 ( co n ce n tr aç ão c e lu la r)  (-) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y 2 ( co n ce n tr aç ão p ro d u to ) (-) 34 0.0 0.2 0.4 0.6 0.8 1.020 22 24 26 28 30  [ te m p e ra tu ra ( o C )]  (-) (a): perfil exato (b): perfil aproximado Figura 6. Perfil de temperatura durante uma fermentação penicilínica não isotérmica O perfil de temperatura (a) mostrado na Figura 6 pode trazer alguma dificuldade prática à sua programação/execução. Nesse contexto, um perfil aproximado, derivado do perfil exato, como o representado pela curva (b) na Figura 6, pode tornar a estratégia de programação de temperatura factível. Esta proposta está de acordo com aquela sugerida por Bailey e Ollis (1986), ou seja, a programação de temperatura delineada pelos cálculos pode ser aproximada na prática industrial com pouco custo agregado. No caso, o perfil aproximado foi construído a partir das seguintes equações baseadas na análise dos dados de temperatura gerados pelo perfil exato:         62511652900 0040 4000 0040 00 ... .. C)(:.. ).().( ).( o               0259040 .C)(:.. o       905002590 9001 0190 9001 90 ... .. C)(:.. ).().( ).( o                   35 5.3.2.2 Simulação do bioprocesso de produção de penicilina em outras condições diferentes da ótima – Caso 2 Uma simulação natural a ser realizada é aquela em condição isotérmica visando verificar se a operação com temperatura constante é de fato menos produtiva em antibiótico do que a operação à temperatura variável. Desta forma, foram realizadas simulações para temperaturas constantes de 20, 25 e 30 oC e os resultados foram confrontados com aqueles obtidos para temperatura variável segundo o perfil ótimo estabelecido pelo método de otimização. A Figura 7 mostra os resultados destas simulações, as quais evidenciam o fato bem conhecido de que altas temperaturas (30 oC) favorecem o crescimento do fungo enquanto que baixas temperaturas (20 oC) favorecem a síntese do antibiótico (BAILEY; OLLIS, 1986). Observa-se também que a operação isotérmica em uma temperatura intermediária às investigadas apresenta resultados bastante semelhantes aos da operação com perfil de temperatura otimizado, configurando-se, assim, como uma alternativa bastante viável para o caso em que a estratégia de temperatura variável não possa ser praticada, embora com produtividade menor como mostram os dados apresentados na Tabela 3. O bom desempenho do bioprocesso a 25 oC deve-se ao fato de que grande parte do perfil ótimo se desenvolve em torno desta temperatura (Figura 6), explicando assim a semelhança entre o comportamento observado para a operação com temperatura variável e aquele com temperatura constante de 25 oC. Apesar de a operação a uma temperatura fixa entre 24 e 25 oC predominar na prática industrial, os benefícios da programação de temperatura durante a fermentação penicilínica são evidentes. Tabela 3. Valores finais da concentração de células (y1) e de produto (y2) em fermentações penicilínicas isotérmicas e não isotérmicas  (oC) y1 (=1) y2 (=1) y2/y1 (=1) ótimo 0.82 1.22 1.49 20 0.53 0.65 1.23 25 0.94 1.18 1.25 30 1.07 0.80 0.75 36 Figura 7. Perfis de concentração de células e de produto durante fermentações penicilínicas isotérmicas a diferentes temperaturas 5.3.2.3 Modelo matemático com restrição – correção da curva de crescimento Revisitando a Figura 3, observa-se um comportamento anômalo da curva de crescimento microbiano que é o declínio da concentração de células ao final do bioprocesso, assemelhando-se a um processo de morte celular, não considerado no modelo matemático. A explicação para esse comportamento é que o valor de b2 representa a máxima concentração celular que pode ser atingida em um tempo infinito, à temperatura constante. Quando se controla o valor de b2 com a equação 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2  ótimo 25 o C 20 o C y 1 ( co n ce n tr aç ão c e lu la r)  (-) 30 o C 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 25 o C  ótimo 30 o C 20 o C y 2 ( co n ce n tr aç ão p ro d u to ) (-) 37 b2=b2(), sob temperatura variável, poderá existir um ponto no domínio temporal no qual b2 pode se tornar menor que a concentração celular (y1) nesse instante, resultando em dy1/dt < 0 e consequente declínio da concentração celular (y1). Este comportamento matemático não é incoerente com eventos reais que podem ocorrer durante o processo tais como morte ou autólise celular, os quais provocam decréscimo na concentração de células. Analisando-se o arquivo de dados que gerou a Figura 3, verifica-se que tal declínio inicia-se no instante =0.7. Se o declínio na concentração celular não for permitido durante a solução do problema, então a restrição y1-b2  0 deve ser incluída na formulação do modelo matemático visando torná-lo coerente com a curva logística postulada para o crescimento fúngico. Quando a restrição é alcançada em um dado ponto da curva de crescimento, a concentração celular atinge seu valor máximo, devendo a partir deste ponto manter-se constante neste valor e não assumir valores inferiores a este sob pena de se consumar um declínio na concentração celular em instantes posteriores do processo, tal como observado na Figura 3. Durante esta fase estacionária, a equação de estado referente à concentração celular (y1) torna-se: 𝑑𝑦1 𝑑𝑡 = 0 (57) Em decorrência da Equação (57), tem-se que f1=0. Assim, o Hamiltoniano do sistema (H) e a equação de estado para a variável adjunta 1 reduzem-se às Equações (58) e (59), respectivamente. 𝐻 = 2𝑏3𝑦1 (58) 𝑑1 𝑑𝑡 = − 𝜕𝐻 𝜕1 = − 2𝑏3 (59) A condição necessária para a otimização do bioprocesso fica então dada por: 𝜕𝐻 𝜕𝜃 = 2 𝑦1 𝜕𝑏3 𝜕𝜃 = 0 (60) 38 Substituindo-se na Equação (60) a expressão de 𝜕𝑏3/𝜕𝜃 dada pela Equação (55), obtém-se: 𝜕𝐻 𝜕𝜃 = 2 𝑦1 [− 2𝑤5𝑤2(𝜃 − 𝑤6) (1 − 𝑤2(25 − 𝑤6)2) ] = 0 (61) Decorre da Equação (61) que o perfil ótimo de temperatura durante a fase estacionária de crescimento será dado por: 𝜃 = 𝑤6 = 20 𝑜𝐶 (62) Do ponto de vista computacional, haverá então dois conjuntos de equações: um válido para a parte não restrita do modelo e outro válido para a parte restrita, quando a restrição é atingida. Nas Figuras 8 a 11 estão apresentados os perfis temporais referentes à resolução do problema de valor no contorno incorporando a restrição sobre a concentração celular. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y 1 ( co n ce n tr aç ão c e lu la r)  (-) Figura 8. Perfil de concentração adimensional de células durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular 39 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y 2 ( co n ce n tr aç ão p ro d u to )  (-) Figura 9. Perfil de concentração adimensional de produto durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular 0.0 0.2 0.4 0.6 0.8 1.0 18 20 22 24 26 28 30   t e m p e ra tu ra ( o C )] (-) Figura 10. Perfil de temperatura durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular 40 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0   ( -)  (-) Figura 11. Perfil da variável adjunta 1 durante uma fermentação penicilínica não isotérmica – modelo incorporando restrição sobre a concentração celular Analisando-se a Figura 8, observa-se que a incorporação da restrição ao conjunto de equações resolveu o problema anteriormente detectado de declínio da concentração celular (Figura 3), gerando uma nova curva de crescimento que é coerente com aquela preconizada pelo modelo matemático. Com relação ao perfil de concentração de produto (Figura 9), este apresentou comportamento similar ao obtido anteriormente (Figura 4) com a diferença de que a partir do instante no qual a restrição é alcançada, em torno de  = 0.7, a concentração de produto aumenta linearmente com o tempo até o final do processo. Quanto ao perfil ótimo de temperatura (Figura 10), observa-se a ocorrência de dois períodos distintos: um primeiro período no qual a temperatura se mantém relativamente alta para favorecer o crescimento celular em detrimento da produção de antibiótico, estendendo-se este período até a restrição de crescimento nulo ser alcançada em  = 0.7, quando então a temperatura ótima sofre um degrau negativo, decrescendo abruptamente a 20 oC e dando início ao segundo período no qual a temperatura é mantida constante neste valor relativamente baixo de 20 oC até o final do processo para favorecer a biossíntese do antibiótico. Analisando-se a Figura 11, verifica-se que ao final do bioprocesso (=1.0), a condição de contorno imposta à variável adjunta 1 é cumprida, isto é, 1(1)=0. 41 6 CONCLUSÕES Neste estudo, demonstrou-se a utilidade do Princípio do Máximo de Pontryagin para a otimização de bioprocessos complexos de produção de antibióticos tais como aqueles conduzidos em modo batelada alimentada ou modo batelada sob condições não isotérmicas. A partir da aplicação deste princípio foi possível determinar, utilizando-se métodos numéricos simples tais como os de Runge-Kutta-Gill e de Newton- Raphson, o perfil ótimo de temperatura em cultivos em batelada e de vazão específica de alimentação de substrato em cultivos em batelada alimentada que maximizam a concentração de antibiótico ao final do processo. O Princípio do Máximo de Pontryagin mostrou ser uma poderosa e adequada ferramenta para a otimização e o controle de bioprocessos visando à máxima produtividade de bioproduto. Entretanto, para a aplicação deste princípio é necessário dispor de um modelo matemático, preferencialmente fenomenológico e representativo do bioprocesso, que permita avaliar a factibilidade da solução encontrada para o problema. No presente estudo, dois modelos fenomenológicos clássicos de bioprocessos de produção de penicilina foram utilizados, juntamente com o Princípio do Máximo de Pontryagin, para determinar as condições operacionalmente ótimas para a produção desse antibiótico e as soluções encontradas são consideradas factíveis, podendo ser implementadas na prática industrial sem grandes dificuldades. Também, constatou-se que para a resolução do problema de controle ótimo em bioprocessos conduzidos em modo batelada alimentada é necessário definir um modelo matemático adequado que permita a determinação da taxa de alimentação singular. 42 6 REFERÊNCIAS BAILLEY, J.E.; OLLIS D.F. Biochemical Engineering Fundamentals. 2nd. ed. New York: McGraw- Hill, 1986. BAYEN, T.; MAZADE, M.; MAIRET, F. Analysis of an optimal control problem connected to bioprocesses involving a saturated singular arc. Discrete and Continuous Dynamical Systems - Series B, American Institute of Mathematical Sciences, v. 20, n. 1, p.39-58, 2015. BERBER, R.; PERTEV, C.; TUERKER, M. Optimization of feeding profile for baker's yeast production by dynamic programming. Bioprocess Engineering, v. 20, p. 263- 269, 1998. BERKHOLZ, R.; GUTHKE, R. Model Based Sequential Experimental Design for Bioprocess Optimization. In: HOFMAN, M.; THONART, P. (Editors). Engineering and Manufacturing for Biotechnology. Heidelberg: Springer Netherlands, 2002. p. 129-141. Disponível em: https://link.springer.com/book/10.1007/0-306-46889-1. BOSCAIN, U.; PICCOLI, B. A short introduction to optimal control. In: SARI, Tewfik (Ed.). Controle Non Linéaire et Applications. Travaux en cours. Hermann, Paris: Les Cours du Cimpa, 2005. P. 19-66. BRYSON, A.E. Optimal control-1950 to 1985. Control Systems Magazine, v. 16, p. 26-33, 1996. CHEN, L. Z.; NGUANG, S. K.; CHEN, X. D. Modelling and Optimization of Biotechnological Processes: Artificial Intelligence Approaches. Berlin: Springer- Verlag, 2006. p.1-6. CHEN, W.C. Medium improvement for β-fructofuranosidase production by Aspergillus japonicus. Process Biochemistry, v. 33, n. 3, p. 267-271, 1998. CIANNELLA, S; ALVES J. Modelagem e Simulação do Processo de Craqueamento Catalítico do Gasóleo Aplicando Cálculo Variacional para Determinação do Perfil 43 Ótimo de Temperatura. In: ENCONTRO NACIONAL DE EDUCAÇÃO, CIÊNCIA E TECNOLOGIA UEPB, v.1, n.1, 2012, Campina Grande, Pb. Anais [...]. Campina Grande: Universidade Estadual da Paraíba, 2012. p. 1-11. CONSTANTINIDES, A.; MOSTOUFI, N. Numerical Methods for Chemical Engineers with MATLAB Applications. Upper Saddle River, NJ: Prentice Hall PTR, 1999. CORON, J. M.; TRELAT, E. Tout est sous controle. Spécial Recherche, p.126–131, 2004. COSTA, A. C.; LIMA, E. L.; ALVES, T. L. M. Otimização de processos fermentativos em batelada alimentada através da teoria de controle singular. In: SIMPÓSIO NACIONAL DE FERMENTAÇÕES, XI, São Carlos, SP, 1996. Anais [...]. São Carlos: UFSCar, v.1, p. 357-362, 1996. COSTA, A.C. Singular control in bioreactors. Rio de Janeiro: Universidade Federal do Rio de Janeiro (UFRJ), 1996. CRISTALDI, M.; GRAU, R.; MARTINEZ, E. Iterative Design of Dynamic Experiments in Modeling for Optimization of Innovative Bioprocesses. Chemical Product and Process Modeling: Modeling and Control, v. 4, n. 2, Article 6, 28 pages, 2009. DI SERIO, M; TESSER, R.; SANTACESARIA, E. A kinetic and mass transfer model to simulate the growth of baker’s yeast in industrial bioreactors. Chemical Engineering Journal, v.82, p.347–354, 2001. FERREIRA, R. A. C.; TORRES, D. F. M. Higher-order calculus of variations on time scales. In: SARYCHEV, A.; SHIRYAEV, A.; GUERRA, M.; GROSSINHO, M.R. (Eds). Mathematical Control Theory and Finance. Berlin: Springer, 2008. 149–159. GOLDSTINE, H.H. A history of the calculus of variations from the 17th through 19th century. New York,N.Y.: Springer-Verlag, 1980. (Studies in the History of Mathematics and Physical Sciences, v. 5). 44 GUTHKE, R; KNORRE W. Optimal Substrate Profile for Antibiotic Fermentations, Biotechnology and Bioengineering, v. 23, n. 12, p. 2771-2777, 1981. KENNEDY, M.; KROUSE, D. Strategies for improving fermentation medium performance. Journal of Industrial Microbiology & Biotechnology, v. 23, p. 456– 475, 1999. LEE, J.; LEE, S.Y.; PARK, S; MIDDELBERG, A.P.J. Control of fed-batch fermentations. Biotechnology Advances. v. 17, p.29–48, 1999. LEI, F.; ROTBOLL, M.; JORGENSEN, S. B. A biochemically structured model for Saccharomyces cerevisiae. Journal of Biotechnology, v. 88, p. 205-221, 2001. LIM, H.C.; LEE, K.S. Process Control and Optimization. In: PONS, M.N. (Ed). Bioprocess Monitoring and Control. München: Hanser Publishers, 1992. p. 159-222. LOBATO, F.S. Determinação do Perfil Ótimo de Alimentação de Substrato no Processo de Fermentação Alcoólica-Influência da Condição Inicial, Matemática Aplicada e Computacional. TEMA Tend. Mat. Apl. Comput. v. 12, n. 1, p. 1-10, 2011. MACKI, J. W.; STRAUSS, A. Introduction to optimal control theory. New York, N. Y.: Springer, 1982. MARINS, N.; TORRES, D.F.M. Calculus of Variations on Time Scales with Nabla Derivatives. Nonlinear Analysis, v. 71, p. e763–e773, 2009. MONROY PEREZ, F. Historical perspective of the Pontryagin Maximum Principle. Lecturas Matemáticas, v. 37, n. 2, p. 117-169, 2016. MONTAÑO, C.D.I. Otimização dinâmica do cultivo semi-contínuo de Pichia pastoris recombinante para produção das enzimas heterólogas alfa amilase e 45 penicilina G Acilase. 2010. 111p. Dissertação (Mestrado em Engenharia Química) Engenharia Química, Universidade Federal de São Carlos, São Carlos, 2010. OLIVEIRA, S.C. Model-based evolutionary operation design for batch and fed-batch antibiotic production bioprocess. In: SILVA, V. (Ed.). Design of experiments applied to industrial chemical process. London: INTECH. 2018. Charpt 4, p. 43-63. PALANKI, S; KRAVARIS C; WANG H. Y. Synthesis of State Feedback Laws for End- Point Optimization in Batch Processer. Chemical Engineering Science, v. 48, n. 1, p. 135-152, 1993. PERTEV, C.; TURKER, M.; BERBER, R. Dynamic modeling, sensitivity analysis and parameter estimation of industrial yeast fermenters. Computers & Chemical Engineering, v 21, p S739-S744, 1997. RAMIREZ, W. F. Process Control and Identification. Boston: Academic Press, 1994. RODRIGUEZ ALFARO, L.H.; ALCORTA GARCIA, E. From Euler-Lagrange system representation to the generalized Hamiltonian one. Revista Electrónica Nova Scientia, v. 7, n. 14, p. 1-23, 2015. RODRIGUEZ FERNANDEZ M. Modelado e Identificaion de Bioprocesos. Tese (Doutorado em Ingeniería Química) Memoria Universidad de Vigo, 2006. RUZ, M. R. El Principio del Máximo de Pontryagin. 2014. 51 p. Treball final de Grau. (Grau de Matemàtiquqes) Facultat de Matemàtiques, Universitat de Barcelona, 2014. SBARCIOG, M.; LOCCUFIER, M.; WOUWER, A. V. An Optimizing Start-up Strategy for a Bio-methanator. Bioprocess and Biosystems Engineering, v. 35, n. 4, p. 565- 578, 2012. 46 SIRCAR, A.; SRIDHAR, P.; DAS, P.K. Optimization of solid state medium for the production of clavulanic acid by Streptomyces clavuligerus. Process Biochemistry, v. 33, n. 3, p. 283-289, 1998. SKOLPAP, W.; SCHARER, J. M.; DOUGL, P. L.; MOO-YOUNG, M. Optimal Feed Rate Profiles for Fed-Batch Culture in Penicillin Production. Songklanakarin J. Sci. Technol., v. 27, n. 5, p. 1057-1064, 2005. SMETS, I.Y.M., VERSYCK, J.E., VAN IMPE, J. F. M. Optimal control theory: A generic tool for identification and control of (bio-)chemical reactors. Annual Reviews in Control. v. 26, n. 1, p 57-73, 2002. SONNLEITNER, B.; KAPPELI, O. Growth of Saccharomyces cerevisiae is controlled by its limited respiratory capacity: Formulation and verification of a hypothesis. Biotechnology and Bioengineering, v. 28, p 927-937, 1986. SONTAG, E. D. Mathematical control theory. New York, N. Y.: Springer, 1990. VANICHSRIRATANA, W. Optimization of a Primary Metabolite Fermentation Process: Effect of Cost Factor on the Optimal Feed Rate Control. Kasetsart Journal (Natural Science), v. 34, p. 495-499, 2000. ZHANG, H.; ZHANG, Z.; LAN, L. H. Evolutionary optimization of a fedbatch penicillin fermentation process. In: INTERNATIONAL SYMPOSIUM ON COMPUTER, COMMUNICATION, CONTROL AND AUTOMATION, Tainan, Taiwan, 5-7 May Proceedings […]. 2010. v. 1, p. 403-406. ISBN: 978-1-4244-5568-3 47 APÊNDICE A.1 Estudo dos fundamentos teóricos do Princípio do Máximo de Pontryagin Para compreender o PMP, primeiramente devem ser apresentados alguns conceitos de cálculo variacional e depois prosseguir com a teoria de controle. O cálculo variacional trata do problema de determinação de extremos, isto é, máximos ou mínimos, de funcionais. Um funcional é uma função real cujo domínio é um espaço de funções. Mais explicitamente, um funcional associa um número real a cada função de um certa classe de funções para as quais o funcional está definido (MONROY PEREZ, 2016; MARINS; TORRES, 2009; FERREIRA; TORRES, 2008; GOLDSTINE, 1980). Por simplicidade e compreensão será considerado inicialmente o problema unidimensional. Sejam (𝑥1, 𝑦1) e (𝑥2, 𝑦2) dois pontos distintos num plano, com 𝑥2 > 𝑥1. Seja 𝑦(𝑥) uma curva ligando esses dois pontos, o comprimento de arco da curva entre os pontos inicial e final é: 𝑙[𝑦] = ∫ √1 + ( 𝑑𝑦 𝑑𝑥 ) 2 𝑑𝑥 𝑥2 𝑥1 (A1) O comprimento de arco 𝑙 é um funcional de 𝑦, isto é, dada uma função continuamente diferenciável 𝑦(𝑥) a ela se faz corresponder um único numero real 𝑙[𝑦] definido pela Equação (A1). O mais simples problema de cálculo variacional pode ser assim formulado: a) Dado o funcional: 𝐽[𝑦] = ∫ 𝑓(𝑦(𝑥), 𝑦′(𝑥), 𝑥) 𝑥2 𝑥1 𝑑𝑥 (A2) onde 𝐽 é uma função real cujo domínio é o espaço das funções, sendo 𝑓: ℝ3 → ℝ uma função conhecida (MONROY PEREZ, 2016; RUZ, 2014). b) Determinar, dentre todas as curvas continuamente diferenciáveis 𝑦(𝑥) que passam pelos pontos fixos (𝑥1, 𝑦1) e (𝑥2, 𝑦2), aquela que minimiza (extremiza) 𝐽. 48 A ideia é testar várias funções no integrando até que uma dada função forneça um valor que extremiza o funcional 𝐽. A fim de encontrar a curva que extremiza 𝐽, considera-se uma curva vizinha �̅� definida a partir de pequenas variações ou mudanças na função 𝑦(𝑥) visando analisar como o valor do funcional 𝐽 é afetado. �̅�(𝑥) = 𝑦(𝑥) + 𝛼𝛿(𝑥) (A3) Desde que haja interesse em realizar pequenas variações em 𝑦(𝑥), introduz- se o parâmetro 𝛼 para controlar tais variações, isto é, o parâmetro 𝛼 possibilita modificar a função 𝑦(𝑥) pela introdução de uma função arbitrária 𝛿(𝑥). O parâmetro arbitrário 𝛼 é um número real e 𝛿 é uma função continuamente diferenciável que se anula em 𝑥 = 𝑥1 e 𝑥 = 𝑥2, isto é: 𝛿(𝑥1) = 𝛿(𝑥2) = 0 (A4) As condições dadas pela Equação (A4) são necessárias para que a curva variada �̅�(𝑥) também passe pelos pontos extremos 𝑥 = 𝑥1 e 𝑥 = 𝑥2. Substituindo-se 𝑦 por �̅� na Equação (A2), obtém-se como resultado uma função de 𝛼, dada por: 𝛷(𝛼) = 𝐽[�̅�] = ∫ 𝑓(�̅�(𝑥), �̅�′(𝑥), 𝑥) 𝑥2 𝑥1 𝑑𝑥 (A5) Uma condição necessária para que �̅�(𝑥) extremize 𝐽 quando 𝛼 = 0, considerando que 𝜕�̅� 𝜕𝛼 = 𝛿 e �̅�(𝛼 = 0) = 𝑦, é: ( 𝜕𝛷 𝜕𝛼 ) 𝛼=0 = ∫ ( 𝜕𝑓 𝜕𝑦 𝛿 + 𝜕𝑓 𝜕𝑦´ 𝛿´) 𝑥2 𝑥1 𝑑𝑥 = 0 (A6) A fim de se obter uma equação para 𝑦(𝑥), é necessário explorar a arbitrariedade de 𝛿. Para isto é preciso eliminar 𝛿´ da Equação (A6), o que pode ser conseguido fazendo-se uma integração por partes e impondo-se as condições 𝛿(𝑥1) = 𝛿(𝑥2) = 0. A Equação (A6) reduz-se então a: ∫ [ 𝜕𝑓 𝜕𝑦 − 𝑑 𝑑𝑥 ( 𝜕𝑓 𝜕𝑦´ )] 𝛿𝑑𝑥 = 0 𝑥2 𝑥1 (A7) 49 De acordo com o Lema Fundamental do Cálculo Variacional, para que seja cumprida a Equação (A7) é necessário que: [ 𝜕𝑓 𝜕𝑦 − 𝑑 𝑑𝑥 ( 𝜕𝑓 𝜕𝑦´ )] = 0 (A8) A Equação (A8) é denominada de Equação de Euler (MONROY PEREZ, 2016; RUZ, 2014), e sua solução fornece a função 𝑦 e consequentemente a função 𝑓 que atribui ao funcional 𝐽 um valor estacionário. A equação diferencial de Euler é uma condição apenas necessária para a existência de um extremo (máximo ou mínimo). Trata-se de uma equação diferencial de segunda ordem cuja solução geral contém duas constantes arbitrárias, permitindo satisfazer as condições de contorno de forma que a curva passe pelas extremidades fixas (BOSCAIN; PICCOLLI, 2005). Pode-se estender o funcional 𝐽 para 𝑛 dimensões e também ampliar a noção de variação e encontrar as equações de Lagrange a partir do princípio variacional, como descrito a seguir. Primeiramente, as seguintes mudanças de notação são feitas: 𝑥 → 𝑡, 𝑦 → 𝑞, 𝑦´ = 𝑑𝑦 𝑑𝑥 → �̇� = 𝑑𝑞 𝑑𝑡 , 𝑓 → 𝐿, 𝐽 → 𝑆 Secundariamente, calcula-se a variação de 𝑆 da mesma forma que foi feito para maximizar 𝐽: 𝛿𝑆 ≡ 𝛿 ∫ 𝐿(𝑞, �̇�, 𝑡)𝑑𝑡 = 0 𝑡2 𝑡1 (A9) Na Equação (A9): 𝑆 = ∫ 𝐿(𝑞1, … , 𝑞𝑛, �̇�1, … , �̇�𝑛, 𝑡)𝑑𝑡 = 0 𝑡2 𝑡1 (A10) Finalmente, para maximizar 𝛿𝑆 = 0, as equações de Lagrange devem ser cumpridas: 𝑑 𝑑𝑡 ( 𝜕𝐿 𝜕𝑞�̇� ) − 𝜕𝐿 𝜕𝑞𝑘 = 0, 𝑘 = 1, … , 𝑛 (A11) 50 A.1.1 Descrição do problema de otimização Um problema de controle ótimo pode ser formulado do seguinte modo: Seja um sistema sob controle, cujo estado num determinado instante é representado por um vetor. Os controles são funções ou parâmetros, geralmente sujeitos a restrições, que agem sobre o sistema alterando sua dinâmica. Um conjunto de equações diferenciais, descrevendo as taxas de variação temporal das variáveis de estado, é usado para modelar o comportamento dinâmico do sistema. Posteriormente, é necessário utilizar a informação presente e as características do problema para construir os controles adequados que permitirão atingir um objetivo específico. São impostas restrições sobre a trajetória ou sobre os controles, sendo a consideração de tais restrições imprescindível para a resolução do problema. Fixa- se um critério, geralmente na forma de um funcional que depende do estado do sistema e dos controles, o qual deve ser minimizado/maximizado simultaneamente à satisfação das condições anteriores. Note-se que a forma da trajetória ótima depende fortemente do critério de otimização, isto é, a trajetória a ser seguida difere se o objetivo for realizar a operação em tempo mínimo ou minimizar/maximizar um outro critério operacional (MONTAÑO, 2010; BOSCAIN; PICCOLLI, 2005; COSTA 1996). Em geral, os bioprocessos a serem otimizados são dinâmicos. Do ponto de vista matemático, considera-se um sistema como dinâmico quando seu estado varia com o tempo (MONTAÑO, 2010; COSTA, 1996). Na ausência de variações espaciais, tais sistemas são regidos por equações diferenciais ordinárias, as quais descrevem o comportamento dinâmico do sistema a ser controlado. Essas equações podem ser dadas na seguinte forma: �̇� = 𝑑𝑥 𝑑𝑡 = 𝑓(𝑥(𝑡), 𝑢(𝑡), 𝑡) (A12) Na Equação (A12), 𝑥(𝑡) é o vetor das variáveis de estado do sistema de dimensão 𝑛, 𝑡 é o tempo, 𝑢(𝑡) é o vetor das variáveis de controle de dimensão 𝑚, que especifica as variáveis de entrada no modelo do sistema, 𝑓 é o vetor de funções (𝑛 × 1) e �̇� é o vetor das derivadas temporais de estado, calculadas pelas equações diferenciais que constituem o modelo. As letras ou símbolos com um traço em baixo e minúsculas indicam vetores e maiúsculas, matrizes (MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). 51 Além das equações de estado, tem-se um conjunto de condições de contorno e restrições sobre as variáveis de estado, no tempo inicial, final, e ao longo da trajetória, conforme Equação (A13). Nesta formulação, o modelo do sistema é reescrito como restrições de igualdade do problema de otimização Equação (A14). 𝑥(𝑡0) = 𝑥 0 , 𝑆 (𝑥(𝑡), 𝑢(𝑡)) ≤ 0, 𝑇(𝑥(𝑡𝑓)) ≤ 0 (A13) 𝑔 = 𝑓(𝑥, 𝑢, 𝑡) − �̇�(𝑡) (A14) Nas Equações (A13) e (A14), 𝑥0 é o vetor de condições inicias, que geralmente é conhecido, 𝑔 é o vetor de restrições de igualdade sobre a trajetória (o modelo), 𝑆 é o vetor de restrições de desigualdade sobre a trajetória e 𝑇 é o vetor de restrições de desigualdade no estado final (MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). No problema de controle ótimo de sistemas dinâmicos, busca-se o extremo de um índice de desempenho 𝐽, o qual é função das variáveis de estado 𝑥 e das variáveis de controle 𝑢, as quais são função do tempo 𝑡 (MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). O índice de desempenho deve ser escolhido de forma que a otimização resulte em uma lei de controle implementável. A forma geral do índice de desempenho é: 𝐽 = ℎ(𝑥(𝑡𝑓), 𝑡𝑓) + ∫ 𝐹(𝑥, 𝑢, 𝑡)𝑑𝑡 𝑡𝑓 𝑡0 (A15) Na Equação (A15) 𝐹: ℝ𝑛 → ℝ é um funcional linear enquanto que o termo ℎ(𝑥(𝑡𝑓), 𝑡𝑓) é a contribuição do estado final para o índice de desempenho, elemento presente com frequência em problemas de controle ótimo (BAYEN et al., 2015; BARCIOG et al. 2012, MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). Na Equação (A15), o primeiro termo é chamado de custo terminal e o segundo termo é conhecido como integral do custo (SMETS et al. 2002). 52 A resolução do problema de controle ótimo visa determinar a política de controle que maximizará ou minimizará um índice de desempenho, sujeito a restrições impostas pela natureza fenomenológica do problema. Uma vez que este é um problema de otimização com restrições, pode-se utilizar o Lagrangeano, dado pela Equação (A16), na qual 𝜆(𝑡) é o vetor das variáveis adjuntas de dimensão 𝑛. 𝐿(𝑥, 𝑢, 𝑡) = 𝐹(𝑥, 𝑢, 𝑡) + 𝜆𝑇(𝑓 − �̇�) (A16) O funcional aumentado (𝐽𝐴), se não considerado em sua definição o termo relacionado ao estado final ℎ, é dado por: 𝐽𝐴 = ∫ 𝐿(𝑥, �̇�, 𝑢, 𝜆, 𝑡) 𝑡𝑓 𝑡0 𝑑𝑡 (A17) A mecânica lagrangeana em geral é baseada em equações diferenciais que descrevem o movimento de qualquer corpo e que dependem das coordenadas generalizadas 𝑞 = (𝑞1, … , 𝑞𝑛) do corpo e suas velocidades �̇� = (�̇�1, … , �̇�2). Essas equações são construídas a partir de uma função lagrangeana 𝐿 = 𝑇 − 𝑉 (onde: 𝑇 é a energia cinética do sistema e 𝑉 o potencial). Então, para maximizar 𝐽𝐴 devem ser cumpridas as equações de Euler-Lagrange dadas pela Equação (A11) (MONROY PEREZ, 2016; RUZ, 2014; BOSCAIN; PICCOLLI, 2005). A.1.2 Princípio do Máximo de Pontryagin (PMP) Utilizando-se o PMP pode-se transformar o problema de controle ótimo em um problema de dois pontos de valores no contorno pela introdução de 𝑛 variáveis adjuntas adicionais. O problema, então, consiste em resolver dois conjuntos de equações diferenciais, sendo um deles referente às equações dinâmicas do modelo matemático e o outro referente às equações diferenciais das variáveis adjuntas, vinculados por meio do operador hamiltoniano H do sistema, o qual é posteriormente otimizado. Para o primeiro conjunto são conhecidas apenas as condições iniciais e para o outro, são fixadas as condições finais (BAYEN et al., 2015; BARCIOG et al. 2012, MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). Uma dedução detalhada do equacionamento resultante pode ser encontrada em Ramirez (1994). 53 Seja um processo com 𝑛 variáveis de estado 𝑥(𝑡), e variáveis de controle 𝑢(𝑡) e comportamento dinâmico descrito por: �̇� = 𝑓(𝑥) + 𝑔(𝑥)𝑢 (A18) Para este processo podem ser especificadas as condições iniciais ou finais. As variáveis de controle estão restritas por: 𝑢𝑚𝑖𝑛 ≤ 𝑢(𝑡) ≤ 𝑢𝑚𝑎𝑥 (A19) Deseja-se encontrar a trajetória das variáveis de controle ao longo do tempo que minimiza (ou maximiza) um dado critério de desempenho, que é uma função escalar 𝜙 dos valores dos estados no tempo final 𝑡𝑓, isto é: 𝐽 = 𝜙(𝑥(𝑡𝑓)) (A20) Na Equação (A20), o índice do desempenho não tem segundo termo como na Equação (A15), tendo neste caso somente o custo terminal. As condições necessárias para a otimização, na presença de restrições sobre o controle contínuo 𝑢(𝑡), são convenientemente expressas em termos do Hamiltoniano 𝐻 (BAYEN et al., 2015; BARCIOG et al. 2012, MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). A mecânica hamiltoniana é baseada em equações diferenciais de movimento escritas em termos de uma função denominada de hamiltoniano. Ao contrário das equações de Euler-Lagrange (Equação (A11)), essas equações diferenciais são de primeira ordem, o que é alcançado pela substituição de variáveis de velocidade generalizadas por outras variáveis conhecidas como momentos conjugados. Uma dedução detalhada do equacionamento resultante pode ser encontrada em Rodriguez-Alfaro (2015). Partindo-se da Equação (A11) e utilizando a transformada de Legendre, define-se a função hamiltoniana ou o Hamiltoniano da seguinte maneira: 𝐻 (𝑞, 𝑝, 𝑡) = �̇�𝑇𝑝 − 𝐿 (𝑞, �̇�, 𝑡) (A21) Interpretando mais facilmente, o Hamiltoniano é igual à soma do integrando da integral do custo da Equação (A15) e do produto do vetor transposto de 𝜆 com o membro direito da Equação (A18) (SMETS et al. 2002). 𝐻 = 𝐹(𝑥, 𝑢, 𝑡) + 𝜆𝑇 [(𝑓(𝑥) + 𝑔(𝑥)) 𝑢] (A22) 54 Então, a função hamiltoniana para o sistema dinâmico dado genericamente pela Equação (A18) e que minimiza/maximiza o critério de desempenho dado pela Equação (A20) é: 𝐻 = 𝜆𝑇 [(𝑓(𝑥) + 𝑔(𝑥)) 𝑢] (A23) Na Equação (A23), como definido anteriormente, 𝜆 é um vetor de variáveis adjuntas de dimensão 𝑛 (BAYEN et al., 2015; BARCIOG et al. 2012, MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). Segundo o PMP, para que o controle contínuo u(t), aplicado no intervalo 𝑡0 < 𝑡 < 𝑡𝑓, produza um mínimo/máximo absoluto de um determinado critério de desempenho (𝐽) e uma resposta ótima 𝑥∗(𝑡) para o sistema dinâmico dado pela Equação (A18), sujeito à condições dadas em 𝑥(𝑡0) ou 𝑥(𝑡𝑓), é necessário que exista um vetor de variáveis adjuntas 𝜆(𝑡), de dimensão n, que satisfaça as seguintes equações: �̇� = 𝑑𝑥 𝑑𝑡 = ( 𝜕𝐻 𝜕𝜆 ) (A24) �̇� = 𝑑𝜆 𝑑𝑡 = − ( 𝜕𝐻 𝜕𝑥 ) (A25) As equações dadas pela Equação (A24) são relações de restrição enquanto que as dadas pela Equação (A25) são as equações de Euler-Lagrange. Para prosseguimento do estudo, será considerada apenas uma única variável de controle u e o Hamiltoniano será reescrito da seguinte forma: 𝐻 = 𝐻0(𝑡) + 𝜙(𝑡)𝑢 (A26) Comparando-se a Equação (A26) com a Equação (A23), tem-se que: 𝐻0(𝑡) = 𝜆𝑇𝑓(𝑥) (A27) 𝜙(𝑡) = 𝜆𝑇𝑔(𝑥) (A28) Para os chamados problemas de ponto terminal livre, as condições dadas para o sistema dado pela Equação (A24) são as condições iniciais (𝑥(𝑡0)), enquanto que as condições finais para o sistema representado pela Equação (A25) são especificadas por 𝜆𝑇(𝑡𝑓) = [ 𝜕𝜙 𝜕𝑥 ] 𝑥=𝑥∗(𝑡𝑓) (condição de transversalidade). Para os 55 problemas de ponto terminal fixo, se são dadas todas as condições finais do sistema de equações dado pela Equação (A24) (𝑥(𝑡𝑓)), especificam-se as condições iniciais para o sistema de equações dado pela Equação (A25) como 𝜆(𝑡0) = 0. Estes dois sistemas de equações constituem o problema de dois pontos de valores no contorno (BAYEN et al., 2015; BARCIOG et al. 2012, MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). Para que o ponto ótimo seja atingido, o PMP requer que o Hamiltoniano assuma um valor máximo absoluto com respeito à variável de controle. Assim, as seguintes equações devem ser satisfeitas: 𝜕𝐻 𝜕𝑢 = 0 𝑝𝑎𝑟𝑎 𝑢𝑚𝑖𝑛 < 𝑢 < 𝑢𝑚𝑎𝑥 , 𝑡𝑜 < 𝑡 < 𝑡𝑓 𝐻 = 𝐻𝑚í𝑛 𝑜𝑢 𝐻𝑚á𝑥 𝑝𝑎𝑟𝑎 𝑢 = 𝑢𝑚í𝑛 𝑜𝑢 𝑢 = 𝑢𝑚á𝑥, 𝑡𝑜 < 𝑡 < 𝑡𝑓 (A29) É possível obter mais informações do PMP derivando o Hamiltoniano com relação ao tempo, obtendo-se: 𝑑𝐻 𝑑𝑡 = ( 𝜕𝐻 𝜕𝑥 ) 𝑇 𝑑𝑥 𝑑𝑡 + ( 𝜕𝐻 𝜕𝑢 ) 𝑇 𝑑𝑢 𝑑𝑡 + ( 𝜕𝐻 𝜕𝜆 ) 𝑇 𝑑𝜆 𝑑𝑡 (A30) Substituindo na expressão da derivada as Equações (A30) e (A25), chega-se a: 𝑑𝐻 𝑑𝑡 = ( 𝜕𝐻 𝜕𝑢 ) 𝑇 𝑑𝑢 𝑑𝑡 (A31) Caso a variável de controle não esteja sobre as restrições, a Equação (A29) é válida e a Equação (A31) transforma-se em: 𝑑𝐻 𝑑𝑡 = 0 (A32) Se a variável de controle estiver situada nas restrições, seu valor é constante e, logo, 𝑑𝑢 𝑑𝑡 = 0 e, assim, a Equação (A32) é também obtida. Conclui-se, então, que, na resposta ótima, 𝐻 = 𝐻∗ = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒, podendo o valor desta constante ser relacionado às variáveis de estado finais (BAYEN et al., 2015; BARCIOG et al. 2012, MONTAÑO, 2010; COSTA, 1996, PALANKI et al., 1993). 56 Nos problemas de tempo final livre, ou seja, aqueles nos quais o tempo final não é fixado a priori, o PMP estabelece, ainda, que, no tempo final 𝐻(𝑡𝑓) = 0. Nos problemas de tempo final livre, o Hamiltoniano deve ser constante ao longo do caminho ótimo, isto é, 𝐻 = 0 para 𝑡𝑜 < 𝑡 < 𝑡𝑓. Como para os sistemas descritos pelas Equações (A24) e (A25), a variável de controle aparece linearmente no Hamiltoniano, torna-se fácil determinar o valor que maximiza esta função pela análise do sinal de 𝜙(𝑡) na Equação (A26), isto é: 𝜙(𝑡) > 0 ⇒ 𝑢 = 𝑢𝑚𝑎𝑥 (A33) 𝜙(𝑡) < 0 ⇒ 𝑢 = 𝑢𝑚𝑖𝑛 (A34) No caso em que 𝜙(𝑡) = 0 durante um intervalo de tempo (𝑡1, 𝑡2), o PMP não consegue prever a solução do problema e a teoria de controle singular deve ser aplicada, sendo este particular inte