Otimização de Processos, Notas de estudo de Engenharia Química

Baixe Otimização de Processos e outras Notas de estudo em PDF para Engenharia Química, somente na Docsity! Escola Politécnica UFBA OTIMIZAÇÃO DE PROCESSOS Prof. Dr. Ricardo de Araújo Kalid – kalid@ufba.br Programa de Engenharia Industrial da Escola Politécnica da UFBA Ricardo Kalid - Otimização de Processos – kalid@ufba.br 2 CURRÍCULO DO INSTRUTOR ÁREAS DE ATUAÇÃO E LINHAS DE PESQUISA Modelagem e simulação em regime estacionário e transiente de processos Identificação de processos Controle de processos Otimização de processos Síntese de redes de transferência de calor e massa OUTROS Professor do Mestrado em Engenharia Química da UFBA e do Mestrado em Produção Limpa Professor (anos 92 e 93) do Curso de Especialização em Instrumentação e Controle (CEINST) promovido pelo Departamento de Engenharia Mecânica da UFBA Professor de Cursos de Educação Continuada (Controle Avançado, Controle Preditivo Multivariável, Identificação de Processos, Otimização de Processos Químicos, Controle de Colunas de Destilação) para DOW, PETROBRAS, GRIFFIN, EDN, CIQUINE, OXITENO, COPENE. Professor (98) do Curso de Especialização em Automação de Sistemas Industriais (CEASI) promovido pelo Depto de Engenharia Elétrica da UFBA Professor e Coordenador (99) do Curso de Especialização em Controle e Automação de Processos Industriais (CECAPI) promovido pelos Depto de Engenharia Química e Elétrica da UFBA Professor e Coordenador (2000 a 2002) do Curso de Especialização em Instrumentação, Automação, Controle e Otimização de Processos Contínuos (CICOP 1ª e 2ª turmas) promovido pelo Depto de Engenharia Química e UFBA e AINST. Coordenador do II e do III Seminário Nacional de Controle e Automação (II SNCA-2001 e III SNCA-2003) PROJETOS COOPERATIVOS E/OU CONSULTORIAS PARA INDÚSTRIAS MONSANTO-GRIFFIN-POLITENO: síntese de redes de transferência de calor e massa DETEN: simulação do reator radial para desidrogenação de parafinas EDN: participou da equipe de desenvolvimento do plano diretor de automação BRASKEM-UNIB: identificação de processos, sintonia de controladores industriais, simulação, controle e otimização do conversor de acetileno da ETENO II (em andamento) PDAI-BA - Programa de Desenvolvimento da Automação Industrial, participantes: UFBA, UNIFACS, CEFET-BA, CETIND-SENAI, FIEB, SEPLANTEC, PETROBRAS, NITROCARBONO, DETEN, OXITENO, OPP, POLIBRASIL, POLITENO, BRASKEM-UNIB GRIFFIN: Sistema de controle de pH. Modelagem e Otimização do Reator de DCA BRASKEM-UNIB-POLITENO-UFBA: Diagnóstico de Malhas de Controle Preditivo Multivariável (MPC) BRASKEM-UNIB -UFBA: projeto de produção + limpa para minimização/reuso de águas industriais INDICADORES DE PRODUÇÃO CIENTÍFICA Trabalhos apresentados em congressos ou seminários: 12 Trabalhos publicados em periódicos: 2 Dissertação de mestrado (1) e tese de doutorado (1) defendidas e aprovadas: 2 Participação de bancas de mestrado (5) e de doutorado (1): 6 Orientação de Iniciação Científica e Tecnológica 17 (concluídas) e 3 em andamento Orientação de Dissertações de Mestrado: 8 (em andamento), 3 concluídas Ricardo de Araújo Kalid, D. Sc. 04/09/64 kalid@ufba.br (0xx71) 203.9811 / 9984.3316 Prof. Depto Engenharia Química da UFBA Graduação em Engenharia Química – UFBA (88) Mestrado em Engenharia Química - UFBA (91) Doutorado em Engenharia Química – USP (99) Ricardo Kalid - Otimização de Processos – kalid@ufba.br 5 INSTRUTOR: RICARDO KALID - kalid@ufba.br Cursos e apostilas sobre MODELAGEM DE PROCESSOS 1. Operações Unitárias em Regime Transiente – Balanços de Massa, Energia e Momentum Aplicados a Processos Químicos. 2. Identificação de Processos Químicos. SIMULAÇÃO DE PROCESSOS 3. Métodos Numéricos e Simulação de Processos. 4. Programação em MATLAB com Aplicação em Reatores Químicos. CONTROLE DE PROCESSOS 5. Sistemas de Controle dos Principais Equipamentos da Indústria de Processos Químicos e Petroquímicos. 6. Controle de Processos Químicos. 7. Definição da Estrutura do Sistema de Controle Multimalha de Processos Multivariáveis. 8. Controle Avançado de Processos – Estratégias Clássicas de Controle. 9. Controle de Coluna de Destilação. 10. Controle Preditivo Multivariável: DMC - Controle por Matriz Dinâmica. 11. Sintonia Ótima de Controladores Industriais OTIMIZAÇÃO DE PROCESSOS 12. Otimização de Processos Químicos sem restrições 13. Otimização de Processos Químicos com restrições 14. Otimização de Processos Químicos a batelada Ricardo Kalid - Otimização de Processos – kalid@ufba.br 6 OTIMIZAÇÃO DE PROCESSOS Índice Lista de Abreviaturas ............................................................................................ 6 Nomenclatura ........................................................................................................ 7 Principais Referências Bibliográficas ................................................................. 7 1. Introdução e Definições .................................................................................... 12 1.1. Objetivos deste Curso ................................................................................... 12 1.2. Programa do Curso ....................................................................................... 13 1.3. Referências Bibliográficas Principais ........................................................... 14 1.4. Por que Otimizar? ......................................................................................... 14 1.5. Exemplos de Aplicação de Otimização ........................................................ 15 1.6. Formulação de um Problema de Otimização ................................................ 15 1.6.1. A Função Objetivo (FO) .......................................................................... 16 1.6.2. As Restrições ........................................................................................... 18 1.6.3. A Região Viável ...................................................................................... 19 1.6.4. As Variáveis de Decisão (VD) ................................................................ 20 1.7. Procedimento Geral para Solucionar um Problema de Otimização ............. 20 1.7.1. Mapeamento da Função Objetivo ............................................................ 21 1.7.2. Obstáculos à Otimização ......................................................................... 25 1.8. Exercícios ...................................................................................................... 27 2. Conceitos Matemáticos ..................................................................................... 28 2.1. Definições ..................................................................................................... 28 2.2. Operações Básicas com Matrizes e Vetores ................................................. 29 2.3. Independência Linear, Matriz Singular e Rank ou Posto de uma Matriz ... 32 2.4. Operadores Linha ou Coluna ........................................................................ 33 2.5. Solução de Sistema de Equações Lineares ................................................... 33 2.6. Graus de Liberdade ....................................................................................... 34 2.7. Autovalores e Autovetores ............................................................................ 34 2.8. Estudo de Função .......................................................................................... 35 Ricardo Kalid - Otimização de Processos – kalid@ufba.br 7 2.9. Continuidade de Funções .............................................................................. 37 2.10. Funções Unimodais e Multinodais ............................................................... 38 2.11. Funções Côncavas e Convexas ..................................................................... 38 2.12. Região Convexa ............................................................................................ 41 2.13. Condições Necessárias e Condições Suficientes para um Extremo de uma Função Irrestrita ....................................................................................................... 44 2.14. Interpretação da Função Objetivo em Termos de uma Aproximação Quadrática 47 2.15. Exercícios ...................................................................................................... 51 3. Formulação Matemática de um Problema de Otimização ........................ 54 3.1. A Função Objetivo (FO) ............................................................................... 55 3.1.1. Tolerância ou Critério de Parada ............................................................. 56 3.1.2. Objetivos Econômicos ............................................................................. 57 3.1.3. Objetivos Operacionais ........................................................................... 64 3.1.4. Combinação de Objetivos Operacionais com Objetivos Econômicos .... 69 3.2. As Funções de Restrição (FR) ...................................................................... 69 3.3. Otimização On-Line ..................................................................................... 70 3.4. Exercícios ...................................................................................................... 72 4. Otimização Unidimensional Sem Restrições (OUSR) ............................... 79 4.1. Métodos Indiretos (MI) para OUSR ............................................................. 80 4.1.1. Método de Newton .................................................................................. 81 4.1.2. Método de Quasi-Newton ........................................................................ 83 4.1.3. Método da Secante .................................................................................. 84 4.2. Métodos Diretos (MD) para OUSR .............................................................. 85 4.2.1. Métodos por Diminuição da Região de Busca ........................................ 85 4.2.2. Métodos por Aproximação Polinomial - Interpolação Quadrática.......... 87 4.2.3. Métodos por aproximação polinomial - Interpolação cúbica .................. 88 4.3. Avaliação dos Métodos Unidimencionais de Otimização ............................ 89 4.4. Exercícios ...................................................................................................... 91 5. Otimização Multidimensional Sem Restrições (OMSR) ............................ 92 5.1. Métodos Indiretos (MI) para OMSR ............................................................ 93 5.1.1. Método do Gradiente ou Método do Gradiente Descendente ................. 94 5.1.2. Método do Gradiente Conjugado ............................................................ 96 5.1.3. Método de Newton .................................................................................. 97 5.1.4. Método de Levenberg-Marquardt ............................................................ 98 Ricardo Kalid - Otimização de Processos – kalid@ufba.br 10 Lista de Abreviaturas CKT - Condições de Kuhn-Tucker CO - Controle Ótimo FO - Função Objetivo FR - Funções de Restrições MD - Método Direto de otimização MI - Método Indireto de otimização MV - Variável Manipulada OMCR - Otimização Multidimensional Com Restrições OMSR - Otimização Multidimensional Sem Restrições OUCR - Otimização Unidimensional Com Restrições OUSR - Otimização Unidimensional Sem Restrições PCO - Problema de Controle Ótimo PL - Programação Linear PNL - Programação Não-Linear PPL - Problema de Programação Linear PPNL - Problema de Programação Não-Linear PPQ - Problema de Programação Quadrática PQ - Programação Quadrática PV - Variável de Processo (pode ser variáveis controladas e/ou medidas) s.a. - sujeito a SEANL - Sistema de Equações Algébricas Não-Lineares SP - SetPoint VA - Variável Auxiliar (qualquer variável ou constante que não é VDep) VD - Variável de Decisão ou de Projeto ou Independente VDep - Variável Dependente VI - Variável Independente Ricardo Kalid - Otimização de Processos – kalid@ufba.br 11 Nomenclatura Y i - valor da média dos experimentos replicados para um certo xi x - valor médio da variável independente ηi - valor calculado da VI para um certo xi bj - parâmetro j de um modelo n - número total de experimentos realizados pi - número de medidas replicadas para um certo xi x - variável determinística, independente e experimental Y - variável estocástica, dependente e experimental Y - variável estocástica, dependente e experimental Principais Referências Bibliográficas As principais referências utilizadas para confeccionar esta apóstila e preparar as aulas são listadas a seguir. Outras referências, que tratam de assuntos específicos, são citadas ao longo do texto: R1. Himmelblau, D. M. and Edgar, T. F.; Optimization of Chemical Process. McGraw-Hill, 1989. Livro essencial para quem quer iniciar e/ou aprofundar seus estudos sobre otimização de processos químicos. Boa parte do conteúdo deste curso e a maioria dos exercícios discutidos/propostos foram retirados deste livro. R2. Himmelblau, D. M.; Process Analysis by Statistical Methods. Jonh Wiley & Sons, 1970. Livro que traz os algoritmos de vários métodos de otimização e aplica esses métodos principalmente ao ajuste de modelos matemáticos a dados experimentais. Livro texto para o Capítulo 6 desta apostila. R3. Beveridge, G. S. and Schehter, R. S.; Optimization Theory and Practice. McGraw-Hill, 1970. Traz uma discussão mais profunda a respeito dos fundamentos matemáticos em que os métodos de otimização são baseados. R4. Reklaitis, G. V.; Ravindran, A.; Ragsdell, K. M.; Engineering Optimization: Methods and Applications. Jonh Wiley & Sons, 1983. Livro importante e complementar ao de Himmelblau e Edgar (R1). Livro texto para os Capítulo 8 e 9 desta apostila. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 12 1. Introdução e Definições Na natureza somente os mais eficientes sobrevivem. A seleção natural elimina as espécies ou os indivíduos menos capacitados, ao mesmo tempo que facilita as coisas para os mais bem dotados, ou seja, existe algum mecanismo que procura maximizar o uso dos recursos ou minimizar o efeito da ineficiência dos processos. O que a natureza sabe fazer com extrema maestria nós devemos imitar, se quisermos suplantar nossos adversários e assim sermos os ganhadores. O problema é que seja qual for o vencedor, a mãe natureza sempre ganha junto, enquanto nós temos que ser mais eficientes que nossos concorrentes. E como a probabilidade de tomar uma decisão errada é muito maior que escolher a opção certa, temos que nos cercar de informações e procedimentos confiáveis. As informações que devemos coletar são, por exemplo: • o(s) objetivo(s) que queremos atingir, • as condições do mercado fornedor e consumidor, • as condições do mercado financeiro, • as condições dos recursos naturais e humandos disponíveis, • as limitações de natureza física e/ou social e/ou psicológicas existentes. Os procedimentos que podemos empregar são: • a experiência acumulada, • o uso de ferramentas matemáticas adequadas, • a escolha de estratégias de atuação oportunas. Se utilizarmos as informações e procedimentos apropriados temos boas chances de sermos os vitoriosos. É somente esse o nosso problema. Neste curso iremos estudar quais são as informações e procedimentos adequados a resolver os os problemas de otimização típicos em engenharia química. 1.1. Objetivos deste Curso Ao final do segundo módulo deste curso seremos capazes de: O1. Entender os princípios de funcionamento dos pricipais algoritmos de otimização. O2. Aplicar corretamente os algoritmos de otimização. O3. Desenvolver um modelo matematico de otimização e resolvê-lo através da utilização de pacotes computacionais. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 15 Em certos casos o uso de métodos matemáticos de otimização não traz benefícios, por exemplo: • quando uma solução razoável pode ser obtida baseada em experiências passadas, • quando existe uma prática padrão em uso, • quando o tempo necessário para avaliar o problema não contribui para o contexto do projeto global, • quando as informações necessárias só podem ser obtidas com grande custo. 1.5. Exemplos de Aplicação de Otimização A otimização pode ser aplicada de inúmeras maneiras em processos e plantas químicas. Típicos projetos onde a otimização tem sido empregada incluem: 1. Determinação do melhor local para construção de uma planta. 2. Escalonamento de tanques para armazenagem de matéria-prima e de produtos. 3. Dimensionamento e layout de pipelines. 4. Projeto de plantas e/ou de equipamentos. 5. Escalonamento de reposição e manutenção de equipamentos. 6. Operação de equipamentos e/ou plantas. 7. Ajuste de modelos a dados experimentais de uma planta. 8. Minimização de inventário. 9. Alocação de recursos ou serviços entre diferentes processos. 10.Planejamento e escalonamento de instalação de plantas. 1.6. Formulação de um Problema de Otimização Podemos definir# a otimização de sistemas das seguintes maneiras: D1. Campo da matemática dedicado ao desenvolvimento de métodos eficientes de determinação de máximos e mínimos de funções de uma ou mais variáveis. D2. A ciência que determina as melhores soluções para certos problemas físicos; problemas que são descritos por modelos matemáticos. # As três definições foram transcritas da apóstila sobre "Otimização de Processos" de autoria de Fernando Pellegrini Pessoa e Marcelo Castier, ambos professores do DEQ-UFRJ Ricardo Kalid - Otimização de Processos – kalid@ufba.br 16 D3. Busca da melhor solução (solução ótima) dentre as diversas soluções possíveis de um problema, segundo um critério estabelecido previamente. Embora nem todos os problemas de otimização possam ser descritos por equações matemáticas (por exemplo, maximizar as propriedades organolépticas de um alimento), apenas os que podem serão tratados neste curso. A formulação de um problema prático de otimização contém duas partes: 1. Ao menos uma função objetivo a ser alcançada. 2. As restrições que devem ser atendidas. Devemos encontrar uma solução (pois podem existir mais) que minimiza/maximiniza a função objetivo e que simultaneamente atenda às restrições, ou seja a solução encontrada deve pertecer à região viável. Para que as condições ótimas sejam alcançadas o sistema deve ter liberdade para manipular as variáveis de decisão, também denominadas variáveis de projeto ou ainda variáveis independentes, isto é, algumas condições operacionais são modificadas de forma que o ponto ótimo viável seja alcançado. 1.6.1. A Função Objetivo (FO) A função objetivo ou critério de desempenho estabelece o alvo a ser alcançado. É uma função matemática cujo máximo ou mínimo se deseja determinar. As FO's podem ser desenvolvidas a partir de três tipos de critérios: C1. Critério estritamente econômico: maximizar o lucro anual, minimizar o custo anual, diminuir o tempo de retorno do investimento, etc. C2. Critério estritamente técnico/operacional: minimizar o consumo de energia ou de matéria-prima, maximizar a produção, etc. C3. Critério técnico-econômico: minimizar a diferença entre o valor desejado e o valor medido numa planta, ao mesmo tempo em que minimiza o custo operacional. O estabelecimento correto da função objetivo é fundamental para o sucesso da otimização. Sua determinação é uma tarefa complexa que requer grande conhecimento do processo/sistema a ser otimizado. A função objetivo pode ser classificada quanto à: • Continuidade: • contínua, por exemplo a temperatura ótima para uma reação reversível • descontínua • discreta, por exemplo o diâmetro ótimo de uma tubulação Ricardo Kalid - Otimização de Processos – kalid@ufba.br 17 • Modalidade: • unimodal (o extremo local é também o global), figura 1-a. x = -100:100; a = 4; b = 2; c = 5; y1= -a*x.^2 + b*x + c; plot(x,y1) -100 -50 0 50 100 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 x 104 Figura 1-A: Função Unimodal Função objetivo (FO): max y y a x b x c 1 1 2= − + +. . Equação 1.6-A max(y1) ans = 5 Obs: expressões nesta fonte (Courier New 10, verde, em negrito) são comandos do MATLAB. expressões nesta fonte (Courier New 10, azul) são as respostas geradas pelo MATLAB. as expressões nesta fonte (Times New Roman 12, preto) são textos em WORD. Comando help nome-da-função permite consultar o manual de referência on-line do MATLAB. Comando help optim permite consular o manual de referência do toolbox de otimização on-line do MATLAB. • multimodal (existem vários extremos locais e um deles é o global) Ricardo Kalid - Otimização de Processos – kalid@ufba.br 20 atenção deve ser dispensada ao estabelecimento das restrições, ou seja, da região viável. 1.6.4. As Variáveis de Decisão (VD) As variáveis de decisão ou de projeto ou independentes correspondem, em número, ao excesso de incógnitas em relação ao número de equações, ou seja, sua quantidade é igual ao número de graus de liberdade do sistema. Se existe apenas uma única solução para o problema, nenhuma otimização é necessária e possível. Portanto, para haver condições de otimizar um processo o mesmo dever ter graus de liberdade maior que zero. As VD caracterizam os possíveis projetos ou condições operacionais do sistema e devem ter uma certa influência sobre a função objetivo. Se uma função objetivo é pouco sensível a uma variável de decisão é conveniente simplificar o problema assumindo um valor fixo para essa variável. Por outro lado, se o critério de desempenho é extremamente sensível a uma determinada variável de projeto, talvez seja difícil reproduzir na prática as condições ótimas calculadas. 1.7. Procedimento Geral para Solucionar um Problema de Otimização Não existe um procedimento ou método que possa ser aplicado eficientemente para todo tipo de problema. A escolha do método depende: • da característica da função objetivo (linear ou não-linear; contínua, discreta ou mista) • da natureza das restrições (linear ou não-linear; contínua, discreta ou mista) • do número de variáveis de decisão. Podemos estabelecer 6 passos principais a serem seguindos na solução de problemas de otimização: P1. Analise o processo e estabeleça as variáveis de decisão (VD) e as auxiliares (VA’s). P2. Estabeleça a função objetivo (FO) em função das variáveis identificadas no item P1 e de coeficientes conhecidos. P3. Estabeleça as restrições: • balanços de massa e energia funções das VD’s e VA’s, • relações constitutivas e/ou empíricas, • limites operacionais máximos e mínimos, • faixa de validades das variáveis, por exemplo, temperaturas e pressões absolutas devem ser positivas, frações molares entre 0 e 1, etc., Ricardo Kalid - Otimização de Processos – kalid@ufba.br 21 • limites externos, por exemplo, capacidade máxima de consumo do mercado. P4. Se o problema é demasiadamente grande: • subdivida em partes e/ou • simplifique a função objetivo e/ou o modelo do processo. P5. Se possível faça o mapeamento da função objetivo, isto é, verifique graficamente como a FO varia com a mudança das variáveis de decisão. P6. Aplique as apropriadas técnicas matemáticas de otimização para o problema. P7. Aplique a análise de sensibilidade da FO, isto é, examine a sensilidade do ponto de mínimo/máximo e o valor da FO às mudanças nos coeficientes das funções e a alterações nas variáveis decisão. Os passos P1, P2, P3 e P4 constituem a representação matemática do problema, exigindo que a equipe responsável pela otimização do sistema tenha: (a) muito conhecimento a respeito do processo, (b) muita atenção e habilidades específicas em modelagem de processos. Deve-se escolher um modelo o mais simples possível (menor quantidade de variáveis, equações constitutivas enxutas) que representa adequadamente o sistema. Executando a etapa P5 temos condições de saber sua ordem de grandeza e de verificar como ela varia com determinada VD. O mapeamento da FO nem sempre é simples de realizar embora seja sempre desejável. A aplicação dos algoritmos de otimização (passo P6) é uma etapa simples, desde que se tenha à mão programas de computador já desenvolvidos e testados. Porém se for necessário implementar ou desenvolver um novo algoritmo será necessário um grande esforço (tempo e recursos humanos). Na etapa P7, validação dos resultados obtidos, é essencial a participação de engenherios e de técnicos que conhecem (bem) o sistema/processo otimizado. 1.7.1. Mapeamento da Função Objetivo Antes de executar algum algoritmo de otimização, é interessante que seja realizado o mapeamento da FO. Veja o exemplo a seguir de uma função Z de duas incógnitas (X,Y): x = -8:.5:8 ; y = x ; [X,Y] = meshgrid(x,y) ; R = sqrt(X.^2 + Y.^2) + eps ; Z = sin(R)./R ; Utilizamos as funções de traçamento de gráficos do MATLAB para mapear de diversas maneiras a função Z = seno(R)/R. Construimos os gráficos bidimensional (figura 1-d), tridimensional (figura 1-e), das curvas de nível Ricardo Kalid - Otimização de Processos – kalid@ufba.br 22 (figura 1-f), combinação do gráfico tridimensional com o das curvas de nível (figura 1-g), curvas de nível pseudo-coloridas (figura 1-h). Ricardo Kalid - Otimização de Processos – kalid@ufba.br 25 (b) plotar um gráfico bidimensional no qual a ordenada é o valor da FO e na abscissa os conjuntos de valores das variáveis de projeto, Figura 1-I. Por exemplo [nl,nc] = size(Z) ; VetorZ = [ ] ; for i = 1 : nl for j = 1 : nc VetorZ = [ VetorZ ; Z(i,j) ] ; end end whitebg plot(VetorZ,'r') 0 2 0 0 4 0 0 6 0 0 8 0 0 1 0 0 0 1 2 0 0 -0 .4 -0 .2 0 0 .2 0 .4 0 .6 0 .8 1 Figura 1-I: Gráfico dos valores alinhados de uma função bidimensional Quando existem muitas VD’s o número de cálculos e/ou de gráficos para mapear a FO fica proibitivo, por exemplo, para 5 variáveis de decisão com 20 pontos para cada, o número de vezes que a FO será avaliada é 205 = 3.200.000. 1.7.2. Obstáculos à Otimização Se a função objetivo (FO) e as funções de restrições (FR) forem "bem comportadas" a otimização não apresenta grandes problemas. Particularmente se a FO e as restrições forem todas lineares está disponível um poderoso método (Programação Linear) que resolve este problema de maneira satisfatória. Entretanto muitos problemas de otimização são não-lineares. Muitos dos problemas práticos de otimização em engenharia química apresentam alguns das dificuldades abaixo: P1. Não disponibilidade de dados ou de um modelo matemático confiável do sistema. P2. Descontinuidades da FO ou das FR’s. Por exemplo, o preço de um compressor ou trocador de calor não varia continuamente como uma função das dimensões, pressão, temperatura, etc. pois o incremento de um Ricardo Kalid - Otimização de Processos – kalid@ufba.br 26 parâmetro não afeta o custo em uma certa faixa de variação, porém para uma outra faixa ocorre um salto no valor do equipamento. P3. Não-lineraridade da FO ou FR. P4. A FO e a FR são definidas através expressões que contém complicadas interações entre as variáveis de decisão. A interação impede a existência de um único ponto ótimo. P5. A FO ou a FR tem comportamento achatado ou exponencial em algumas faixas de variação das variáveis de decisão. Isto significa que o problema é pouco ou extremamente sensível, respectivamente, a mudanças dessas variáveis. P6. A FO apresenta muitos extremos locais perto da região que contém o extremo global. Quando não é possível a aplicação dos métodos de otimização podemos: (a) fazer um estudo de caso, isto é, escolher criteriosamente um número limitado de opções e analisar qual a melhor alternativa; (b) fazer um estudo de sensibilidade, semelhante ao estudo de caso, apenas mais sistematizado e com um número maior de casos a analisar. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 27 1.8. Exercícios E1.1. Crie uma nova função de duas VD’s acrescentando à função Z do ítem 1.7.1 (página 21) um ruido randômico (a função no MATLAB é rand). Construa gráficos bidimensionais, tridimensionais e de contorno. Analise os resultados obtidos assumindo que a nova função Z é uma FO. E1.2. Desenvolva uma função objetivo para otimizar seus rendimentos. Não esqueça de estabelecer as restrições impostas pelo sistema (patrão, esposa(o), filhos, etc). E1.3. Dada a função ( ) 222121 5.123 xxxxxf ++= construa seu gráficos e ache o ponto de mínimo. E1.4. Dada a função ( ) ( )[ ]22121231 10exp xxxxxxf −−−= , graficamente obtenha o ponto de máximo. Estude a sensibilidade. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 30 Multiplicação de uma matriz por um escalar: Cada elemento da matriz é multiplicada pelo número escalar: ijij basBAs =⇔= . Equação 2-G Transposta de um produto de matrizes: ( ) TTT ABBA = Equação 2-H Produto interno entre dois vetores: Seja x e y dois vetores de dimensão n, então ∑ = == n i ii T yxyxyx 1 , Equação 2-I ∑ = = n i i T xxx 1 2 é um escalar Equação 2-J Se o resultado do produto interno entre dois vetores é igual a zero, então esses vetores são ortogonais, em vetores bi ou tridimensionais isto significa que os mesmos são perpendiculares entre si. Inversa de uma matriz: Não existe a versão matricial da divisão escalar. Por definição a inversa de uma matriz, necessariamente quadrada, é a matriz tal que IAAAA == −− 11 Equação 2-K onde I é a matriz identidade Um uso freqüente da inversa de uma matriz é para expressar um conjunto de variáveis em termos de um outro conjunto, uma operação importante em problemas de otimização com restrições. Por exemplo: xzAxAz =⇔= −1 Equação 2-L Ricardo Kalid - Otimização de Processos – kalid@ufba.br 31 Determinante de uma matriz: A função determinante, ou simplismente o determinante de uma matriz (necessariamente quadrada), é denotado por det[A] ou |A|. É um número escalar indispensável na investigação de sistemas de equações lineares. Não iremos apresentar a definição do determinante, mas apenas como calculá-lo: • para uma matriz de dimensão 1, |A1x1| = a11 • para uma matriz de dimensão 2, |A2x2| = a11.a22 - a12.a21 • para uma matriz de dimensão 3, |A3x3| = a11.(a22.a33 - a23.a32) - a12.(a21.a33 - a23.a31) + a13.(a21.a32 - a22.a31) • para uma matriz de dimensão n, utilize o comando do MATLAB: det(A) Derivada de uma função escalar de um campo vetorial: Seja uma função escalar de um vetor de n variáveis f(x) , define-se vetor gradiente ou simplismente gradiente de f(x) ao operador ( ) ( ) ( ) ( ) ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ==∇ Δ nx xf x xf x xfxf ∂ ∂ ∂ ∂ ∂ ∂ 1 Equação 2-M Define-se a matriz Hessiana ou simplesmente Hessiana(o) de f(x) , H(x) , à matriz das derivadas segunda ordem de de f(x) : ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ == Δ 2 2 2 2 1 2 2 2 2 2 2 12 2 1 2 21 2 2 1 2 2 2 nnn n n x xf xx xf xx xf xx xf x xf xx xf xx xf xx xf x xf x xfxH ∂ ∂ ∂∂ ∂ ∂∂ ∂ ∂∂ ∂ ∂ ∂ ∂∂ ∂ ∂∂ ∂ ∂∂ ∂ ∂ ∂ ∂ ∂ Equação 2-N Seja o vetor b , n x 1 , de coeficientes constantes então ( ) ( ) b x xf x xfb b x xb TTT ⎥ ⎦ ⎤ ⎢ ⎣ ⎡∂ =⇔= ∂∂ ∂ ∂ ∂ Equação 2-O Seja os vetores x e z , n x 1 , e a matriz A , n x n , a partir da definição dada pela equação 2-m demonstra-se que: Ricardo Kalid - Otimização de Processos – kalid@ufba.br 32 ( ) zA x zAxT = ∂ ∂ Equação 2-P e ( ) xAxA x xAx T T += ∂ ∂ Equação 2-Q Se A é simétrica vale a seguinte expressão ( ) xA x xAxT 2= ∂ ∂ Equação 2-R e a partir da equação 2-n, se A além de simétrica for real, então ( ) A x xAxT 22 2 = ∂ ∂ Equação 2-S 2.3. Independência Linear, Matriz Singular e Rank ou Posto de uma Matriz Matriz cujo determinante é igual a zero é denominada de matriz singular. Isto acontece quando todos os elementos de uma ou mais linhas (ou colunas) são nulos ou quando uma ou mais linhas (ou colunas) da matriz tem uma dependência linear com outra(s) linha(s) [ou coluna(s)], isto é, podemos escrever uma linha (ou coluna) como uma função linear de uma ou mais linhas (ou colunas) da matriz. Portanto se o determinante de uma matriz é zero, esta matriz é singular e suas linhas ou colunas são linearmente dependentes. Para matrizes quadradas linhas dependentes implicam em colunas dependentes. Por definição as colunas da matriz A, aj, são linearmente independentes se d aj j j n = ∑ = 1 0 se e somente se dj = 0 qualquer que seja o j Equação 2-T A dependência linear ocorre quando para algum valor não nulo de dj a equação 2-t é satisfeita. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 35 vvA α= Equação 2-Z Para cada autovalor corresponde um autovetor. Os autovalores e autovetores fornecem informações a respeito da natureza das funções em uma otimização. Se todos os autovalores de uma matriz A são positivos (maiores que zero), então A é positiva definida e tem inversa. Se todos os autovalores de uma matriz A são negativos (menores que zero), então A é negativa definida. Rearranjando a equação 2-z obtemos ( ) 0=− vIA α Equação 2-AA cujas variáveis desconhecidas são os autovalores α e os autovetores v . Devido o lado direito da equação 2-aa ser zero, duas soluções são possíveis: • a solução é trivial: v = 0 • ou existe mais de uma solução, então det[A - αI] = 0. 2.8. Estudo de Função Seja uma função de uma única variável, f(x) , a condição necessária, mas não suficiente, para que um ponto x*, que não esteja na fronteira do domínio da função, ou simplesmente ponto interno, seja ponto de máximo ou de mínimo é que a derivada da função neste ponto seja nula. Seja o ponto x* o ponto crítico ou ponto estacionário isto é f'(x*) = 0. Este ponto pode ser um ponto de máximo, de mínimo ou de inflexão, veja a figura 2- a. 1500 2000 2500 Figura 2-A: Função monovariável A condição suficiente para que x* seja um mínimo local é que f'(x*) = 0 e f''(x*) > 0, ou seja, se f'(x*) = 0 e f''(x*) > 0 então x* é ponto de mínimo. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 36 Porém existem pontos de mínimo que não atendem a essa condição, daí a condição suficiente. Por exemplo em f(x).= x4/3 , x* = 0 é o ponto de mínimo, embora f''(x*) não seja definido neste ponto. A condição suficiente para que x* seja um máximo local é que f'(x*) = 0 e f''(x*) < 0, ou seja, se f'(x*) = 0 e f''(x*) < 0 então x* é ponto de máximo. Porém existem pontos de máximos que não atendem a essa condição, daí a condição suficiente. A condição necessária para que o ponto x* de uma função multivariável seja estacionário é que todas as derivadas parciais da função em relação às suas variáveis sejam nulas naquele ponto, isto é, f'(x*) = 0 A condição suficiente para que o ponto x* de uma função multivariável seja mínimo local é que a matriz das derivadas parciais de segunda ordem da função em relação a todas as suas variáveis seja positiva definida naquele ponto, isto é, H(x*) > 0. A condição suficiente para que o ponto x* de uma função multivariável seja máximo local é que a matriz das derivadas parciais de segunda ordem da função em relação a todas as suas variáveis seja negativa definida naquele ponto, isto é, H(x*) < 0. Na seção 2.11 será apresentado o conceito de matriz positiva (negativa) definida. Por ora basta entender que existe um paralelo entre esse conceito e o de derivada positiva (negativa) de segunda ordem de uma função monovariável. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 37 2.9. Continuidade de Funções Em problemas de otimização é melhor que as funções e suas derivadas sejam contínuas. Na figura 2-b vemos um gráfico de uma função descontínua, enquanto que na figura 2-c temos uma função contínua, mas com derivada de 1a ordem descontínua. 0 5 10 15 20 25 30 35 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x 104 Figura 2-B: Função descontínua 0 5 10 15 20 25 30 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Figura 2-C: Função contínua com derivada descontínua Uma função monovariável é continua no ponto xo se: a) f(xo) existe b) )(lim xf oxx→ existe c) )()(lim oxx xfxfo = → Uma função de duas variáveis é contínua no ponto (xo,yo) se: a) f(xo,yo) existe b) ( ) ( ) ),(lim ,, yxf oo yxyx → existe c) ( ) ( ) ( ) ( )ooyxyx yxfyxfoo ,,lim ,, =→ Quanto mais afastada estiver a descontinuidade do ponto de máximo (mínimo) da função, mais facilmente este ponto é encontrado. Se uma função não é continuamente diferenciável, veja figura 2-c, métodos de otimização que utilizam derivadas não podem ser empregados, pois a derivada não é definida nos pontos de descontinuidade de f(x). Um exemplo de descontinuidade acontece no projeto de tubulações. Existe um número finito de diâmetros de tudo disponíveis para compra, veja figura 2- d, portanto a escolha do diâmetro recairá entre uma dessas opções. Para resolver esse problema existem duas alternativas: Ricardo Kalid - Otimização de Processos – kalid@ufba.br 40 A equação 2-bb não é uma maneira conveniente de testar a concavidade ou convexidade de uma função. Em seu lugar utiliza-se a segunda derivada de f(x), 2f(x), freqüentemente denominada de matriz de Hessian ou matriz Hessiana de f(x) e denotada por H(x), uma matriz simétrica das derivadas de segunda ordem de f(x). Por exemplo, se f(x) é uma função quadrática de 2a ordem ( ) 222221122111 xhxxhxhxf ++= Equação 2-CC a Hessiana é dada por ( ) ( ) ( ) ( ) ( ) ( ) ⎥⎦ ⎤ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =∇= 2212 1211 2 2 2 12 2 21 2 2 1 2 2 2 2 hh hh x xf xx xf xx xf x xf xfxH ∂ ∂ ∂∂ ∂ ∂∂ ∂ ∂ ∂ Equação 2-DD Analogamente às funções monovariáveis, a matriz das derivadas segundas de uma função multivariável, ou seja a Hessiana, determina se uma função tem um máximo ou um mínimo em uma região, ou seja, a matriz H(x) estabelece se a função é côncava ou convexa numa dada região, respectivamente. Resumindo, dada uma função objetivo f(x) a sua matriz Hessiana H(x) determina a sua natureza: 1. Para f(x) ser estritamente convexa H deve ser positiva definida. Para f(x) ser convexa H deve ser positiva semi-definida. H é positiva definida se e somente se xTHx é > 0 para todo x ≠ 0. H é positiva semi-definida se e somente se xTHx é > 0 para todo x ≠ 0. Se todos os autovalores de H forem positivos (>0) então H é positiva definida. Se todos os autovalores de H forem não-negativos (>0) então H é positiva semi-definida. 2. Para f(x) ser estritamente côncava H deve ser negativa definida. Para f(x) ser côncava H deve ser negativa semi-definida. H é negativa definida se e somente se xTHx é < 0 para todo x ≠ 0. H é negativa semi-definida se e somente se xTHx é < 0 para todo x ≠ 0. Se todos os autovalores de H forem negativos (<0) então H é negativa definida. Se todos os autovalores de H forem não-positivos (<0) então H é negativa semi-definida. 3. H é indefinida se xTHx < 0 para alguns x e xTHx > 0 para outros. Se uma função tem um ponto estacionário onde a matriz Hessiana tem autovalores negativos e positivos, a função não é côncava e nem convexa. 4. Para f(x) ser simultaneamente côncava e convexa H = 0, e neste caso f(x) é uma função linear. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 41 Algumas vezes utiliza-se o fato de que a soma de funções convexas (côncavas) também é uma função convexa (côncava). Ou seja, podemos analisar f(x) em termos das funções que a compõe. Se uma função é separável, isto é, f(x) = f1(x) + f2(x) , e se f1(x) e f2(x) são convexas (côncavas) então f(x) também é convexa (côncava). 2.12. Região Convexa São raros os problemas práticos sem restrições. As restrições podem ser impostas aos valores que as variáveis de decisão podem assumir ou sobre a própria função objetivo. Em geral as restrições são expressas em termos de equações e/ou inequações ( ) ( ) η,,2,10 ,,2,10 … … ++=≤ == kparaxg e kparaxg k k Equação 2-EE É importante o tipo da região de busca definida pelo conjunto de restrições do problema. Regiões convexas são mostradas na figura 2-g (a) e (b). Figura 2-G: Regiões convexas Ricardo Kalid - Otimização de Processos – kalid@ufba.br 42 Na figura 2-h vemos regiões viáveis não-convexas. Figura 2-H: Regiões não-convexas Se uma região é completamente limitada por funções côncavas do tipo gi(x) > 0 , então é uma região convexa fechada. Por outro lado se a região de busca é restrita por inequações da forma gi(x) < 0 e convexas, então é uma região convexa aberta. Encontramos muitas variações quanto ao tipo e número de extremos de uma FO e das suas FR's. No caso geral, é necessário encontrar todos os extremos locais e compará-los para determinar o extremo global. Este procedimento nem sempre é exigido quando certas combinações de FO's e FR's estão presentes. Considere uma função objetivo côncava restrita por uma região convexa. Como a FO é côncava ela admite máximo local que também será o máximo global, pois se a região é convexa, na figura 2-i(a) e (b) observa-se que apenas um ponto de máximo é encontrado. Este, portanto é o máximo local e global. Observe que se a região não for convexa existem vários máximos locais, veja figura 2-i(c), e a comparação entre os mesmos é necessária. Similarmente para uma função convexa numa região convexa, existirá apenas um ponto de mínimo. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 45 valores muito pequenos de Δx ,os termos das derivadas de ordem igual ou superior a 2 são desprezíveis em relação ao elemento das derivadas de 1a. Portanto para x* ser ponto estacionário o gradiente de f(x) deve ser nulo em x* , ∇Tf(x*) = 0 . Essa condição é apenas necessária pois o ponto x* pode ser ponto de sela ou ponto de inflexão. Se a condição necessária (∇Tf(x*) = 0) é atendida, então ∇2f(x*) deve ser maior que zero para que x* seja ponto mínimo, pois a equação 2-gg deve ser satisfeita mesmo para valores muito pequenos de Δx , e nesse caso os termos de ordem igual ou superior a 3 da equação 2-ff são desprezíveis em relação ao termo da derivada 2a. Raciocínio análogo com a equação 2-hh pode ser feito para interpretar se um ponto estacionário é máximo. Logo para x* ser ponto de mínimo, ∇Tf(x*) = 0 e ∇2f(x*) > 0 . Essa condição é apenas suficiente pois pode ser que a 2a derivada não exista no ponto x* , apesar deste ser ponto de mínimo. Similarmente raciocina-se para o ponto de máximo. Se ∇f(x*) = 0 e ∇2f(x*) = 0 então x* é ponto de sela. Resumindo as condições necessárias (CN1 e CN2 abaixo) e as condições suficientes (CS3 e CS4) que garantem que x* é um extremo são as seguintes: CN1. f(x) seja uma vez diferenciável no ponto x* . CN2. ∇f(x*) = 0 , isto é, x* seja um ponto estacionário. CS1. f(x) seja duas vezes diferenciável no ponto x* . CS2. ∇2f(x*) = H(x*) seja positiva definida para que um mínimo exista em x* , e seja negativa definida para que um máximo exista em x* . Na tabela 2-a vê-se um resumo das condições discutidas nesta seção. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 46 Tabela 2-A ∇f(x*) H(x*) = ∇2f(x*) ΔxT.∇2f(x*).Δx Próximo de x* f(x) -f(x) Pode-se afimar que ≠ 0 não é ponto estacionário = 0 positiva definida > 0 > 0 existe mínimo = 0 positiva semi- definida > 0 possivelmente > 0 possivelmente existe mínimo = 0 negativa definida < 0 < 0 existe máximo = 0 negativa semi- definida < 0 possivelmente < 0 possivelmente existe máximo = 0 indefinida ambos > 0 e > 0, dependendo de Δx > 0, < 0 ou nenhum dos dois nada se pode afirmar Observações: • ∇f(x) que é o gradiente da função f(x). • H(x) = ∇2f(x) é denominada de matriz Hessiana de f(x), ou simplesmente Hessiano(a) de f(x). • Alguns autores denominam o determinante de ∇2f(x) de Hessiano de f(x), ou seja, H(x) = det[∇2f(x)] , porém esta não é a nomeclatura adotada neste trabalho. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 47 2.14. Interpretação da Função Objetivo em Termos de uma Aproximação Quadrática Se uma função de duas variáveis é ou pode ser aproximada por uma função quadrática do tipo ( ) 21122222211122110 xxbxbxbxbxbbxf +++++= Equação 2-II então os autovalores da matriz Hessiana de f(x) podem ser utilizados para interpretar a natureza de f(x) nos pontos estacionários x* . A tabela 2-b lista as conclusões que podem ser alcançadas examinando os autovalores de H(x*), e as ilustram as diferentes superfícies correspondentes a cada tipo de função quadrática. Tabela 2-B: Interpretação geométrica de uma função quadrática Caso Relação entre os autovalores Sinal α1 Sinal α2 Tipo de contorno Interpretação geométrica Ponto estacionário Problema de otimização Fig. 1 α1 = α2 - - Círculo Monte circular Máximo bem-comportado e raro na prática 2.11 2 α1 = α2 + + Círculo Vale circular Mínimo bem-comportado e raro na prática 2.11 3 α1 > α2 - - Elipse Monte elíptico Máximo bem-comportado e mais freqüente 2.12 4 α1 > α2 + + Elipse Vale elíptico Mínimo bem-comportado e mais freqüente 2.12 5 |α1| = |α2| + - Hipérbole Sela simétrica Ponto de sela degenerado 2.13 6 |α1| = |α2| - + Hipérbole Sela simétrica Ponto de sela degenerado 2.13 7 α1 > α2 + - Hipérbole Sela alongada Ponto de sela degenerado 2.13 8 α2 = 0 - Reta Condilheira estacionária Infinitos degenerado 2.15 9 α2 = 0 + Reta Vale estacionário Infinitos degenerado 2.15 10 α2 = 0 - Parábola Cordilheira crescente No ∞ degenerado 2.16 11 α2 = 0 + Parábola Cordilheira decrescente No ∞ degenerado 2.16 Uma função objetivo é dita bem-comportada quando seus contornos formam uma região convexa. Como mostrado na tabela 2-b os autovalores da matriz Hessiana de f(x) indicam a forma da superfície formada por f(x). Por sua vez os autovetores de H(x) correspondem às direções dos principais eixos dos contornos de f(x). Esta informação pode ser utilizada para definir direções de pesquisa mais eficientes. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 50 Figura 2-M: Geometria de FO de 2a ordem de duas VD’s: ponto de sela Figura 2-N: Geometria de FO de 2a ordem de duas VD’s: vale Figura 2-O: Geometria de FO de 2a ordem de duas VD’s: cordilheira decrescente Ricardo Kalid - Otimização de Processos – kalid@ufba.br 51 2.15. Exercícios E2.1. Utilizando a definição de inversa de uma matriz, expresse o vetor x em função do vetor z. z1 = x1 + x2 z2 = 2x1 + x2 E2.2. Utilizando a técnica da matriz aumentada (operadores-linha) calcule a inversa de A: ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 12 41 A Compare a matriz A-1 com a inversa calculada pelo comando inv(A) do MATLAB. E2.3. Calcule a solução dos seguintes sistemas de equações lineares: A x = b (a) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 0 1 ; 21 12 bA (b) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 5 6 ; 11 22 bA (c ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 3 6 ; 11 22 bA E2.4. Utilizando a equação det[A - αI] = 0 calcule os autovetores e autovalores da matriz ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 12 21 A . Compare os resultados com os obtidos pelo comando [v,d] = eig(A) E2.5. Determine os pontos estacionários da função f x x x x ( ) = + + + 1 2 14 3 2 e classifique-os. E2.6. Verifique a região na qual f(x) e f'(x) são contínuas: a) ( ) x xf 1= b) ( ) ( )xxf ln= E2.7. Estude a variação da função f(x) = 2x2 - x3 , isto é, estabeleça as regiões de concavidades e/ou convexidades, calcule os pontos estacionários e obtenha o valor de de f(x) nesses pontos. Trace gráficos de f(x). E2.8. Determine a convexidade/concavidade das funções abaixo: a) f(x) = 2x12 - 3x1x2 + 2x22 b) f(x) = x12 + x1x2 + 2x2 + 4 c) f(x) = 2x1 + 3x2 + 6 d) f(x) = 2x12 + 2x1x2 + 1,5x22 + 7x1 + 8x2 + 24 Ricardo Kalid - Otimização de Processos – kalid@ufba.br 52 E2.9. Determine a região viável e verifique, através da Hessiana e de gráficos, se a mesma é convexa ou não: a) -x12 + x2 > 1 ; x1 - x2 > -2 b) -x12 + x2 > 1 ; x1 - x2 < -2 c) x1 < 6 ; x2 < 6 ; x1 > 0 ; x2 > 0 ; x1 + x2 < 6 E2.10.Classifique os pontos estacionários e verifique se as condições necessárias e as suficientes são atendidas: a) f(x) = x4/3 b) f(x) = 4x3 c) f(x) = x13+ x22 - 3x1 + 8x2 + 2 d) f(x) = 4 + 4,5x1 - 4x2 + x12 + 2x22 - 2x1x2 + x14 - 2 x12x2 E2.11.Happel e Jordan (Chemical Process Economics, Marcel Dekker, New York, 1975, p. 178) desenvolveram uma função objetivo (custo) para uma coluna de destilação: f(n,P,R) = 14720(100 - P) + 6560R - 30,2PR + 6560 - 30,2P + 19,5n(5000R - 23PR +5000 - 23P)0,5 + 23,2(5000R - 23PR + 5000 - 23P)0,62 onde n é o número de estágios teóricos, R é a razão de refluxo e P é o percentual de recuperação da corrente do fundo. Qual o ponto ótimo? Neste ponto a FO é convexa? Existe alguma região não convexa nas vizinhanças do ponto ótimo? E2.12.Uma reação homogênea converte 2 compostos orgânicos (A e B) num produto P pelo aquecimento do meio reacional. Os reagentes podem ser injetados no reator, enquanto que vapor passa por uma serpentina para aquecer o meio reacional. O produto P pode ser vendido a 50$/kg-mol. Para 1 kg-mol de alimentação, o custo da alimentação (em $/kg-mol) é uma função da fração molar de A (xA) e é dado por f(xA) = 2 + 10xA + 20xA2 . O custo do vapor (em dolares) é uma função de S (kg de vapor/kg-mol de alimentação) e é g(S) = 1 + 0,003S + 2x10-6S2 . O rendimento a P é dado por yP(xA,S) = 0,1 + 0,3xA + 0,001S + 0,0001xAS , onde as unidades de yP são kg-mol de P/kg-mol da alimentação. Pede-se a) Obtenha a função lucratividade (base de 1 kg-mol de alimentação) em função de xA e S b) Maximize a FO sujeita às seguintes restrições: 0 < xA < 1 e S > 0 c) Demonstre matematicamente se f é uma função côncava ou convexa? d) A região viável é convexa? Por que? Ricardo Kalid - Otimização de Processos – kalid@ufba.br 55 extremo local diferente do global, mas que apresente características mais apropriadas para o sistema. Como foi dito no Capítulo 1, a expressão matemática de um problema de otimização tem duas partes: P1. a função objetivo (FO), ou seja, o critério de desempenho a ser atingido, cuja especificação é obrigatória; P2. as restrições ou funções de restrição (FR), que estão quase sempre presentes. De maneira geral as restrições podem ser escritas sob a forma de equações (algébricas ou diferenciais) e/ou inequações (algébricas ou diferenciais). Genericamente o problema de otimização pode ser formulado da seguinte forma ( ) ( ) 0 e/ou 0 a sujeito )(min ≤ = xh xg xf x Equação 3-A Para resolver os problemas linearres ou não-lineares de otimização com ou sem restrições existem inúmeros métodos numéricos. Nos capítulos 4, 5, e 6 descreveremos alguns dos algoritmos mais utilizados nos problemas de engenharia química. Nas seções 3.1, 3.2 e 3.3 estudaremos como definir a FO e suas FR's. 3.1. A Função Objetivo (FO) Devemos ser capazes de traduzir expressões verbais do tipo maximização do lucro ou minimização dos custos em termos matemáticos. Devemos expressar a FO em termos de unidades monetárias ou em unidades quantificáveis, e omitir expressões filosóficas do tipo "construir um mundo melhor" ou "desenvolver uma sociedade mais humana". Também não trataremos de problemas que envolvem múltiplas funções objetivos. Na definição da FO podemos considerar apenas objetivos econômicos (maximizar a lucratividade, por exemplo) ou apenas objetivos operacionais (diminuir a diferença entre o valor desejado e o valor medido na operação de um equipamento) ou combinar os dois tipos de objetivos numa única FO. Todos os procedimentos numéricos de otimização requerem que seja definido um critério de parada ou tolerância, pois a solução exata nunca é encontrada, mas apenas uma aproximação da mesma. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 56 3.1.1. Tolerância ou Critério de Parada Como dissemos anteriormente, o ponto ótimo obtido está numa vizinhança do ótimo verdadeiro. Esta vizinhança pode ser tão pequena quanto a precisão das informações utilizadas, mas sempre será necessário arbitrar tal aproximação, isto é, a tolerância. Um dos seguintes critérios podem ser utilizados para interromper os processos iterativos de procura do ponto ótimo: Erro absoluto na FO: ( ) ( ) 11 ε<−+ kk xfxf Equação 3-B Erro relativo na FO: ( ) ( ) ( ) 2 1 ε<− + k kk xf xfxf Equação 3-C Erro absoluto no ponto ótimo: nixx ki k i ≤∀<− + 3 1 ε Equação 3-D ou ( ) 4 1 211 ε<−=− ∑ = ++ n i k i k i k i k i xxxx Equação 3-E Erro relativo no ponto ótimo: ni x xx k i k i k i ≤∀< −+ 5 1 ε Equação 3-F Erro absoluto na direção de busca: ( ) ( ) ni x xfxf n i x k k ≤∀< ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ =∇ ∑ = 6 2 1 ε ∂ ∂ Equação 3-G Na verdade devemos utilizar no mínimo 2 critérios para caracterizar que as tolerâncias especificadas foram alcançadas. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 57 3.1.2. Objetivos Econômicos A FO pode ser escrita em relação ao lucro e/ou custo, ou simplesmente em termos da lucratividade de um dado projeto. O termo lucratividade é utilizado para mensurar a quantidade de lucro gerado, tal medida pode ser expressa de diferentes maneiras. Qual a forma mais adequada depende do ambiente onde esta inserido o problema. Para desenvolver uma FO que leve em conta a lucratividade de um dado empreendimento é necessário compreender conceitos tais como depreciação, fluxo de caixa, custos de capital, custos operacionais, valor presente e valor futuro, inflação, vida econômica, investimento inicial, capital de giro, taxa de retorno de investimento, valor presente líquido, tempo de retorno, etc. Infelizmente não temos tempo para discutir com profundidade tais temas e por isso não os trataremos, com excessão para os custos de capital e operacional que abordaremos a seguir. Podemos estabeler os seguintes objetivos econômicos: OE1. A minimização dos custos ou maximização dos lucros operacionais. OE2. A minimização dos custos de investimento (ou de capital). OE3. A minimização de uma função que inclua os custos de capital e operacional. A primeira categoria de funções objetivo (OE1) envolvem apenas custos variáveis e a receita com as vendas e é encontrada em situações onde o custo de capital é uma quantia fixa (os equipamentos já estão comprados e funcionando). A segunda categoria (OE2) aparece em situações onde não existem custos variáveis, por exemplo, em muitos problemas de projeto mecânico de equipamentos. A terceira categoria (OE3) inclue simultaneamente custo de capital e custo operacional. Por exemplo na decisão de implantar uma nova fábrica ou ampliar uma já existente. Este tipo de problema é mais complexo e difícil de equacionar. O Exemplo 3.1 apresenta um problema de otimização cuja FO é função apenas dos custos operacionais. Exemplo 3.1: FO função apenas dos custos e receitas operacionais. Tempo ótimo de campanha de um reator catalítico. Considere um ciclo de operação de um reator catalítico que sofre regeneração periódica de seu leito. Seja x1 o número de dias nos quais o catalisador é utilizado no reator e x2 o número de dias necessários para a regeneração. Portanto cada ciclo dura x1 + x2 dias. As seguintes hipóteses podem ser consideradas: H.1. A operação do reator apenas pode ser inicializada no turno da manhã, portanto x1 + x2 deve ser um número inteiro. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 60 custo = S x (peso do material empregado na fabricação do vaso) e a expressão matemática da FO é: ( ) $ em custo de unidade em 4 2,,,,onde min 2 1 1,,,, tDLDSLDtSf f LDtS ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ += ππρρ ρ Como S, ρ e t são constantes então a FO também pode ser ( ) peso de unidades em 2 ,,,onde min 2 2 2,, tDLDLDtf f LDt ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ += ππρρ ou ( ) area de unidades em 2 ,onde min 2 3 3, ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ += DLDLDf f LD ππ Note que as três FO's diferem apenas de fatores multiplicativos constantes, portanto a inclusão ou não desses fatores apenas afetará o valor da FO, mas não alterará o valor das variáveis de projeto. Logo, por simplicidade, utilizaremos a FO em unidades de área. Como o volume do vaso cilíndrico é constante existe uma relação entre altura e diâmetro do mesmo: VD LLDV 2 2 4 4 π π =⇔= Assim existe apenas uma variável de decisão e podemos reescrever a FO da seguinte forma ( ) area de unidades em 4 2 onde min 2 4 4 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ += D VDDf f D π Identificação da(s) restrições: não existem restrições. Para encontrar o ponto ótimo temos que diferenciar f4 em relação a variável de projeto igualar a zero, resultando em Ricardo Kalid - Otimização de Processos – kalid@ufba.br 61 1 D L sejaou 4 ementeconsequent e, 4 3/1 3/1 =⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛= ⎟ ⎠ ⎞ ⎜ ⎝ ⎛= opt opt opt VL VD π π Logo a razão (L/D)opt obtida é significativamente diferente da esperada (L/D = 3). Esta diferença pode ser devida a escolha de hipóteses inapropriadas. Brummerstedt (1944)# e Happel e Jordan (1975)@ desenvolveram hipóteses mais realistas para o problema do dimensionamento de vasos de pressão: H1. Topo e fundo tem o formato de elipses 2:1, com uma área dada por 2(1,16D2) = 2,32D2 H2. Custo de fabricação das extremidades é maior que para a lateral; Happel sugere um fator de 1,5. H3. O custo de fabricação por unidade de peso S ($/unidade de peso) e massa específica (ρ) são constantes. H4. A espessura t é função do diâmetro do vaso, do esforço a que o aço será submetido, da pressão e da corrosão admissíveis. Por exemplo, para uma pressão de projeto de 250 psi (17 atm) e uma corrosão prevista de 1/8 in, a espessura t (em polegadas) em função do diâmetro D (em pés): t = 0,0108D + 0,125 Identificação da(s) variáveis de decisão: considerando as discussões anteriores a única VD é o diâmetro D do vaso. # Brummerstedt, E. F., Natl. Pet. News, 36, R282, R362, R497 (1944). @ Happel, J. e D. G. Jordan, Chemical Process Economics, 2a. edição, Marcel Dekker, New York (1975). Ricardo Kalid - Otimização de Processos – kalid@ufba.br 62 Identificação do objetivo: é o mesmo anterior, ou seja, minimizar o custo de fabricação do vaso de pressão As novas hipóteses forçam a FO a ser expressa em dólares, pois a área ou o peso não são mais diretamente proporcionais ao custo, assim a nova FO é ( ) ( )[ ] $ em custo de unidade 32,25,1,,,,onde min 2 5 5,,,, tDLDSLDtSf f LDtS πρρ ρ += Substituindo t(D) em f5 e lembrando que S e ρ são constantes e que ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ += 34 2 DLDV π obtemos ( ) $ em custo de unidade 0263,03041,05,00432,0 min 32 6 6 DD D VVDfonde f D +++= Identificação da(s) restrições: não existem restrições. Ao solucionar este problema para diferentes níveis de pressão Happel apresentou a seguintes soluções: Ponto ótimo: (L/D)opt Pressão de projeto (psi) Capacidade (ft3) 100 250 400 2500 1,7 2,4 2,9 25000 2,2 2,9 4,3 Observe que o ponto ótimo se aproxima da razão 1 à medida que a pressão e a capacidade do tanque diminuem e que o valor empírico (L/D = 3) pode estar errado em ±50 % em relação ao (L/D)opt. Ainda é necessária uma interpretação cautelosa desses resultados pois não foram consideradas as perdas de material durante a fabricação. Consideremos agora um exemplo que leva em conta os custos de capital e operacionais. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 65 Exemplo 3.4: FO função apenas de um objetivo operacional/técnico. Reconciliação de dados de Processos. Suponha que as vazões de alimentação e descarga de um processo são medidas periodicamente. Determine o melhor valor da vazão mássica da corrente A (MA) em kg/h para o processo mostrado na figura 3-p. As correntes B e C foram medidas em intervalos de 1 hora. Figura 3-P: Diagrama de blocos do processo do Exemplo 3.4. Solução: Identificação da(s) variáveis de decisão: obviamente MA Identificação do objetivo: A FO pode ser definida como: minimizar a soma dos quadrados dos desvios entre MA medido e esperado Para matematizar a expressão da FO temos que inicialmente estabelecer algumas hipóteses simplificadores do sistema: H.1. O processo está em estado estacionário. H.2. As medições das correntes B (MB) e C (MC) estão corretas. Modelo do processo (balanço de massa): o que entra = o que sai MA + MC = MB (a) então a FO pode ser escrita como: ( ) ( ) ( )[ ]∑ = −−= exp. 1 2 min no i CB rec AA AM ii A MMMMfonde Mf (b) para os dados deste problema: f(MA) = [MA - (92,4 - 11,1)]2 + [MA - (94,3 - 10,8)]2 + [MA - (93,8 - 11,4)]2 (c) Identificação da(s) restrições: não existem restrições. MB1 = 92,4 kg/h MB2 = 94,3 kg/h MB3 = 93,8 kg/h A B C MC1= 11,1 kg/h MC2= 10,8 kg/h MC3= 11,4 kg/h Planta Ricardo Kalid - Otimização de Processos – kalid@ufba.br 66 Diferenciando (c) em relação a MA e igualando a zero, obtemos MAopt = 82,4. Neste caso MAopt = Mbmed - MCmed O valor mínimo que f(MA) pode alcançar é zero. Outros métodos de reconciliação de balanços de massa e energia são discutidos por Tamhane e Mah (1985)@. No Exemplo 3.5. temos outro caso de FO com critério apenas técnico/operacional. Exemplo 3.5: FO função apenas de um objetivo operacional. Tempo de residência ótimo de um reator a batelada. O esquema reacional abaixo acontece num reator a batelada bem agitado. Para as concentrações iniciais e constantes das taxas CAo = 50 g-mol/L k1 = 2,0 h-1 CBo = 5 g-mol/L k1 = 1,0 h-1 CCo = 0 g-mol/L k1 = 0,2 h-1 CDo = 0 g-mol/L k1 = 0,6 h-1 calcule o tempo de residência ótimo para obter um máximo rendimento a CB. Solução: Identificação das variáveis de decisão: a priori k1, k2, k3, k4, CA, CB, CC, CD e tempo de residência Identificação do objetivo: A FO pode ser definida como maximizar a concentração do produto desejado B ou ( )max C k k k k C C C C tB A B C D1 2 3 4, , , , , , , , (a) Para matematizar a expressão da FO temos que inicialmente estabelecer algumas hipóteses simplificadores do sistema: H1. O sistema não troca massa com o meio. H2. Reator opera isotermicamente. H3. Equações das taxas de 1a ordem. @ Tamhane, A. C., e R. S. H. Mah, "Data Reconciliation and Gross Error Detection in a Chemical Process Network", Technometrics, 27: 409 (1985). B A C D k1 k2 k3 k4 Ricardo Kalid - Otimização de Processos – kalid@ufba.br 67 H4. Reação em fase líquida, portanto podemos admitir volume constante. Modelo do processo balanço de massa por componente: acúmulo = formado - consumido AB A CkCk dt dC 12 −= (b) BA B CkkkCk dt dC )( 4321 ++−= (c) B C Ck dt dC 3= (d) B D Ck dt dC 4= (e) balanço de massa global: massa inicial = massa final ou, lembrando que os coeficientes estequiométricos das reações são todos iguais, em termos do número de moles, o balanço global pode ser escrito como: no moles inicial = no moles final Lembrando ainda que o volume reacional é constante, o balanço global pode ser escrito em termos da concentração dos reagentes e produtos como: concentração inicial reag. + prod. = concentração final reag. + prod. ou CAo + CBo + CCo + CDo = CA + CB + CC + CD (f) Queremos apenas encontrar o tempo de residência que maximiza CB, portanto tentaremos escrever a CB como uma função apenas do tempo. Resolvendo a equação (c ) para CA ( ) ⎥⎦ ⎤ ⎢⎣ ⎡ +++= B B A Ckkkdt dC k C 432 1 1 (g) Diferenciando (g) em relação ao tempo ( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ +++= dt dC kkk dt Cd kdt dC BBA 4322 2 1 1 (h) Substituindo (g) e (h) em (b) obtemos uma equação de segunda ordem linear homogênea em CB: ( ) ( )[ ] 021432143212 2 =−+++++++ B BB Ckkkkkk dt dC kkkk dt Cd (i) cuja solução, para este caso em que as raizes da equação característica são reais e distintas, é dada por Ricardo Kalid - Otimização de Processos – kalid@ufba.br 70 R2. física, por exemplo limites máximo para temperaturas a depender dos materiais presentes no processo, máxima capacidade produtiva da unidade, fechamento dos balanços de massa e energia, fração molar deve estar entre 0 e 1, temperatura e pressão absoluta maior que 0, etc; R3. mercadológica, por exemplo deve atender a uma demanda mínima e/ou consumir uma quantidade máxima de matéria-prima e/ou de energia, etc. Portanto para que o problema de otimização tenha uma solução correta tais restrições devem ser explicita e matematicamente escritas. A formulação do problema de otimização só estará completa se estabelecemos a FO e as suas FR's. No Exemplo 3.1. a hipótese H1. na verdade é uma restrição operacional. Não a levamos em conta no momento da solução do problema matemático, mas para validar a resposta tivemos que arredondar a mesma. Não raro este procedimento é observado na prática, porém devemos ter sempre em mente que encontraremos apenas uma solução sub-ótima. No Exemplo 3.2. a consideração de volume fixo é uma restrição física (capacidade necessária para o equipamento). Neste caso foi possível incorporar diretamente a restrição na FO, o que simplifica a resolução do problema. Também nos problemas de controle ótimo devemos incorporar restrições quanto a variação das MV's e aos valores extremos que as mesmas podem atingir. Assim a equação 3-l deve ser reescrita da maneira que se segue ( )[ ] ( )[ ] ( )[ ] ( )[ ] ( )( )[ ]{ } ( ) ( ) ( ) ( ) ( ) ( )maxmin maxmin 2 2 T 1 T a sujeito min tMVtMVtMV tMVtMVtMV tMVtMVWtMVtEWtE Δ≤Δ≤Δ ≤≤ +ΔΔ+ θ Equação 3-M 3.3. Otimização On-Line Neste ponto vale a pena comentar um pouco sobre a otimização on-line de processos químicos. Como o próprio nome já diz, é um problema de otimização, como os vistos anteriormente, cujas variáveis são atualizadas automativamente com as atuais condições operacionais do processo e resolvido, isto é, o valor das variáveis de decisão é encontrado, com uma freqüência relativamente elevada, por exemplo 20 vezes por hora. Porém só podemos pensar na implementação da otimização on-line do processo, ou seja, na implantação de uma metodologia que periodicamente defina qual os set-points para que uma determinada unidade ou planta opere nas condições economicamente ótimas, uma vez estabilizada a operação da Ricardo Kalid - Otimização de Processos – kalid@ufba.br 71 unidade, por exemplo através da implementação de um algoritmo tipo MPC♦ (multivariable predictive control ou controle preditivo multivariável, tais como DMC, DMCL, LDMC, QDMC, GPC). Essa otimização em linha pode ser implementada de duas formas: a) Um MPC seguindo setpoints enviados por um otimizador hierarquicamente superior, figura 3-q, neste caso é necessário um sistema supervisor que calcule e envie para o MPC os SP’s. O T IM IZA D O R P LA N T A M V’sP V’S M P C SP ’s das P V’sM V’s e P V’s Figura 3-Q: MPC Supervisionado por Otimizador On-Line b) MPC com otimização intrínseca, figura 3-r, neste caso a condição operacional ótima é obtida simultaneamente com o atendimento das restrições operacionais. Esta configuração é mais moderna e aparenta ser mais eficiente que a anterior, mas ainda se encontra em fase de desenvolvimento nos meios acadêmicos. MPC COM OTIMIZAÇÃO PLANTA MV’sPV’S Figura 3-R: MPC com Otimização Intrínseca ♦ Maiores informações sobre controle preditivo podem ser encontradas em DMC: Controle por Matriz Dinâmica. Apóstila de Curso Ministrado na Copene de autoria de Ricardo de Araújo Kalid e Darci Odloak, setembro de 1994. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 72 3.4. Exercícios E3.1. Uma central de insumos básicos produz, armazena e distribui várias matérias-primas para outras empresas, conforme a figura 3-s. Processo 3 Processo 1 x1,1 x3,1 x2,1 x4,1 x1,2 x3,2 x2,2 x4,2 x1,3 x3,3 x2,3 x1,4 x3,4 x2,4 T1 T4 T3 T2 T6 T5 T7 x1 x4,4 _ _ _ x5,2 Figura 3-S: Diagrama de blocos do parque de tancagem onde Ti (i = 1, ... , 7) são tanques de armazenamento. Para a distribuição dos produtos são feitos contratos por períodos pré-determinados, com demandas mínimas a serem necessariamente atendidas pela empresa. Desenvolva uma Função Objetivo com suas restrições que minimizem o custo de armazenamento. Escreva a FO em função dos custos Ci e da massa armazenada em cada tanque mi . onde: Ci = custo do produto armazenado; mi = massa do produto armazenado. Dados referentes ao processo: Restrições de processamento : 10 ton/h < x1 < 50 ton/h 6 ton/h < x3 < 30 ton/h. Capacidade de armazenamento dos tanques: 100 ton < mi < 500 ton, i = 1, ..., 7 Inicialmente todos os tanques estão com 120 ton. O balanço de massa tem que ser atendido: 11,41,41,4 33,33,33,311,31,31,3 33,23,23,211,21,21,2 31,33,13,111,11,11,1 xbax xbaxxbax xbaxxbax xbaxxbax += +=+= +=+= +=+= As demandas dos seguintes produtos são conhecidas: Ricardo Kalid - Otimização de Processos – kalid@ufba.br 75 Figura 3-U: Diagrama de blocos de uma refinaria Tabela 3-A: Dados de rendimento da refinaria e do mercado Porcentagem de rendimento volumétrico Produção máxima Petróleo tipo 1 Petróleo tipo 2 (bbl/dia) Gasolina 80,0 44,0 24000 Querosene 5,0 10,0 2000 Óleo combustível 10,0 36,0 6000 Resíduo 5,0 10,0 Custo de processamento ($/bbl) 0,5 1,0 Pede-se: a) Desenvolva uma função objetivo linear. b) Estabeleça as restrições. c) Desenhe num gráfico as curvas das restrições. d) Identifique no gráfico a região viável. e) Calcule nas interseções entre as curvas das restrições o valor alcançado pela FO e verifique qual o ponto ótimo. Ref: Edgar, T. F. & Himmelblau, D. M. "Optimization of Chemical Processes". E3.4. Programação de uma unidade de craqueamento térmico. A figura 3-v mostra os vários componentes da alimentação e a corresponde conversão aos vários produtos devido ao craqueamento térmico. As possíveis cargas incluem etano, propano, gasolina natural debutanizada (DNG) e gasóleo, algumas das quais podem ser alimentadas simultaneamente. Oito produtos podem ser produzidos em várias proporções de acordo com a matriz mostrada na tabela 3-b. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 76 Figura 3-V: Diagrama de blocos do craquemento térmico Tabela 3-B: Rendimento aos produtos (fração mássica) Alimentação Produto Etano Propano Gasóleo DNG Metano 0,07 0,25 0,10 0,15 Etano 0,40 0,06 0,04 0,05 Etileno 0,50 0,35 0,20 0,25 Propano --- 0,10 0,01 0,01 Propileno 0,01 0,15 0,15 0,18 Butadieno 0,01 0,02 0,04 0,05 Gasolina 0,01 0,07 0,25 0,30 Óleo combustível --- --- 0,21 0,01 A capacidade de processamento do craqueador é de 90,7 ton/h (vazão total baseada numa mistura média). Etano usa o equivalente a 1,1 ton da capacidade por ton de etano; propano usa o equivalente a 0,9 ton da capacidade por ton de propano; gasóleo usa o equivalente a 1,1 ton da capacidade por ton de gasóleo; e DNG utiliza uma razão de 1,0. As unidades a jusante tem uma capacidade máxima de processamento limitada a 22,7 ton/h de etileno e 9,1 ton/h de propileno. O combustível necessário ao processamento de cada tipo de substância presente na carga é apresentado na tabela 3-c. Tabela 3-C: Consumo de combustível por tipo de carga Componente da carga Combustível requerido para processamento Etano 4,65x103 kcal/kg Propano 2,79 x103 kcal/kg Gasóleo 2,17 x103 kcal/kg DNG 2,53 x103 kcal/kg Craqueamento térmico Óleo combustível Metano Óleo combustível DNG Gasóle Propan o Etano Etano Propan Etileno Propilen o Butadieno Gasolina Ricardo Kalid - Otimização de Processos – kalid@ufba.br 77 O metano e óleo combustível produzido pela unidade são reciclados como combustível. Todo o etano e propano não convertidos são reciclados para a alimentação. A energia fornecida por cada tipo de combustível é mostrada na tabela 3-d. Tabela 3-D: Calor fornecido pelos combustíveis Tipo de combustível Energia liberada na combustão Gás natural 11,96 x103 kcal/kg Metano 11,96 x103 kcal/kg Óleo combustível 10,00x103 kcal/kg Devido a perda de calor a energia necessária para a pirólise, uma quantidade fixa de energia igual a 5,04 kcal/h, deve ser fornecida ao forno. O preço de compra dos combustíveis, da carga e dos produtos são mostrados na tabela 3-e. Tabela 3-E: Preço de compra e venda dos insumos Subtância Preço ($/kg) Etano (preço de compra) 0,1444 Propano (preço de compra) 0,2145 Gasóleo (preço de compra) 0,2756 DNG (preço de compra) 0,2235 Metano (preço de venda como combustível) 0,1186 Etileno (preço de venda) 0,3913 Propileno (preço de venda) 0,3040 Butadieno (preço de venda) 0,5873 Gasolina (preço de venda) 0,2189 Óleo combustível (preço de venda como combustível) 0,0992 Pede-se: a) Desenvolva uma função objetivo linear que maximize o lucro da unidade. b) Estabeleça as restrições operacionais e/ou físicas e/ou de merdado a que o craqueador esta submetido. Ref: Edgar, T. F. & Himmelblau, D. M. "Optimization of Chemical Processes". E3.5. Discutas as duas FO's abaixo, compare-as apontando suas igualdades e diferenças. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 80 ( ) ( ) ( ) ( ) ( ) ( ) 0 0 0 2 1 2 1 == == == xf x xf xf x xf xf x xf nx n x x ∂ ∂ ∂ ∂ ∂ ∂ Equação 4-A Cada derivada parcial igualada a zero constitui uma equação não-linear, assim o problema de OMSR transformou-se num sistema algébrico de equações não- lineares, um problema tão difícil de resolver quanto o original. Além disso, enquanto que a matriz Hessiana H(x) de f(x) é simétrica, o Jacobiano das equações das condições necessárias tem uma estrutura arbitrária, o que dificulta a resolução do problema. Assim, muitos engenheiros preferem atacar o problema de minimização/maximização por MD's em lugar de utilizar MI's. A seguir, na seção 4.1 apresentaremos os MI's mais utilizados em problemas de OUSR da engenharia química, enquanto na seção 4.2 estudaremos os MD's. 4.1. Métodos Indiretos (MI) para OUSR Existem basicamente 4 procedimentos para resolver problemas de OUSR envolvendo a aplicação das condições necessárias de otimalidade para uma função: MI1. Método de Newton MI2. Aproximação por diferenças finitas do método de Newton, também chamado de método de quasi-Newton MI3. Método da secante Para comparar a eficiência desses métodos é conveniente examinar a taxa de convergência de cada um, isto é, a velocidade com que a solução é alcançada. A taxa de convergência pode ser expressa de várias maneiras, mas uma classificação comum é a seguinte: Linear: 10 1 ≤≤≤ − − ∗ ∗+ αα xx xx k k Equação 4-B Onde o super-escrito k se refere a iteração k e não é uma potência de x. Na resolução de problemas práticos se o método apresenta uma taxa de convergência linear sua velocidade é lenta. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 81 O símbolo ||.|| indica a norma Euclidiana de um vetor. Para x n- dimensional: ∑ = = n i ixx 1 2 Equação 4-C Ordem p: 10 1 ≤≤≤ − − ∗ ∗+ ααpk k xx xx Equação 4-D Na resolução de problemas práticos se o método apresenta uma taxa de convergência de ordem p > 1 sua velocidade é rápida. Se p = 2 , a ordem de convergência é dita quadrática. Superlinear: ∞→→<= − − ∗ ∗+ ∞→ kquandoeou xx xx kkk k k 00lim 1 αα Equação 4-E Na resolução de problemas práticos se o método apresenta uma taxa de convergência superlinear sua velocidade é rápida. O algoritmo de qualquer um dos MI's para minizar f(x) pode ser resumido da seguinte forma: P1. Escolha um intervalo de busca que contenha o ponto de mínimo e que a FO seja unimodal neste intervalo. P2. Defina uma tolerância. P3. Aplique o método até que a tolerância seja alcançada, certificando-se a cada iteração que f(xk+1) < f(xk) . Quanto mais positivo for f"(xk) , ou sua aproximação, mais rapidamente f(x) diminuirá. Para maximizar f(x) , minimize -f(x) . 4.1.1. Método de Newton Lembremos sempre que o objetivo é minimizar (ou maximizar) uma certa função f(x), lembremos também que no momento estamos estudando apenas o caso monovariável. Vamos rever o método de Newton para resolução de uma equação não- linear. Este método está baseado na expansão da função em série de Taylor e truncamento da mesma no segundo termo. Assim expandindo f(x) em torno de um ponto xk obtemos: Ricardo Kalid - Otimização de Processos – kalid@ufba.br 82 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )kkk x kk xfxxxfxf dx xdfxxxfxf k ′−+=⇔−+= . Equação 4-F como qualquer função pode ser escrita da forma f(x) = 0 , então podemos reescrever a equação 4-f, se a aproximação e truncamento executado acima forem válidos, da seguinte maneira: ( ) ( ) ( ) ( ) ( )( )k k kkkk xf xfxxxfxxxfxf ′ −=⇔=′−+= 0. Equação 4-G Ao trocar x por xk+1 na equação 4-g obtemos o método de Newton para resolução de uma equação não-linear. Agora podemos resolver o problema de OUSR, ou seja resolver a equação não-linear f'(x) = g(x) = 0 pelo método de Newton, utilizando a fórmula de recorrência: ( ) ( ) ( ) ( )k k kk k k kk xf xfxx xg xgxx ′′ ′ −=⇔ ′ −= ++ 11 Equação 4-H Como a convergência deste método não é garantida temos que nos assegurar que a cada iteração f(xk+1) < f(xk) para que o mínimo seja alcançado [f(xk+1) > f(xk) para o máximo]. As vantagens do método de Newton são: V1. O procedimento é tem convergência quadrática local (p = 2 na equação 4- d) se f "(x) ≠ 0 no ponto ótimo. V2. Para uma FO quadrática o mínimo é obtido em uma iteração, pois a expansão e truncamento da série de Taylor da função f(x) até o termo da 2a derivada (inclusive) é exata. V3. Converge rapidamente quando a estimativa inicial é boa. As desvantagens do método de Newton são: D1. Só se aplica a funções onde existam f '(x) e f "(x) . D2. Deve-se calcular a cada iteração f '(x) e f "(x) . D3. Sensível à estimativa inicial. D4. Se f "(x) → 0 a convergência é lenta. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 85 para um problema monovariável, apesar dessa convergência ser menor que no método quasi-Newton, termina por ser mais eficiente pois precisa avaliar um número menor de funções. 4.2. Métodos Diretos (MD) para OUSR Podemos subdividir os MD's em dois grandes grupos: G1. Métodos por diminuição da região de busca. G2. Métodos por aproximação polinomial da FO. 4.2.1. Métodos por Diminuição da Região de Busca O algoritmo desses métodos para minimização de uma f(x) pode ser assim descrito: P1. Escolha uma região viável (a,b) que contenha uma função unimodal (função convexa). P2. Calcule f(x) em dois ou mais pontos. P3. Compare os valores de f(x) eliminando as regiões que apresentam valores de f(x) maiores. Na figura 4-a(a), a região a esquerda de x1 contém os valores de f(x) > f(x1) e por isso deve ser eliminada na busca do mínimo de f(x) . Analogamente a região à direita de x2 na figura 4-a(c), f(x) > f(x2) , também deve ser descartada. Porém na figura 4-a(b) a função é suave e o ponto de mínimo pode estar em qualquer lugar do intervalo (a,b). Na verdade quando estamos bem próximos do ponto ótimo este caso ocorre e b - a ≈ 0 . Figura 4-A: 3 possibilidades para f(x) após dividido o intervalo (a,b) por 2 pontos. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 86 Conforme seja escolhida a localização dos pontos x1 e x2 este método torna-se mais ou menos eficiente: M1. Três intervalos iguais de busca, isto é, x1 - a = b - x2 = x2 - x1 neste caso o intervalo é reduzido de 1/3 a cada iteração, assim se L0 é o intervalo de busca original (b - a) e Lk é o intervalo após k interações, então 0 3 2 LL k k ⎟ ⎠ ⎞ ⎜ ⎝ ⎛= Equação 4-Q Observe que é necessário avaliar duas vezes a função f(x) a cada iteração. M2. Método da seção áurea, neste caso o intervalo eliminado manterá sempre uma mesma proporção com o intervalo original. Na figura 4-b e na equação 4-r observamos as distâncias utilizadas neste procedimento. Figura 4-B: Seção áurea 618,0 1 1 1 =⇒ − ==⇒= + =+ Δ Δ y y y x y yx y y yx yx Equação 4-R e ( ) 0618,0 LL kk = Equação 4-S Seja ak e bk são os contornos do intervalo de busca na iteração k , então o algoritmo para este método é o seguinte: P1. determine ak e bk , e a tolerância ε P2. Calcule Lk = bk - ak , x1k = ak + 0,382.Lk e x2k = bk - 0,618.Lk P3. Calcule f(x1k) e f(x2k) P4. Se f(x1k) > f(x2k) então ak+1 = x1k , bk+1 =bk , x1k+1 = x2k Se f(x1k) < f(x2k) então bk+1 = x2k , ak+1 =ak , x2k+1 = x1k P5. Se a tolerância não foi atingida volte ao passo P2. Veja a figura 4-c. a bx y Ricardo Kalid - Otimização de Processos – kalid@ufba.br 87 Figura 4-C: Seleção dos pontos interiores pela seção áurea 4.2.2. Métodos por Aproximação Polinomial - Interpolação Quadrática Outros tipos de MD são os métodos por aproximação polinomial, entre os quais destacamos o da interpolação quadrática. Nesses métodos a função f(x) é aproximada por extrapolação ou interpolação utilizando um polinômio. Coggins (1964)♣ discutiu várias técnicas envolvendo um polinômio aproximador, chegando à conclusão que esses são levemente mais eficientes que o método da seção áurea. Admitindo que f(x) é unimodal e conhecendo um intervalo de busca que contenha o ponto de mínimo, avaliamos f(x) em x1, x2 e x3 e obtemos f(x1) , f(x2) e f(x3) . Suponha que a função f(x) possa ser aproximada, numa vizinhança do ponto ótimo xf(x)* , por uma função quadrática g(x) então ( ) ( ) ( ) ( ) ( ) ( ) 23333 2 2222 2 1111 .. .. .. xcxbaxgxf xcxbaxgxf xcxbaxgxf ++=≅ ++=≅ ++=≅ Equação 4-T Mas o mínimo da função quadrática g(x) é dado por é dado por ( ) c bx xg 2 −= ∗ Equação 4-U Resolvendo a equação 4-t para obter os valores de b e c e substituindo o resultado na equação 4-u obtemos ♣ Coggins, G. F., Univariate Search Methods, Imperial Chemical Industries Ltd., Central Instr. Lab. Res. Note 64/11, 1964 Ricardo Kalid - Otimização de Processos – kalid@ufba.br 90 o julgamento do engenheiro é essencial na aceitação ou não das soluções encontradas pelos pacotes computacionais. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 91 4.4. Exercícios E4.1. Utilizando variados métodos indiretos obtenha o extremo das funções a) y1 da Equação 1-A b) f(x) do exercício E2.5 c) f(x) do exercício E2.6.a d) f(x) do exercício E2.6.b e) f(x) do exercício E2.7 f) f(x) do exercício E2.10.a g) f(x) do exercício E2.10.b E4.2. No Exemplo 3.2, verifique se os pontos de mínimo da equação f6(D) para uma pressão de projeto de 250 psi, estão corretos. Resolva este problema utilizando um método indireto e um método direto. Compare o número de iterações necessário para cada procedimento. E4.3. No Exemplo 3.5, verifique se o ponto de máximo da equação (j) está correto. Resolva este problema utilizando o método de Newton, o quasi- Newton e o da secante. Compare o número de iterações necessárias para cada procedimento. E4.4. Utilizando a função fmin('f',x) do toolbox de otimização do MATLAB resolva os problemas E4.2 e E4.3. Modifique os parâmetros do método numérico e verifique o comportamento do mesmo. Se necessário digite help fmin para ver como é a sintaxe deste comando, ou recorra ao manual de referência do toolbox de otimização do MATLAB. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 92 5. Otimização Multidimensional Sem Restrições (OMSR) A otimização numérica de FO multivariáveis não-lineares requer técnicas eficientes e robustas. A eficiência é importante pois os procedimentos para solucionar o problema são sempre iterativos e a tentativa-e-erro torna-se impraticável quando estão envolvidas mais de 3 ou 4 variáveis de decisão. A robustez (habilidade de encontrar a solução) é uma propriedade desejável pois as FO's não-lineares tem um comportamento imprevisível, pois podem ter múltiplos extremos, pontos de sela, regiões convexas e não convexas, etc. Mas felizmente estão disponíveis algoritmos numéricos poderosos e que solucionam a maioria dos problemas práticos. Neste capítulo discutiremos a solução dos problemas OMSR que tem a seguinte formulação geral encontre x* = [x1 x2 ... xn]T que minimiza f(x) = f(x1, x2, ..., xn) A maioria dos procedimentos interativos tem duas etapas que se alternam até que a tolerância desejada seja alcançada, assim podemos descrever os algoritmos da seguinte forma P1. escolha uma estimativa inicial para o vetor das variáveis de decisão x0 = [x1 x2 ... xn]T , P2. calcule a direção de busca sk , P3. minimize, seguindo a direção sk, a FO, isto é calcule xk+1 = xk + Δxk , onde Δxk = λsk, P4. verifique se a tolerância foi atingida, se não volte ao passo P2. Os vários métodos de programação não-linear (NLP - nonlinear programming) diferem uns dos outros principalmente na maneira como obtem a direção de busca sk . Alguns NLP utilizam as derivadas para definir a direção de busca (métodos indiretos) ou apenas a FO (métodos diretos). Os métodos que utilizam derivadas analíticas em geral usam menos tempo de computação, porém requerem maior tempo do engenheiro para a análise do problema. Também é importante lembramos que a natureza da FO influencia na escolha e no desempenho do algoritmo de otimização. A seguir, na seção 4.1 apresentaremos os MI's mais utilizados em problemas de OMSR da engenharia química, enquanto na seção 4.2 estudaremos os MD's. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 95 D1. Se os eixos das curvas de nível da FO são paralelos aos eixos formados pelas variáveis de decisão, por exemplo quando não existe interação entre as variáveis de decisão [ f(x1,x2) = x12 + x22 ], a direção de busca apontará diretamente para o ótimo, figura 5-be figura 5-c. Porém quando esses eixos não são paralelos, a taxa de convergência cai drasticamente, por exemplo quando existe interação entre as variáveis de busca [ f(x1,x2) = x12 + x22 + x1.x2 ] , figura 5-d e figura 5-e. Podemos resumir esta dificuldade dizendo que este método é muito sensível a forma das curvas de nível. D2. O gradiente descendente determina um ponto estacionário, não distinguindo se o mesmo é um extremo ou ponto de sela. Se for um ponto de sela é necessário empregar um método que não usa o gradiente para mover o ponto para fora dessa região. A estacionaridade do ponto pode ser testada examinando a matriz Hessiana da FO, conforme descrito no capítulo 2. Ricardo Kalid - Otimização de Processos – kalid@ufba.br 96 D3. x1 = -2:0.1:2; x2 = -2:0.1:2; [X1,X2] = meshgrid(x1,x2) ; fx1 = X1.^2 + X2.^2; mesh(x1,x2,fx1) -2 -1 0 1 2 -2 -1 0 1 2 0 2 4 6 8 Figura 5-B: Gráfico de fx1 v1 = [0.1,0.5,1,1.5,2,3,4,6]; c = contour(x1,x2,fx1,v1); clabel(c); -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0.1 0.5 1 1.5 2 3 4 6 6 6 6 Figura 5-C: Curvas de nível de fx1 x1 = -2:0.1:2; x2 = -2:0.1:2; [X1,X2] = meshgrid(x1,x2) ; fx2 = X1.^2 + X2.^2 + X1.*X2; mesh(x1,x2,fx2) -2 -1 0 1 2 -2 -1 0 1 2 5 10 15 20 Figura 5-D: Gráfico de fx2 v2 = [0.1,0.5,1,1.5,2,3,4,6]; c = contour(x1,x2,fx,v2); clabel(c); -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0.1 0.5 1 1.5 2 3 4 4 4 4 6 6 Figura 5-E: Curvas de nível de fx2 5.1.2. Método do Gradiente Conjugado Este método, também de 1aordem, combina (linearmente) o gradiente da iteração atual com o da anterior definindo a nova direção de busca, daí o nome gradiente conjugado. Também tem taxa de convergência quadrática e tem um desempenho superior ao método do gradiente descendente sem aumentar significativamente o esforço computacional. A demonstração desta última afirmação é encontrada em Himmelblau e Edgar (1980), página 219. O algoritmo do método do gradiente conjugado pode ser assim resumido: P1. Defina uma estimativa inicial xo P2. Calcule f(xo) e ∇f(xo) , faça so = -∇f(xo) Ricardo Kalid - Otimização de Processos – kalid@ufba.br 97 P3. Compute x1 = xo + λo.so . λo é o escalar que minimiza f(x) na direção so , ou seja é obtido através de uma busca unidimensional. P4. Calcule f(x1) e ∇f(x1) . A nova direção de busca é uma combinação linear entre so e ∇f(x1) : ( ) ( )[ ] ( )[ ]( )[ ] ( )[ ]00 11 011 xfxf xfxfsxfs T T ∇∇ ∇∇ +−∇= Equação 5-G Para a k-ésima iteração: ( ) ( )[ ] ( )[ ]( )[ ] ( )[ ]kTk kTk kkk xfxf xfxfsxfs ∇∇ ∇∇ +−∇= ++ ++ 11 11 Equação 5-H P5. Teste se a FO esta sendo minimizada, se não volte ao passo P4. P6. Termine o algoritmo quando ||sk|| for menor que a tolerância desejada. Um problema com este método é que se a razão entre os produtos internos dos gradientes entre as iterações k+1 e k for pequena recaimos no algoritmo do gradiente descendente. 5.1.3. Método de Newton Este é um método indireto de 2a ordem. Vamos aproximar a FO por uma série de Taylor truncada no 3o termo. ( ) ( ) ( )[ ] ( ) ( )( )kkTkkTkk xxHxxxfxfxf ΔΔ+Δ∇+≅ .. 2 1 Equação 5-I O mínimo de f(x) na direção xk é obtido derivando a equação 5-i em relação a x e igualando a zero ( ) ( )( ) 0. =Δ+∇ kkk xxHxf Equação 5-J ou ( )[ ] ( )kkkkk xfxHxxx ∇−=−=Δ −+ 11 Equação 5-K onde [H(xk)]-1 é o inverso da matriz Hessiana. Observe que a equação 5-k define simultaneamente o sentido e a magnitude {-[H(xk)]-1∇f(xk)} do passo de atualização (Δxk) das variáveis de decisão (xk) . Se f(x) for realmente quadrática então o mínimo é alcançada em uma iteração. Entretanto, em geral para uma FO não-linear são necessários várias