OTIMIZAÇÃO MULTIOBJETIVO NO CONTROLE DA DENGUE Bettina Fiorini Bannwart Dissertação apresentada à Universidade Es- tadual Paulista “Júlio de Mesquita Filho” para a obtenção do t́ıtulo de Mestre em Bio- metria. BOTUCATU São Paulo - Brasil Dezembro - 2013 OTIMIZAÇÃO MULTIOBJETIVO NO CONTROLE DA DENGUE Bettina Fiorini Bannwart Orientadora: Profa. Dra. Helenice de Oliveira Florentino Silva Co-orientadora: Profa. Dra. Daniela Renata Cantane Dissertação apresentada à Universidade Es- tadual Paulista “Júlio de Mesquita Filho” para a obtenção do t́ıtulo de Mestre em Bio- metria. BOTUCATU São Paulo - Brasil Dezembro - 2013 ii Dedicatória Aos meus queridos e amados pais, Dimas e Sandra. Agradecimentos Agradeço a Deus pelo dom da vida, pela Sua grandeza ao me mostrar que sou protegida, guiada e iluminada. “Agradeço por me dar abrigo na tempestade, por endireitar o que está torto, por criar sáıdas onde parece não haver escapatória. Agradeço, Senhor, pela Sua compaixão, pela Sua graça, pela Sua bondade, que estão sempre presentes, me sustentando nos momentos mais dif́ıceis.” Agradeço aos meus pais, Sandra e Dimas, por todo amor, compreensão, força, amizade, apoio, incentivo e conforto que sempre me proporcionaram. Também às minhas irmãs, Ĺıvia e Paloma, que são meus maiores tesouros, minhas amigas, minha base. À minha querida vó/mãe Zoraide que de tudo fez por mim durante esta jornada e por quem serei eternamente grata. Agradeço à querida profa Dra Helenice de Oliveira Florentino Silva pela orientação ı́mpar destinada a mim, pelos conselhos e apoio também na vida pessoal, pela atenção, pela disponibilidade e principalmente pela amizade que constrúımos durante o desenvolvimento deste trabalho, alguém por quem sempre terei admiração. Agradeço à profa Dra Daniela Renata Cantane pela co-orientação e pelas valiośıssimas dicas e sugestões. Aos professores e funcionários do departamento de Bioestat́ıstica, em especial ao Marcos e Arthur pela amizade e imensa colaboração. Às amizades que a Biometria me proporcionou: querida Ĺıvia, Luiz e Angelo por toda ajuda e apoio destinados a mim. Agradeço ao Renan por todo apoio, atenção, carinho, alegrias e con- selhos que, apesar de exagerados, contribuiram para esta formação e com certeza, contribuirão para as próximas. Obrigada pela amizade verdadeira e por tudo o que v fez por mim e pela nossa amizade. À Fapesp e à CAPES pelo apoio financeiro. Agradeço à todos que de alguma forma contribuiram e apoiaram as minhas escolhas e decisões, pela paciência e pelo est́ımulo, aos que estão perto e aos que estão longe. Sumário Página LISTA DE FIGURAS viii LISTA DE TABELAS x RESUMO xi SUMMARY xiii 1 INTRODUÇÃO 1 2 DENGUE 3 2.1 A transmissão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Sorotipos do v́ırus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Ciclo de vida do mosquito . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 Controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4.1 Controle biológico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4.2 Controle qúımico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.3 Controle mecânico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.4 Controle genético . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 MODELO MATEMÁTICO 12 4 OTIMIZAÇÃO MULTIOBJETIVO 20 4.1 Modelo Multiobjetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2 Problema de Otimização Multiobjetivo (MOP) . . . . . . . . . . . . . . . 22 vii 4.3 Dominância . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5 ALGORITMO GENÉTICO MULTIOBJETIVO 25 5.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.2 População Inicial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.3 Avaliação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.4 Elite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.5 Seleção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.6 Crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.7 Mutação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.8 Migração . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.9 Atualização, Finalização e Critério de parada . . . . . . . . . . . . . . . 36 6 EXPERIMENTOS COMPUTACIONAIS 39 7 CONCLUSÃO 49 REFERÊNCIAS BIBLIOGRÁFICAS 51 8 APÊNDICE 55 A MÉTODOS NUMÉRICOS 55 A.1 Método Numérico de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . 55 A.1.1 Runge-Kutta de ordem 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 57 A.1.2 Método de Runge-Kutta de ordem 4 para Sistemas de Equações Dife- renciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 A.2 Integração Numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A.2.1 Regra 1/3 de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 A.2.2 Regra 1/3 de Simpson Generalizada . . . . . . . . . . . . . . . . . . . 61 Lista de Figuras Página 1 Disseminação da dengue em escala mundial em 2007. . . . . . . . . . . . 3 2 Ciclo de vida do Aedes Aegypti. . . . . . . . . . . . . . . . . . . . . . . . 6 3 Dinâmica da população de mosquitos sem o uso de controle. . . . . . . . 13 4 Dinâmica da população de mosquitos com a inserção de mosquitos machos estéreis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 5 Dinâmica da população de mosquitos utilizando controle qúımico u1 e controle genético u2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 6 Fluxograma representativo de uma estrutura simples do AG. . . . . . . . 27 7 Estrutura do cromossomo definida por matriz de 2 linhas e tf + 1 colunas. 28 8 Nı́vel de dominância dos indiv́ıduos representados no espaço objetivo. . . 30 9 Ilustração do método da Roleta Viciada . . . . . . . . . . . . . . . . . . 31 10 Probabilidade acumulada usando dados do gráfico da Figura 9. . . . . . . 32 11 Crossover de 1 ponto de posição aleatória. . . . . . . . . . . . . . . . . . 33 12 Crossover de 2 pontos de posição aleatória. . . . . . . . . . . . . . . . . 33 13 Crossover uniforme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 14 Ilustração de como o operador Crossover age sobre os indiv́ıduos. . . . . 34 15 Exemplo de um processo de mutação sofrida em dois genes de cada linha da matriz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 16 Conjunto de valores de u1 otimizados em que cada linha é representada por uma solução do problema. . . . . . . . . . . . . . . . . . . . . . . . . 41 17 Conjunto de valores de u2 otimizados em que cada linha é representada por uma solução do problema. . . . . . . . . . . . . . . . . . . . . . . . . 42 ix 18 Curva de Pareto obtida pelo Algoritmo Genético Multiobjetivo. . . . . . 43 19 Variável de controle u1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 20 Variável de controle u2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 21 Evolução da população na Fase Aquática A. . . . . . . . . . . . . . . . . 45 22 Evolução da população de Fêmeas Imaturas I. . . . . . . . . . . . . . . . 46 23 Evolução da população de Fêmeas Fertilizadas F. . . . . . . . . . . . . . 46 24 Evolução da população de Machos Naturais M. . . . . . . . . . . . . . . 47 25 Evolução da população de Machos Estéreis S. . . . . . . . . . . . . . . . 47 26 Área S relacionada à integração I . . . . . . . . . . . . . . . . . . . . . . 60 Lista de Tabelas Página 1 Vocabulário do AG. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2 Pseudocódigo do AG multiobjetivo. . . . . . . . . . . . . . . . . . . . . . 37 3 Parâmetros utilizados no sistema de otimalidade (6). . . . . . . . . . . . 39 4 Parâmetros utilizados no AG. . . . . . . . . . . . . . . . . . . . . . . . . 40 5 Pseudocódigo do método RK de ordem 4. . . . . . . . . . . . . . . . . . 58 OTIMIZAÇÃO MULTIOBJETIVO NO CONTROLE DA DENGUE Autora: BETTINA FIORINI BANNWART Orientadora: Profa. Dra. HELENICE DE OLIVEIRA FLORENTINO SILVA Co-orientadora: Profa. Dra. DANIELA RENATA CANTANE RESUMO A dengue é uma doença infecciosa febril aguda causada por um v́ırus da famı́lia Flaviridae e é transmitida através de mosquito, comumente do gênero Aedes Aegypti. Esta doença tem sido atualmente um problema mundial de saúde pública, pois em 45 dias de vida um único mosquito pode contaminar até 300 pessoas e existe uma estimativa da Organização Mundial da Saúde (OMS) que anualmente de 50 a 100 milhões de pessoas se infectam em mais de 100 páıses em quase todos os conti- nentes, cerca de 550 mil doentes necessitam de hospitalização e 20 mil morrem em conseqüência da dengue. Desta forma, atualmente a dengue é um assunto de intensa pesquisa, seja na busca de vacinas e tratamentos para a doença ou nas formas eficien- tes e econômicas de controle do mosquito. Assim, este trabalho apresenta um estudo dos processos biológicos para formulação de modelos matemáticos que descrevem a dinâmica populacional do mosquito transmissor da dengue, visando a investigação de xii controles otimizados destes mosquitos. É proposta a utilização de controles qúımico com uso inseticida e genético com liberação de machos estéreis no ambiente natu- ral. Tal problema de controle ótimo visa a minimização dos investimentos com o inseticida e com a produção de machos estéreis, minimizando também a quantidade de fêmeas fertilizadas e o efeito do inseticida sobre os machos estéreis inseridos na população. É proposto um algoritmo genético para resolução do modelo de controle ótimo aplicado a problemas de combate ao mosquito da dengue e são discutidos os resultados computacionais obtidos. MULTIOBJECTIVE OPTIMIZATION IN DENGUE CONTROL Author: BETTINA FIORINI BANNWART Adviser: Prof. Dr. HELENICE DE OLIVEIRA FLORENTINO SILVA Co-adviser: Prof. Dr. DANIELA RENATA CANTANE SUMMARY Dengue is an febrile infectious disease caused by a virus of the Flavi- ridae family and is transmitted throught mosquito, the Aedes Aegypti genus is the most common. This disease has been currently a worldwide problem of public health, because about 45 days of life a single mosquito can infect up to 300 people and is an estimate of the World Health Organization (WHO) annually from 50 to 100 million people are infected in over 100 countries from all continents, about 550 thousand patients require hospitalization and 20 thousand die as a result of dengue. Thus, dengue is currently a subject of intense research, whether in the search for vaccines and treatments for the disease or the efficient and economical forms of mosquito control. Thus, this paper presents a study of biological processes for formulation of mathematical models that describe the population dynamics of the mosquito that transmits dengue, aimed at investigating optimized controls these mosquitoes. It is xiv proposed to use controls with chemical insecticide use and genetic-releasing sterile males into the natural environment.Such optimal control problem aims at minimi- zing the investments of the insecticide and the production of male-sterile, while also minimizing the amount of fertilized females and the effect of the insecticide on the inserted male-sterile population. We propose a genetic algorithm to solve the model applied to optimal control problems to combat dengue mosquito and discusses the computational results. 1 INTRODUÇÃO Atualmente considerada a arbovirose mais comum que atinge o ho- mem, a dengue tem sido motivo de diversas pesquisas no que diz respeito à criação de vacinas, tratamentos eficazes e métodos mais econômicos de controle do mosquito transmissor. O v́ırus é principalmente transmitido pelo mosquito Aedes Aegypti especialmente nas regiões tropicais e subtropicais do mundo onde as condições ambi- entais favorecem o desenvolvimento e a proliferação do mesmo. Durante os peŕıodos quentes e chuvosos o número de mosquitos no ambiente cresce, resultando assim no aumento de casos de dengue, (Thomé, 2007). Ainda não existe vacina para a dengue e uma forma eficaz de trata- mento da doença ainda permanece em estudos. Assim, o mais importante neste momento é aperfeiçoar as formas de controle do mosquito transmissor. Os métodos de controle do mosquito mais conhecidos e abordados atualmente são: controle mecânico (em geral, visitas de agentes de saúde nas re- sidências), controle qúımico (aplicação de inseticidas), controle biológico (uso de predadores e larvicidas biológicos) e controle genético, cuja técnica baseia-se na in- trodução de mosquitos machos estéreis na população visando reduzir a fertilidade da mesma. Esteva & Yang (2006) propuseram um modelo matemático que des- creve a dinâmica da população de mosquitos (dividida em fase aquática e alada) utilizando o controle genético. Thomé (2007) adaptou este modelo de forma que, além do controle genético, foi inserido também o controle qúımico aplicado essen- cialmente na fase adulta do mosquito. Além disso, Thomé (2007) apresentou um funcional que analisa o ı́ndice de performance do sistema, criando um problema de 2 otimização que visa diminuir o custo dependido com os controles. Os objetivos foram: diminuir a quantidade de fêmeas fertilizadas, pre- servar a quantidade de machos estéreis inseridos na população e diminuir o custo com os controles aplicados no sistema, ou seja, minimizar o custo com inseticidas e machos estéreis. O presente trabalho propõe, assim, uma abordagem multiobjetivo para o problema de otimização proposto por Thomé (2007) e propõe um algoritmo genético para resolvê-lo. O trabalho está divido da seguinte forma: no Caṕıtulo 2 estão descritas as caracteŕısticas do mosquito transmissor e da dengue, bem como sua forma de transmissão e métodos de controle. No Caṕıtulo 3 são descritos modelos matemáticos que descrevem a dinâmica da população dos mosquitos transmissores da dengue na presença do controle qúımico com inseticida e controle genético com a inserção de mosquitos machos estéreis propostos pelos autores Esteva & Yang (2006) e Thomé (2007). No Caṕıtulo 4 é proposto um modelo de otimização multiobjetivo para minimização dos investimentos com controles qúımico e genético, minimização do número de fêmeas fertilizadas no sistema e maximização a quantidade de machos estéreis que permanece no ambiente (ou minimização do efeito do inseticida sobre os machos estéreis). No Caṕıtulo 5 é proposto um algoritmo genético multiobjetivo para resolução do modelo de otimização apresentado no Caṕıtulo 4. No Caṕıtulos 6 são apresentados os resultados computacionais e no Caṕıtulo 7 são apresentadas as conclusões obtidas no trabalho. 2 DENGUE A dengue é uma doença infecciosa causada por um arbov́ırus cujo desenvolvimento e proliferação são favorecidos pelas condições ambientais de páıses de climas tropicais e subtropicais, como ilutrado na Figura 1. Fonte: Adaptado de Dengue (2013) Figura 1 - Disseminação da dengue em escala mundial em 2007. Segundo a Organização Mundial da Saúde (2013), a dengue se tornou uma grande preocupação na saúde pública mundial. Estima-se que 2,5 bilhões de pessoas habitam em mais de 100 páıses endêmicos e áreas suscept́ıveis à transmissão 4 da dengue. A distribuição geográfica dos mosquitos vetores e dos v́ırus levou ao ressurgimento global da epidemia de dengue nos últimos 25 anos com o desenvolvi- mento da hiperendemicidade em muitos centros urbanos dos trópicos, (Organização Mundial da Saúde, 2013). Diversos fatores contribúıram para produzir condições epidemiológicas dos páıses endêmicos: o rápido crescimento populacional, a migração da zona rural para a região urbana, a infra-estrutura básica inadequada e o aumento de reśıduos sólidos que proporcionam habitats para as larvas em áreas urbanas. No Brasil, a disseminação da doença acontece principalmente no verão até o mês de maio, quando há maior incidência de chuvas. Até o ano de 2005, a epidemiologia da dengue apresentou dois padrões distintos: ondas epidêmicas loca- lizadas em grandes centros urbanos (1986 a 1993) e epidemias/circulação endêmica do v́ırus em todas as regiões do páıs (1994 a 2005). Já no ano de 2006 marcou o ińıcio de um novo cenário na epidemiologia da doença, caracterizada pela migração de gravidade para crianças, com ápice nas epidemias registradas no ano de 2008, em especial no Estado do Rio de Janeiro, (Teixeira et al., 2008). Hoje há uma grande preocupação com a dengue hemorrágica, pois de acordo com a Organização Mundial da Saúde (2013), em sua forma mais severa a dengue hemorrágica pode ser uma complicação potencialmente fatal devido ao vazamento de plasma, acúmulo de ĺıquido, dificuldade respiratória, hemorragia grave ou ataque a órgãos como f́ıgado, timo, baço e gânglios linfáticos. 2.1 A transmissão O Aedes Aegypti é um transmissor do v́ırus com importância epide- miológica. Outro transmissor em potencial é o mosquito Aedes Albopictus, no en- tanto, a transmissão por este vetor é menos comum já que ele não é doméstico e prefere os ocos de árvores para depositar seus ovos, possui hábitos antropof́ılicos1 e 1Adaptado para parasitar ou infectar os seres humanos. 5 zoof́ılicos2 diurnos e não costuma frequentar domićılios. Já o Aedes Aegypti é um mosquito essencialmente doméstico, vive em locais bastante frequentados por pessoas e possui hábitos preferencialmente diurnos, (Teixeira et al., 1999). Segundo Thomé (2007), a transmissão da doença nos seres humanos é causada pela fêmea adulta do mosquito através do ciclo homem-Aedes Aegypti - homem e isso ocorre porque apenas elas são hematófagas (alimentam-se de sangue). A fêmea pica uma pessoa infectada, conserva o v́ırus cujo peŕıodo de incubação é de 8 a 12 dias e o retransmite para outras pessoas em todas as picadas que realizar durante sua vida. Copulando com o macho uma vez a fêmea já está apta a realizar entre 150 e 200 posturas de ovos em diversos criadouros, garantindo a dispersão e a preservação da espécie, (Instituto Oswaldo Cruz, 2013). Em 45 dias de vida, um único mosquito é capaz de contaminar até 300 pessoas e o peŕıodo de incubação do v́ırus no ser humano varia de 3 a 15 dias após a picada do mosquito. Não há transmissão por contato direto de um doente ou de suas secreções para uma pessoa sadia, nem através da água ou alimento. No entanto, há relatos de transmissão através de doações de sangue e órgãos e também através da gestação, (Organização Mundial da Saúde, 2009). 2.2 Sorotipos do v́ırus De acordo com Marinho (2013) e Organização Mundial da Saúde (2013), existem quatro sorotipos diferentes do v́ırus da dengue, denominados DEN-1, DEN-2, DEN-3 e DEN-4. Estes sorotipos pertencem ao gênero Flavivirus e são so- rologicamente relacionados, porém antigenicamente distintos e um não concede imu- nidade ao outro, apenas a ele mesmo, assim pessoas que vivem em áreas endêmicas podem contrair a doença até quatro vezes, (Thomé, 2007; Teixeira et al., 1999). 2Infectam preferencialmente os animais. 6 2.3 Ciclo de vida do mosquito Assim como ocorre com outros mosquitos, o ciclo do Aedes aegypti passa por uma metamorfose completa, composta pelas fases: ovo, larva, pupa e adulto, Figura 2. Pode-se dizer também que este ciclo é dividido em dois estágios: fase aquática (ovos, larvas e pupas) e fase alada (mosquitos adultos). Fonte: Adaptado de Casa das Ciências (2013) Figura 2 - Ciclo de vida do Aedes Aegypti. Os mosquitos machos e fêmeas se acasalam logo após emergirem dos 7 casulos de pupa. Sendo que os machos vivem poucos dias após o cruzamento e alimentam-se de néctar de plantas até morrerem. Já as fêmeas fertilizadas alimentam-se de sangue até que seus ovos amadureçam e, em seguida, utilizam depósitos de águas naturais ou artificiais para a oviposição, preferencialmente locais escuros e não polúıdos. Os ovos ainda podem sobreviver vários meses na ocorrência de dessecação e, quando úmidos, eclodem em larvas de 2 a 3 dias conforme a variação da temperatura e do clima. As larvas possuem quatro estágios (́ınstares3) e durante esse peŕıodo se alimentam de algas e microorganismos filtrados da água. Os estágios larvais passam rapidamente com exceção do último no qual a larva permanece por mais tempo até alcançar a fase pupal. Nesta fase, as pupas não se alimentam e levam até três dias para amadurecerem. Logo após adquirirem a forma alada, os mosquitos estão prontos para acasalarem e iniciar o ciclo novamente. 2.4 Controle A incidência crescente e amplo alcance geográfico da dengue fazem com que o desenvolvimento de um combate eficaz contra essa doença seja considerado uma prioridade de saúde internacional, (Guy et al., 2011). De acordo com a Organização Mundial da Saúde (2013), não há vacinas preventivas contra a dengue, no entanto estas estão sendo pesquisadas em âmbito mundial e, embora tenha havido recentes progressos e ocorrência de testes em todo mundo, sua conclusão ainda é um desafio. Dessa forma, a OMS considera o controle e a prevenção, atualmente, as melhores formas de combate à doença. Existem muitas medidas preventivas para impedir que a dengue seja disseminada da população. Com sucesso ou não, elas vêm sendo intensamente aplicadas e/ou estudadas. São exemplos de algumas des- tas medidas: sistema de alerta de surtos e casos, capacitação para o manejo dos casos cĺınicos, mobilização social, controles f́ısico (ou mecânico), genético, qúımico e biológico do mosquito e muitos outros. 3Estágio larval de alguns artrópodes atingido após uma muda ou ecdise. 8 Para controle da doença, têm sido estudados programas de controle do Aedes Aegypti. Atualmente os controles f́ısicos e qúımicos são os mais intensamente utilizados na prática. Os meios mais utilizados para o controle f́ısico são a diminuição dos criadouros e o alerta da população, e do qúımico é a aplicação de inseticidas. Segundo Donaĺısio & Glasser (2002), o desenvolvimento de resistência de mosquitos aos inseticidas é um processo inevitável. A exposição a dosagens matam os indiv́ıduos suscet́ıveis e os resistentes sobrevivem e transferem essa capacidade a seus descendentes. Além dos vários mecanismos de resistência presentes nos insetos, que permitem sua sobrevivência após contato com o inseticida, há outra forma cha- mada de comportamental, que define o processo de seleção de indiv́ıduos com aptidão para evitar total ou parcialmente o contato com doses que seriam letais, (Donaĺısio & Glasser, 2002). Assim, nas últimas décadas, vem sendo reiterada a recomendação do controle integrado do Aedes Aegypti com implementação descentralizada, envolvendo o poder público e a sociedade. Esse tipo de estratégia teria maior sustentabilidade que aquelas verticais centralizadas e baseadas em um único método. No controle integrado do Aedes Aegypti, as medidas preventivas são direcionadas principalmente aos criadouros, constituindo-se de ações simples e eficazes, especialmente aquelas que consistem em cuidados a serem adotados pela população. A tecnologia atual- mente em estudo abrange medidas de controle f́ısico, genético, qúımico e biológico, (Donaĺısio & Glasser, 2002). 2.4.1 Controle biológico Dentre as diversas medidas de controle biológico, citamos o uso de predadores e de larvicidas biológicos. Os predadores do tipo peixes larvófagos são os mais recomendados por sua fácil obtenção e manutenção, especialmente para bebedouros de grandes animais, fossos de elevador de obras, espelhos d’água/fontes ornamentais, piscinas abandonadas e depósitos de água não potável. Segundo a EMBRAPA (2010), o 9 peixe Barrigudinho, conhecido também como Lebiste, Guaru ou Guppy, possuindo apenas quatro cent́ımetros pertencente à famı́lia Poecilidae é uma posśıvel arma na guerra biológica para o controle do mosquito Aedes Aegypti, transmissor da dengue, (EMBRAPA, 2010). O peixe em questão resiste à variações de ambiente, tem grande potencial de reprodução facilitando a disseminação e sobrevivência da espécie. O Bacillus thuringiensis israelensis, Bti, é um larvicida constitúıdo de uma bactéria patogênica a um grande número de pragas. A pesquisa foi desenvolvida em Israel, 1978, com o patroćınio da OMS, e muitos produtos com base nesta bactéria tem sido utilizados em programas de controle de moscas e mosquitos, (Ministério da Saúde, 2013). Ainda de acordo com informações do Ministério da Saúde (2013), existe uma rede de monitoramento que avalia o estágio de resistência do Aedes Aegypti ao uso de inseticidas. Ao ser detectada a resistência ao uso de organofosforados no munićıpio, desencadeia-se o processo de substituição pelo Bti. 2.4.2 Controle qúımico No Brasil, o controle qúımico é um dos mais utilizados, usando inse- ticidas para controle da população de mosquitos adultos, levando-os a morte. No entanto, o uso intenso de inseticidas eleva o custo do controle, pode afetar a saúde pública e desenvolver resistência nos mosquitos, (Donaĺısio & Glasser, 2002). Além disso, é indispensével o uso racional e seguro desses larvicidas já que podem causar grandes impactos ambientais. Os inseticidas podem atingir diferentemente as fases do mosquito: o larvicida atinge as fases de larva e pupa do mosquito, enquanto o adulticida é responsável pela indução da mortalidade do adulto através de equi- pamentos portáteis (dentro das casa) e de equipamentos pesados (pulverização nas ruas). 2.4.3 Controle mecânico Segundo o Ministério da Saúde (2013), o controle mecânico consiste em práticas capazes de impedir a procriação do mosquito transmissor, tendo como prin- 10 cipais providências a proteção, a destruição ou a destinação adequada de criadouros, que devem ser executadas sob a supervisão do ACE (agente de combate de endemias) ou ACS (agende comunitário de saúde), prioritariamente pelo próprio morador e/ou proprietário da residência. Esse controle é feito através de coleta de pneumáticos, de reśıduos sólidos proṕıcios ao aparecimento de criadouros com destino posterior adequado e de vedação de depósitos armazenadores de água, (Ministério da Saúde, 2013). 2.4.4 Controle genético Quanto aos métodos de controle genéticos, a utilização de machos estéreis visando reduzir a fertilidade da população local foi usada em testes em campo e tem sido objeto de muitos estudos, apresentando resultados promissores. Segundo Thomé (2007), insetos estéreis são liberados no ambiente na- tural de tal maneira que os acasalamentos resultem na produção de ovos inviáveis, ou seja, não serão capazes de atingir a fase adulta, o que deve reduzir a população do mosquito a tal ńıvel que controle a transmissão da dengue. Os insetos se tornam estéreis devido ao uso de agentes que causam mutações tais como a radiação gama. Esta técnica, conhecida como Sterile Insect Technique (SIT), foi de- senvolvida pelo entomólogo americano Edward Knipling em meados da década de 50, e tem se mostrado muito eficiente no controle de pragas agŕıcolas, aparecendo como uma medida alternativa para a técnica usual de aplicação de inseticida, que além de promover resistência do inseto ao produto qúımico, é não seletiva, ou seja, prejudica outras espécies que vivem no mesmo habitat do mosquito, (Thomé, 2007). Devido à distribuição espacial heterogênea dos criadouros do vetor da dengue, há dificuldade em aplicar na prática este tipo de controle genético, visto que quando as indústrias fabricam estes indiv́ıduos em grande escala, elas precisam libertá-los imediatamente na natureza, devido ao curto tempo de vida do vetor, (Thomé, 2007). No Brasil, em 2005, foi realizada com sucesso a liberação de moscas estéreis com objetivo de combater a mosca-das-frutas Ceratitis Capitata, conhecida 11 em todo o mundo como mosca-do-mediterrâneo ou moscamed. Os custos por hectare para utilização da moscamed esterelizada é em média 20% superior do que o uso do agrotóxico e a liberação dos machos estéreis é feita através de aeromodelos monito- rados por controle remoto. A pesquisa com o Aedes Aegypti transgênico começou no páıs em 2010, com a adaptação do mosquito em laboratório da Universidade de São Paulo. O Projeto Aedes Transgênico (PAT) realizado em parceria com a empresa britânica Oxitec, desenvolveu a primeira linhagem do Aedes modificado. Entre os anos de 2011 e 2012, foi realizado um projeto piloto em dois bairros de Juazeiro (BA) - Mandacaru e Itaberaba - cujo ı́ndice de proliferação do mosquito é considerado alto. A técnica mostrou-se bem sucedida apresentando uma redução de 90% da população do mosquito em aproximadamente 6 meses, (Ministério da Saúde, 2012). Ainda está em estudo a aplicação desta técnica em cidades de médio e grande porte a fim de mensurar a redução da doença na população, verificar a melhor maneira de adaptação do mosquito no ambiente e adequar o transporte e a loǵıstica do mesmo. Pode-se notar que a dengue é um assunto de intensa pesquisa, seja na busca de vacinas, tratamentos eficazes ou em uma forma econômica de controle do mosquito. Estes estudos podem muitas vezes serem facilitados com o uso de modelos matemáticos, os quais podem auxiliar na obtenção de estimativas e no entendimento da dinâmica populacional dos mosquitos. Assim, o caṕıtulo 3 apresenta alguns estudos relativos a esses tipos de modelos. 3 MODELO MATEMÁTICO Esteva & Yang (2006) e Thomé (2007) propuseram modelos que ana- lisam o impacto da aplicação de inseticida e da inserção dos machos estéreis na natureza. Os machos estéreis não apresentam qualquer risco aos seres humanos, já que apenas as fêmeas podem picar e transmitir a doença. Para descrever a dinâmica da população com os controles aplicados, os autores consideram o ciclo de vida dos mosquitos divididos em fase aquática (ovos, larva, pupa) e fase alada (mosquitos adultos). Consideram também o controle genético com introdução de mosquitos machos estéreis na fase alada esterilizados pela técnica de irradiação, e a inserção de inseticidas na população Aedes Aegypti. O resultado esperado é a redução da população do inseto e/ou até o seu desaparecimento do ambiente. O modelo proposto por Esteva & Yang (2006) descreve a dinâmica dos mosquitos nas fases aquática e alada utilizando as seguintes notações para tamanhos das populações: • A(t)→ População de mosquitos na Fase Aquática, no tempo t; • I(t)→ População de mosquitos Fêmeas Imaturas no tempo t; • F (t)→ População de mosquitos Fêmeas Fertilizadas no tempo t; • U(t)→ População de mosquitos Fêmeas Não-Fertilizadas no tempo t; • M(t)→ População de mosquitos Machos Naturais, no tempo t; • S(t)→ População de mosquitos Machos Estéreis, no tempo t 13 São considerados no modelo as taxas de mortalidade dadas por µA (fase aquática), µI (fêmea imatura), µF (fêmea fertilizada), µU (fêmea não-fertilizada), µM (macho natural) e µS (macho estéril), as taxas de oviposição ϕ das fêmeas fertilizadas, C mede a capacidade do meio (nutrientes, espaço, criadouros), taxa γ com que o mosquito evolui para fase alada, taxa βS de acasalamento entre machos estéreis e fêmeas imaturas. As populações de mosquitos fêmeas imaturas I(t) e machos naturais M(t) são originadas da evolução da fase aquática para a alada a uma taxa γ a uma proporção r e (1− r) respectivamente. A taxa de encontro dos machos naturais com as fêmeas é β e do cruzamento destes surge a população de fêmeas fertilizadas F (t). Esta população alimenta a fase aquática a uma taxa ϕ(1 − A/C)F , em que ϕ é a taxa de oviposição intŕınseca e C é a capacidade do meio relacionada com o número de nutrientes, espaço, etc. A Figura 3 ilustra esta dinâmica. Fonte: Adaptado de Barsante et al. (2011) Figura 3 - Dinâmica da população de mosquitos sem o uso de controle. 14 O modelo matemático (1) descreve esta dinâmica populacional e está apresentado a seguir.  dA dt = ϕ(1− A/C)F − γA− µAA dI dt = rγA− βI − µII dF dt = βI − µFF dM dt = (1− r)γA− µMM (1) Com a inserção de mosquitos estéreis no meio ambiente a uma taxa α, cria-se uma população de machos estéreis S e a probabilidade de encontro entre fêmeas imaturas com machos naturais M torna-se igual a ( M (M + S) ) . A taxa per capita com que as fêmeas são fertilizadas é dada por ( βM (M + S) ) . Segundo Esteva & Yang (2006) a probabilidade de encontro de machos estéreis com fêmeas imaturas não depende apenas do número de mosquitos machos irradiados e é dada por ( pS (M + S) ) , em que p, 0 < p < 1, é a proporção com que os mosquitos estéreis são colocados nos locais adequados. A taxa de acasalamento efetiva dos mosquitos estéreis é dada por qβ, com 0 ≤ q ≤ 1. A taxa per capita com que as fêmeas imaturas cruzam com mosquitos estéreis S é dada por ( βSS (M + S) ) , em que βS = pqβ. Assim, a dinâmica populacional com a inserção de machos estéreis está esquematizada na Figura 4 e o modelo matemático para esta dinâmica proposto por Esteva & Yang (2006) está apresentado em (2). 15 Fonte: Adaptado de Barsante et al. (2011) Figura 4 - Dinâmica da população de mosquitos com a inserção de mosquitos machos estéreis.  dA dt = ϕ(1− A/C)F − γA− µAA dI dt = rγA− βMI M + S − βSSI M + S − µII dF dt = βMI M + S − µFF dM dt = (1− r)γA− µMM dS dt = α− µSS , (2) 16 No sistema 2 acima a equação da variável de estado U (População de Mosquitos Fêmeas Não-Fertilizadas) está desacoplada do sistema dinâmico, ou seja, as demais equações não dependem dela, e é descrita por: dU dt = βSSI M + S − µUU Thomé (2007) fez uma adaptação no modelo (2) e além do controle genético com inserção de machos estéreis, utilizaram também o controle qúımico com uso de inseticida aplicado na fase adulta do mosquito. A Figura 5 esquematiza a dinâmica de mosquitos apresentada por Thomé (2007) utilizando o controle qúımico com inseticida e o controle genético com inserção de machos estéreis representados por u1(t) e u2(t), respectivamente, cujas funções são variantes no tempo t. Fonte: Adaptado de Barsante et al. (2011) Figura 5 - Dinâmica da população de mosquitos utilizando controle qúımico u1 e controle genético u2. 17 O modelo matemático apresentado para esta dinâmica de mosquitos está apresentado em (3).  dA dt = ϕ(1− A/C)F − γA− µAA dI dt = rγA− [ βM M + S + βSS M + S + (µI + u1) ] I dF dt = βMI M + S − (µF + u1)F dM dt = (1− r)γA− (µM + u1)M dS dt = u2 − (µS + u1)S (3) As condições de equiĺıbrio para o sistema (3) foi amplamente discutida e está apresentada em Thomé (2007). Para o estudo da qualidade do controle qúımico e genético otimizado, Thomé (2007) propôs o uso do ı́ndice de performance J apresentado em (4). J = J [u1, u2] = 1 2 ∫ T 0 (c1u 2 1 + c2u 2 2 + c3F 2 − c4S 2)dt (4) Em que c1, c2, c3 e c4 representam respectivamente a importância re- lativa ao uso de inseticidas u1, importância relativa à liberação de mosquitos estéreis u2, importância relativa à presença de fêmeas fertilizadas F e importância relativa à preservação de machos estéreis S no ambiente. Os custos podem admitir um caráter econômico, social, ambiental, entre outros, dependendo do interesse no uso do modelo. A discussão está em torno do alto custo do controle com mosquitos estéreis comparado ao custo com inseticida. No entanto, o uso de inseticidas pode causar efeitos como danos à saúde e ao meio ambiente, além de desenvolver a re- sistência nos mosquitos. Dessa forma, as variáveis u1 e u2 representam os custos no 18 processo do controle do mosquito: o custo de aplicação de inseticida e de inserção de machos estéreis. O problema de controle ótimo consiste em determinar u1 e u2 do sis- tema (3) de forma que minimizem o funcional (4), ou seja minimizem o custo com inseticida e com mosquitos estéreis, minimizando também a quantidade de fêmeas fertilizadas preservando o máximo posśıvel a quantidade de mosquitos estéreis intro- duzidos. A variável de controle considerada no problema de controle ótimo u1(t) está associada ao investimento em inseticidas no instante t e u2(t) relacionada ao investimento com a produção e liberação de mosquitos estéreis no instante t. Assim, a parcela c1u 2 1 mede o custo com inseticida, c2u 2 2 mede o custo com a produção e liberação de mosquitos estéreis, c3F 2 mede o custo social, peso dado ao número de fêmeas fertilizadas e a ultima parcela −c4S2 está associada ao custo de perda de machos estéreis devido ao uso de inseticida. Os coeficientes ci, i = 1, 2, 3 e 4, funcionam como pesos de cada parcela e os quadrados nas variáveis ampliam os efeitos de grandes variações. Thomé (2007) utilizou as condições de equiĺıbrio (5), para o sistema (3), por ele estudadas como condições iniciais, visto que estas são atingidas na pior situação posśıvel do ponto de vista da dengue e o controle é assim aplicado partindo da situação onde o mosquito se alastrou na natureza.  A(0) = A0 = C(R− 1) R I(0) = I0 = rγA0 µI + β F (0) = F0 = γ + µACA0 ϕ(C − A0) M(0) = M0 = (1− r)γA0 µM S(0) = S0 = 0 (5) em que R = ϕrγβ (γ + µA)(β + µ1)µF . 19 A taxa básica de reprodutividade representada por R é uma das gran- dezas mais importantes na epidemiologia e é definida como o número médio de casos secundários decorrentes de um processo primário em uma população totalmente sus- cept́ıvel e mede o potencial máximo de reprodução para uma doença infecciosa, (Keeling & Rohani, 2008). Ainda de acordo com Keeling & Rohani (2008), se R > 1 compreende-se que a infecção atingiu seu peŕıodo endêmico. Para resolução do problema de controle ótimo, que consiste em mini- mizar o funcional (4), sujeito ao sistema (3) e às condições iniciais (5), Thomé (2007) utilizou o Prinćıpio Máximo de Pontryagin4 para determinar a formulação do con- trole ótimo u1 e u2. Desta forma ele obteve um problema de otimalidade estendido formado pelas restrições já existentes que são sistema de estado (3) e as condições iniciais (5), juntamente com um sistema adjunto, cujas variáveis adjuntas dificultam a solução numérica deste problema, pois estas apresentam comportamento instável. Assim propomos uma outra abordagem multiobjetivo para este pro- blema e um algoritmo genético para resolução deste, os quais estão discutidos nos caṕıtulos a seguir. 4Utilizado na teoria de Controle Ótimo, pode ser interpretado como um teorema de multiplica- dores de Lagrange em espaços de dimenção finita. 4 OTIMIZAÇÃO MULTIOBJETIVO 4.1 Modelo Multiobjetivo Problemas de otimização que apresentam dois ou mais objetivos são chamados problemas multiobjetivos. Neste caṕıtulo propomos uma abordagem mul- tiobjetivo para o problema de controle ótimo do mosquito da dengue, utilizando controle com inseticida e com mosquitos machos estéreis. Este problema consiste em determinar u1 = u1(t) e u2 = u2(t) relacionados, respectivamente, aos controles qúımico e genético no peŕıodo 0 ≤ t ≤ T , em que T é o peŕıodo total de controle. Tal problema está definido na forma: minimize J(u1, u2) = [J1, J2, . . . , Jk] sujeito a Sistema (3); Condições iniciais (5); Condições finais sobre algumas variáveis do sistema (3); Condições de não-negatividade: 0 ≤ u1(t) ≤ 1 e u2(t) ≥ 0, para todo 0 ≤ t ≤ T . Os funcionais J1, J2,...,Jk (k inteiro e maior ou igual a 2) estão rela- cionados com os objetivos propostos. Alguns objetivos que podem ser considerados são: 1. Minimizar o custo com inseticidas; 2. Minimizar o custo despendido com machos estéreis; 3. Minimizar a quantidade de fêmeas fertilizadas no sistema; 21 4. Maximizar o número de mosquitos estéreis que permanecem no ambiente; Observação: Podem ser considerados outros objetivos ou mais de um objetivo no mesmo funcional. Como exemplo, um problema de determinar o controle ótimo utili- zando dois objetivos seria o problema mostrado em (6): min J1 = J1[u1, u2] = 1 2 ∫ T 0 (c1u 2 1 + c2u 2 2 + c3F 2)dt max J2 = J2[u1, u2] = 1 2 ∫ T 0 (c4S 2)dt sujeito a dA dt = ϕ(1− A/C)F − γA− µAA dI dt = rγA− [ βM M + S + βSS M + S + (µI + u1) ] I dF dt = βMI M + S − (µF + u1)F dM dt = (1− r)γA− (µM + u1)M dS dt = u2 − (µS + u1)S A(0) = A0 = C(R− 1) R I(0) = I0 = rγA0 µI + β F (0) = F0 = γ + µACA0 ϕ(C − A0) M(0) = M0 = (1− r)γA0 µM S(0) = S0 = 0 em que R = ϕrγβ (γ + µA)(β + µ1)µF Não negatividade, 0 ≤ u1(t) ≤ 1 e u2(t) ≥ 0,∀t. (6) 22 O item a seguir dedica-se a mostrar os elementos básicos da teoria mul- tiobjetivo e posteriormente é proposto algoritmo genético para resolução do problema de otimização (6). 4.2 Problema de Otimização Multiobjetivo (MOP) Segundo Reeves & Beasley (1993) uma técnica heuŕıstica de otimização é inspirada em processos intuitivos que procuram boas soluções sem se comprometer com a otimalidade. As heuŕısticas ainda podem ter seu desempenho melhorado desde que combinadas com caracteŕısticas particulares do problema, assim chamadas de metaheuŕısticas. Para Freitas et al. (2009), metaheuŕısticas representam um conjunto de algoritmos heuŕısticos genéricos estudados desde 1970, baseando-se em ideias de diversas fontes para realizar a busca pela solução dos problemas de otimização. Ele não busca a melhor solução para o problema - a solução ótima - no entanto, o uso destes métodos são oportunos em se tratando de problemas com mais de uma função de satisfação: problemas multiobjetivo, como é o caso deste trabalho. Dentre as vantagens das metaheuŕısticas, destaca-se a facilidade na implementação além de possúırem ótima performance em problemas combinatoriais e fornecerem soluções de ótima qualidade em problemas de alto grau de dificuldade à um baixo custo operacional. Dentre as várias metaheuŕısticas difundidas, temos Simulated Annealing (SA), Colônias de Formigas, Colônia de Abelhas, Busca Tabu e o Algoritmo Genético (AG) - ou Genetic Algorithm. Um Problema de Otimização Multiobjetivo (MOP-Multi-objective Op- timization Problem), como o próprio nome sugere é um problema que envolve mais de uma função objetivo a ser minimizada ou maximizada simultaneamente. Segundo Deb & Kalyanmoy (2001), um MOP pode ser definido como um problema de procurar um vetor de variáveis de decisão que satisfaz certas res- trições e otimiza um vetor cujos elementos são funções objetivos. Tais funções for- necem objetivos usualmente conflitantes. 23 Considerando n funções-objetivo que formam o vetor f(x) = [f1(x), f2(x), . . . , fn(x)] T , a forma geral do problema de otimização multiobjetivo é: min/max f(x) sujeito a gk(x) ≥ 0, k = 1, . . . , K; hm(x) = 0,m = 1, . . . ,M ; x (L) i ≤ xi ≤ x (U) i , sendo x um vetor de variáveis de decisão tal que x = [x1, x2, . . . , xn] T representando a solução, K e M são, respectivamente, o número de restrições de desigualdade e de igualdade do problema, x (L) i e x (U) i representam os limites inferior e superior da variável xi. Neste tipo de problema, há duas classes posśıveis de soluções: • As que apresentam desempenhos ruins comparadas às demais, sob os objetivos simultaneamente considerados; e • As que, quando comparadas às demais, são melhores em um ou mais objetivos, também chamadas soluções eficientes. De acordo com Coello (2010), raramente acontece o caso de um único ponto otimizar simultaneamente todos os objetivos. Portanto, a noção de otimização adotada neste caso é o proposto por Francis Ysidro Edgeworth em 1881 e generalizado por Vilfredo Pareto em 1896, dáı o nome Pareto-Ótimo. Basicamente, emprega-se a relação de dominância de Pareto para com- parar as soluções fact́ıveis do problema. Para isso, definimos: 4.3 Dominância Definição: A solução x1 domina a solução x2 se ocorrer (considerando um problema de minimização): 24 • fm(x1) ≤ fm(x2) em todas as funções-objetivo, (∀m, m = 1, 2, . . . ,M). • fm(x1) < fm(x2) em pelo menos uma função-objetivo, (para algum m ∈ 1, . . . ,M). Uma solução que não é dominada por nenhuma outra é chamada solução não dominada ou solução eficiente formando o conjunto Pareto-Ótimo. Já a Fronteira de Pareto é o conjunto dos valores das funções objetivos cujas soluções pertencem ao conjunto Pareto-Ótimo. A busca pela solução de problemas multiobjetivos tem gerado estu- dos importantes dentro da Pesquisa Operacional. Segundo Coello (2010), apesar da grande quantidade de métodos de programação dispońıveis na literatura para a solução de MOP, a aplicabilidade tende a ser limitada (ao se deparar com funções objetivas diferenciáveis ou então com a Fronteira de Pareto), dáı vem a motivação para o uso de abordagens alternativas, como os algoritmos evolucionários. Estes algoritmos são considerados uma boa opção pois adotam uma po- pulação de soluções, o que permite encontrar elementos do conjunto Pareto-Ótimo em uma única geração enquanto em outros métodos, uma única solução não-dominada é gerada por vez, (Coello, 2010). Segundo Coello (2010), os algoritmos evolucionários tendem a ser menos suscept́ıveis à descontinuidade, adquirindo uma vantagem sig- nificativa frente aos demais métodos de programação. 5 ALGORITMO GENÉTICO MULTIOBJE- TIVO 5.1 Introdução Neste caṕıtulo propomos um algoritmo genético multiobjetivo para a solução do problema (6). Goldberg (1953) define Algoritmo Genético como um algoritmo evolu- tivo de busca baseado em mecanismos de seleção e genética naturais, ou seja, uma metaheuŕıstica que possui analogia com a natureza. A ideia central do algoritmo é explorar informações históricas sobre pontos antigos sobre os novos a fim de melhorar o desempenho da busca. Isso é feito através de processos iterativos, em que cada iteração é chamada geração. O AG, desenvolvido pioneiramente nos trabalhos de John Henry Hol- land (década de 70), inspira-se na teoria de evolução das espécies de Charles Darwin, privilegiando a sobrevivência dos indiv́ıduos mais aptos procurando preservar, no processo de seleção de espécies, aqueles que possuem caracteŕısticas mais desejáveis. Nos anos 80, David Goldberg, aluno de Holland, consegue o primeiro sucesso em aplicação industrial de Algoritmos Genéticos. Desde então, estes algoritmos vêm sendo aplicados com sucesso nos mais diversos problemas de otimização. De acordo com a teoria Darwinista, “os que sobrevivem e geram des- cendentes são aqueles selecionados e adaptados ao meio devido às relações com os de sua espécie e também ao ambiente onde vivem. A cada geração, a seleção natural favorece a permanência das caracteŕısticas adaptadas, constantemente aprimoradas e constantemente melhoradas. É a evolução das espécies”. O AG busca exatamente 26 a simulação desde processo seletivo que ocorre na natureza: a reprodução, diversi- ficação e a sobrevivência do mais apto, (Aliano Filho, 2012). O método diferencia-se dos demais algoritmos devido a quatro carac- teŕısticas: • Trabalha com múltiplas soluções; • Busca uma população de pontos, e não um único ponto; • Utiliza informações de recompensa - função objetivo - e não derivadas, con- tinuidade, convexidade ou outro conhecimento auxiliar. Isso faz do AG um método mais canônico do que muitos sistemas de busca; • Usa regras de transição probabiĺısticas e não determińısticas. Tais caracteŕısticas, em conjunto, conferem robustez ao algoritmo e consequentemente vantagem sobre outras técnicas comumente usadas, (Goldberg, 1953). Além disso, o AG constitui um método de busca aleatória funcionando melhor que heuŕısticas comuns quando se dispõe de pouca informação sobre o espaço de busca. De acordo com Goldberg (1953), um algoritmo genético simples que confere bons resultados na prática é composto pelos operadores: seleção, crossover (cruzamento) e mutação que fazem a solução evoluir no sentido da otimização. No caso do crossover, os indiv́ıduos escolhidos cruzam dando origem a novos indiv́ıduos através da recombinação de partes dos cromossomos dos pais. E, por fim, a mutação confere diversidade aos indiv́ıduos prevenindo a convergência para ótimos locais. O vocabulário a seguir, Tabela 1, descrito por Goldberg (1953) será necessário para na compreensão do método: 27 Tabela 1. Vocabulário do AG. Linguagem Natural AG Indiv́ıduo Solução População Conjunto de indiv́ıduos Geração Iteração Cromossomo Estrutura que representa a solução Gene Unidade do cromossomo Alelo Valor do gene Aptidão do indiv́ıduo Valor da função objetivo de cada indiv́ıduo Goldberg (1953) considera o algoritmo dividido em partes, como vemos a seguir na Figura 6. Figura 6 - Fluxograma representativo de uma estrutura simples do AG. 28 5.2 População Inicial Os indiv́ıduos, cromossomos ou soluções na população inicial podem ser gerados construtiva ou aleatoriamente. Para isto é necessário definir o tamanho da população e uma estrutura para estes cromossomos. Normalmente, utiliza-se a ale- atoriedade para gerar a população garantindo a diversidade genética dos indiv́ıduos e abrangendo um maior espaço de solução. A forma construtiva é usada quando se almeja uma melhora no desempenho do método, pois ela garante o recebimento de indiv́ıduos com caracteŕısticas desejadas, (Aliano Filho, 2012). De acordo com Aliano Filho (2012), o tamanho da população é esco- lhido de acordo com a necessidade do problema, no entanto, um baixo número de indiv́ıduos pode convergir a um ótimo local e, por outro lado, muitos indiv́ıduos fazem com que o AG perca sua eficiência computacional. Para resolução do problema de otimalidade (6) propomos um algoritmo genético onde é gerada uma população inicial com n indiv́ıduos. Os indiv́ıduos da população, também chamados de cromossomos ou soluções, são definidos por matrizes contendo duas linhas e (tf + 1) colunas. A primeira linha é referente aos valores da variável de decisão do problema de controle u1 e a segunda referente a u2. Cada coluna é referente a um valor do tempo discretizado em t0, t1, t2, t3, . . .,tf , ou seja, do tempo inicial t0 até o tempo final tf . Assim, cada elemento (gene) (i, j), i = 1, 2 e j = 0, 1, . . . , f , desta matriz (cromossomo) representa o valor da variável ui no tempo tj como ilustra a Figura 7 Os cromossomos são ilustrados na Figura 7. Figura 7 - Estrutura do cromossomo definida por matriz de 2 linhas e tf+1 colunas. Esta população evolui de acordo com operadores genéticos de seleção, 29 cruzamento (crossover) e mutação, de forma que promova uma tendência dos in- div́ıduos representarem soluções cada vez melhores à medida que o processo evolutivo continua. 5.3 Avaliação A avaliação tem a finalidade de estabelecer uma medida de qualidade para o indiv́ıduo. Cada indiv́ıduo é avaliado de acordo com sua aptidão ou fitness. No algoritmo proposto, a avaliação dos indiv́ıduos da população é feita usando um valor atribúıdo ao ńıvel de dominância em que está a solução. Visto que o problema (6) é multiobjetivo, pode-se determinar, dentre as soluções da população, as soluções dominadas e não-dominadas da seguinte forma: Sejam u* e u** soluções da população, se J1(u*) < J1(u**) e J2(u*) > J2(u**) então diz-se que u* domina u**. Se u* não for dominada por qualquer outra solução da população, diz-se que u* é uma solução não-dominada ou eficiente. Todas as soluções não-dominadas na população recebem o ńıvel de dominância igual a 1. Para determinar as soluções com ńıvel de dominância igual a 2, retira-se todas as soluções de ńıvel 1 da população e determina quais são as soluções não-dominadas do conjunto restante. Para determinar as de ńıvel 3, faz o mesmo procedimento, retirando as soluções de ńıveis 1 e 2 da população. Desta forma pode-se classificar todas as soluções da população a partir do ńıvel de dominância, conforme mostra a Figura 8. A aptidão ou fitness do indiv́ıduo é inversamente proporcional ao ńıvel de dominância. Bons indiv́ıduos tem ńıvel de dominância iguais ou próximo de 1. Assim é sugerido o cálculo do fitness de cada indiv́ıduo na forma: fi = 1 nivel(i) em que fi é o fitness do indiv́ıduo i, nivel(i) é o ńıvel do indiv́ıduo i. As soluções com fitness igual a 1 comporão a elite em cada geração. Uma solução aproximada para o conjunto de soluções não-dominadas 30 do problema (6) será dada pela elite da última geração. Figura 8 - Nı́vel de dominância dos indiv́ıduos representados no espaço objetivo. 5.4 Elite Um dos operadores de maior importância no AG, a elite possui o obje- tivo de armazenar os melhores cromossomos impedindo que os processos de seleção, crossover, mutação e migração modifiquem totalmente a caracteŕıstica inicial da po- pulação. Este operador é atualizado a cada geração, se necessário, ao encontrar soluções cada vez melhores. No AG proposto, todos os indiv́ıduos que apresentam o fitness igual a 1 farão parte da elite. A elite da última geração será a melhor aproximação da Fronteira de Pareto. 5.5 Seleção A seleção elege os cromossomos, dando maior chance aos indiv́ıduos de melhor fitness para realizarem o próximo passo: crossover. A seleção pode ser feita de várias maneiras, o método utilizado neste problema é apresentado a seguir. 31 Roleta Viciada O método emprega o prinćıpio de probabilidade de sobrevivência do indiv́ıduo mais apto. É usada uma probabilidade de seleção proporcional ao seu fitness, (Reeves & Beasley, 1993). Os indiv́ıduos de melhor fitness possuem melhores chances de serem escolhidos, enquanto os de pior fitness tem suas chances reduzidas. A probabilidade do i-ésimo indiv́ıduo ser selecionado é medida da forma: pi = fi∑n i=1 fi Os indiv́ıduos considerados de melhor aptidão recebem uma maior fatia na roleta de forma a aumentar suas chances de seres escolhidos, enquanto os outros recebem uma fatia menor. Na Figura 9, os indiv́ıduos a, b, c e d estão dispostos na roleta com suas respectivas aptidões. Observe que o indiv́ıduo a possui maiores chances de ser selecionado: Figura 9 - Ilustração do método da Roleta Viciada Depois de associadas as fatias a cada indiv́ıduo, sorteia-se um número r no intervalo [0, 1] e em seguida compara seu valor com a probabilidade acumulada qj, com q0 = 0 e qj = Σj k=1pk, j = 1, 2, . . . , n, em que n é o número total de indiv́ıduos. 32 Assim, se qj−1 ≤ r ≤ qj, o indiv́ıduo selecionado deve ser j. Para ilustrar, sejam sorteados dois valores de r, r = 0, 127 e r = 0, 862. Ao compararmos seus valores na faixa da probabilidade acumulada, Figura 10, os indiv́ıduos escolhidos seriam a e c. Figura 10 - Probabilidade acumulada usando dados do gráfico da Figura 9. No algoritmo proposto utilizou-se o método da roleta para fazer a seleção dos indiv́ıduos que participarão do crossover. 5.6 Crossover O crossover, operador primário da reprodução, acontece quando dois indiv́ıduos da população intermediária são selecionados (pais) e pareados para formar novos indiv́ıduos (filhos) a partir da troca dos genes dos indiv́ıduos selecionados, (Goldbarg & Luna, 2005). Não há garantias de que o crossover gere soluções melhores que os in- div́ıduos pais, no entanto, as soluções são aquelas que sobreviveram a algum processo de seleção e que possivelmente transmitirão boas caracteŕısticas à nova geração. De acordo com a literatura, o crossover é aplicado segundo uma certa probabilidade, chamada taxa de crossover, que varia entre 60% a 90%. Os três tipos de crossover que mais se destacam são: 1. Crossover de 1 ponto: é o tipo mais comum e o mais presente em algoritmos. Dois pais são pareados, um ponto de corte é escolhido aleatoriamente e então 33 dois filhos são obtidos a partir da troca de informações. Na Figura 11, os pais P1 e P2 geram dois filhos F1 e F2 a partir de um ponto de corte aleatório. Figura 11 - Crossover de 1 ponto de posição aleatória. 2. Crossover de n-pontos: extensão do crossover de 1 ponto. Neste caso, são sorteados n pontos de corte para a troca do material genético. Exemplo para n = 2 encontra-se na Figura 12. Figura 12 - Crossover de 2 pontos de posição aleatória. 3. Crossover uniforme: Para cada par de pais é gerada uma sequência de números aleatórios 0− 1. Se os números forem 1’s então copia-se a informação gênica do PAI-1 e se forem zeros copia-se do PAI-2. O outro filho é obtido inversamente, Figura 13. 34 Figura 13 - Crossover uniforme. Para a implementação do algoritmo genético proposto neste trabalho foi usado o método de crossover de 1 ponto. A ilustração de como foi realizado o operador está na Figura 14. Figura 14 - Ilustração de como o operador Crossover age sobre os indiv́ıduos. 35 5.7 Mutação Finalizado o processo do crossover, a população é submetida a outro operador, a mutação, que modifica um ou mais genes do cromossomo. De acordo com Goldberg (1953) apud Aliano Filho (2012), a importância deste processo está em garantir alternativas de exploração, mantendo um ńıvel mı́nimo de abrangência na busca e impedir a convergência prematura, além de manter a diversidade da população. O procedimento da mutação é o estabelecimento de uma probabilidade de ocorrência, de preferência uma baixa probabilidade (pm ≤ 5%), para que não se perca as boas informações genéticas da população. Neste trabalho, a mutação foi feita seguindo os passos: 1. Define-se aleatoriamente, através da probabilidade de mutação (pm), quais indiv́ıduos da atual geração sofrerão mutação; 2. Sorteia-se se a mutação em cada indiv́ıduo será realizada em u1 ou u2 ou em ambos (u1 e u2); 3. Sorteia-se quantas e quais colunas serão trocadas por um novo valor; 4. Sorteia-se o valor que cada gene vai receber, ou seja, o valor a ser inserido nas linhas e colunas sorteadas na matriz (indiv́ıduo). Na Figura 15, os segundos e quintos genes foram sorteados para sofrer mutação. 36 Figura 15 - Exemplo de um processo de mutação sofrida em dois genes de cada linha da matriz. No final deste processo, é obtida uma nova população. Esta população passará pelo processo da avaliação e terá os indiv́ıduos da elite atualizados. 5.8 Migração O operador migração é executado quando existir a necessidade de re- novação de uma subpopulação. A fim de evitar a convergência à um ótimo local, a migração é utilizada com o objetivo de diversificar a população, assim como a mutação. No AG proposto, a migração é realizada substituindo os indiv́ıduos repetidos ou com valores dos objetivos muito próximos, por novos indiv́ıduos cria- dos aleatoriamente. Além de evitar repetições, esta prática oferece diversidade aos indiv́ıduos. 5.9 Atualização, Finalização e Critério de parada Segundo Goldberg (1953), os indiv́ıduos que passaram pelos operadores genéticos crossover, mutação e migração serão inseridos na população inicial e os melhores cromossomos serão preservados (os piores serão exterminados) de acordo com o tamanho da população criada inicialmente, ou seja, a população retorna ao seu tamanho original. Essa nova população corresponderá à população inicial da próxima geração. Este processo continua até que um critério de parada pré-estabelecido seja 37 satisfeito, que pode ser um número fixo de gerações ou uma finalização para casos onde não há melhoria significativa na população. Neste trabalho, o critério de parada do AG corresponde ao final do processo das G gerações estabelecidas para o algoritmo. Um pseudocódigo do AG é ilustrado na Tabela 2. Tabela 2. Pseudocódigo do AG multiobjetivo. Algoritmo AG 1 g ← 0; 2 Gere a população inicial popINICIAL; 3 Avalie popINICIAL; 4 Separe a elite = u*; 5 enquanto (os critérios de parada não estiverem satisfeitos) faça 6 g ← g + 1; 7 Aplique a Seleção; 8 Aplique o Crossover ; 9 Aplique a Mutação; 10 Crie a nova população pop(g); 11 Avalie os indiv́ıduos da pop(g); 11 Atualize a elite; 12 Aplique a Migração; 13 Analise o critério de parada; 14 fim enquanto 15 Apresente a elite como solução do problema 16 fim AG Uma grande vantagem da utilização do algoritmo genético para oti- mização do problema (6) é a facilidade de inserção de restrições, como por exemplo 38 a inserção de condição final sobre as variáveis de estado do sistema (3). As res- trições, em geral, dificultam muito a utilização de métodos tradicionais de resolução de equações diferenciais. Neste trabalho é proposto a seguinte restrição para o número de fêmeas fertilizadas: F (t) ≤ Ffixo, para todo t ≥ tfixo em que Ffixo pode ser um valor de referência para F (t) e tfixo um valor qualquer desejado para t. 6 EXPERIMENTOS COMPUTACIONAIS O algoritmo genético multiobjetivo foi implementado utilizando o soft- ware MATLAB 7.14.0.739 (R2012a) em um computador Intelr CoreTM2 Duo 2GB de memória RAM do Laboratório Cient́ıfico de Informática (LCI) do Departamento de Bioestat́ıstica do Instituto de Biociências de Botucatu. Os parâmetros a seguir, Tabela 3, retirados de Thomé (2007), foram utilizados para implementar o sistema de otimalidade (6). Neste trabalho as integrais no problema (6) foram resolvidas utilizando o Método 1/3 de Simpson e o sistema (4) com as condições iniciais (5) foi resolvido com o uso do Método de Runge-Kutta de ordem 4, estes métodos estão descritos no Apêndice A.1.2 deste trabalho. Tabela 3. Parâmetros utilizados no sistema de otimalidade (6). ϕ C γ r β βS µA µI µF µM µS tfixo Ffixo 0, 50 13 0, 07 0, 50 1 0, 70 0, 05 0, 05 0, 05 0, 10 0, 10 60 0, 10 Utilizando os dados da Tabela 3 nas equações (5), foram determinadas as condições iniciais, que estão apresentadas a seguir:  A(0) = 8, 3200 I(0) = 0, 2773 F (0) = 5, 5467 M(0) = 2, 9120 S(0) = 0 40 Na Tabela 4 encontram-se os parâmetros utilizados na implementação do algoritmo genético. Tabela 4. Parâmetros utilizados no AG. n G T P pm pc 500 250 120 10−16 0, 05 0, 80 em que n é o número de indiv́ıduos da população; G é o número de gerações; T é o peŕıodo total de aplicação do controle; P é a penalização utilizada nos indiv́ıduos infact́ıveis; pm é a probabilidade de um indiv́ıduo sofrer mutação, e pc é a porcentagem de indiv́ıduos selecionados para o crossover A população inicial foi constrúıda aleatoriamente com u1(t) variando no intervalo [0; 1] e u2(t) variando no intervalo [0; 10∗M(0)], com t = 0, 1, 2, . . . , 120. Esta variação escolhida para u2(t) baseia-se na prática da técnica SIT que utiliza uma faixa de inserção de mosquitos machos estéreis variando de 5 a 10 vezes o número de machos naturais na população, para ter uma grande vantagem dos machos estéreis sobre os machos naturais. Foi escolhido M(0) pois no tempo inicial ocorre o maior número de machos naturais na população. A variação de u1(t) como definida, permite infactibilidades no sistema, mas por outro lado esta variação faz com que haja maior investigação no espaço de busca de soluções. Assim, o AG proposto utiliza tal intevalo para u1(t) e insere penalidades sobre as soluções infact́ıveis. O problema multiobjetivo (6) foi resolvido utilizando o AG proposto com os dados anteriormente apresentados. As Figuras 16 e 17 apresentam o conjunto de todos os valores otimizados de u1(t) e u2(t), respectivamente. 41 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Tempo (dias) u1 Figura 16 - Conjunto de valores de u1 otimizados em que cada linha é representada por uma solução do problema. 42 0 20 40 60 80 100 120 0 5 10 15 20 25 30 Tempo (dias) u2 Figura 17 - Conjunto de valores de u2 otimizados em que cada linha é representada por uma solução do problema. Observa-se na Figura 16 que o AG sugere uma maior aplicação do controle u1 nos primeiros 20 dias e, em seguida, reduz esta quantidade para ńıveis muito baixos. Na Figura 17, observa-se que para a inserção de mosquito machos estéreis foi sugerida uma grande quantidade do controle na faixa de 20 a 100 dias. Com todos os valores otimizados de u1(t) e u2(t) foi constrúıda a Curva de Pareto e esta está apresentada na Figura 18. 43 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 4 0 0.5 1 1.5 2 2.5 3 3.5 x 10 6 J1 J2 Figura 18 - Curva de Pareto obtida pelo Algoritmo Genético Multiobjetivo. A Curva de Pareto é constrúıda utilizando o conjunto das soluções não-dominadas obtidas na otimização do problema multiobjetivo (6) pelo algoritmo genético proposto neste trabalho. Esta curva, obtida pelo processo multiobjetivo, fornece uma variedade de soluções que permitem explorar as mais diversas opções de investimentos otimizados. A t́ıtulo de ilustração, são apresentados os valores encontrados para as variáveis de controle u1(t) e u2(t) de uma solução eficiente espećıfica, representada por “*” na Curva de Pareto da Figura 18. 44 0 20 40 60 80 100 120 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time (days) V ar iá ve l d e co nt ro le u 1 Figura 19 - Variável de controle u1. 0 20 40 60 80 100 120 0 5 10 15 20 25 Time (days) V ar iá ve l d e co nt ro le u 2 Figura 20 - Variável de controle u2. 45 Pode ser observado nas Figuras 19 e 20 que a estratégia otimizada sugerida pelo algoritmo genético com este investimento espećıfico em J1 e J2 foi utilizar grande quantidade de inseticida no ińıcio do processo de controle e posterior- mente, baixar esta quantidade de inseticida a um ńıvel muito baixo. Com relação ao controle genético, a estratégia foi aumentar gradativamente a inserção de mosquitos machos estéreis no vigésimo dia e manter em ńıvel alto em quase todo processo. Esta estratégia otimizada para u1(t) e u2(t) permitiu um decrescimento em todos os seg- mentos da população de mosquitos naturais, Figuras 21-24, e permitiu a preservação dos machos estéreis no ambiente, Figura 25. A seguir são apresentadas as Figuras 21-25 que exemplificam a evolução respectivamente das populações de mosquitos na fase aquática, de fêmeas imaturas, fêmeas fertilizadas, machos naturais e machos estéreis, para um ponto representativo na Curva de Pareto, representado por “*” na Figura 18. 0 20 40 60 80 100 120 0 1 2 3 4 5 6 7 8 9 Tempo (dias) E vo lu çã o da p op ul aç ão n a fa se a qu át ic a A Figura 21 - Evolução da população na Fase Aquática A. 46 0 20 40 60 80 100 120 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Tempo (dias) E vo lu çã o da p op ul aç ão d e fê m ea s im at ur as I Figura 22 - Evolução da população de Fêmeas Imaturas I. 0 20 40 60 80 100 120 0 0.5 1 1.5 2 2.5 3 Tempo (dias) E vo lu çã o da p op ul aç ão d e fê m ea s fe rt ili za da s F Figura 23 - Evolução da população de Fêmeas Fertilizadas F. 47 0 20 40 60 80 100 120 0 1 2 3 4 5 6 Tempo (dias) E vo lu çã o da p op ul aç ão d e m ac ho s na tu ra is M Figura 24 - Evolução da população de Machos Naturais M. 0 20 40 60 80 100 120 0 10 20 30 40 50 60 70 80 90 100 Tempo (dias) E vo lu çã o da p op ul aç ão d e m ac ho s es té re is S Figura 25 - Evolução da população de Machos Estéreis S. 48 Os resultados obtidos pelo AG multiobjetivo apresentam soluções efi- cientes, as quais são calculadas considerando diversos ńıveis de investimento nos objetivos propostos, variando de pequenos a grandes investimentos, dando ao to- mador de decisões liberdade para optar pelo investimento mais apropriado ao seu interesse, seja ele econômico, social ou ambiental. 7 CONCLUSÃO Este trabalho abordou o estudo dos processos biológicos e matemáticos de um modelo matemático que descreve a dinâmica populacional do mosquito trans- missor da dengue e do modelo multiobjetivo para auxiliar na otimização da aplicação dos controles qúımico (utilizando inseticida) e genético (com liberação de machos estéreis no ambiente natural) deste mosquito, visando a diminuição dos custos com inseticida e com a produção de machos estéreis, minimizando a quantidade fêmeas fertilizadas e preservando os machos estéreis. Foram propostas técnicas numéricas utilizando algoritmos genéticos para resolução do modelo de controle ótimo multi- objetivo aplicado a problemas de combate ao mosquito da dengue. Por fim, foram discutidos os resultados computacionais obtidos com a aplicação do algoritmo pro- posto. Na abordagem mono-objetivo (Thomé (2007)) a análise do comporta- mento dos controles sobre a função objetivo e as variáveis de estado do sistema é feita considerando diferentes situações de acordo com a necessidade de cada peso na função objetivo, necessitando assim de um número elevado de simulações. Na abordagem multiobjetivo, as soluções eficientes são calculadas considerando diversos ńıveis de investimento, variando sua importância no objetivo, de pequenos a grandes investimentos, dando ao tomador de decisões liberdade para optar pelo investimento mais apropriado ao seu interesse, seja nos custos econômicos ou nos impactos sociais e ambientais de cada controle. Neste trabalho, para ilustração da técnica multiobjetivo foi implemen- tado o problema de otimização multiobjetivo com 2 funcionais, mas o AG proposto pode ser facilmente implementado com mais de 3 funcionais, de acordo com o inte- 50 resse do usuário, sem nenhuma mudança em sua estrutura. O AG proposto apresenta grande flexibilidade na inserção de restrições no problema de controle ótimo, o que em geral dificulta a resolução dos sistemas de equações diferenciais. Isto permitiu um grande ganho em termos de controle das variáveis de estado devido a inserção de condições finais para F. Outra vantagem é a penalização sobre soluções infact́ıveis que permitiu maior exploração da região de factibilidade. O AG é de fácil implementação e mostrou-se uma boa ferramenta para determinação de controles otimizados para o problema em questão. REFERÊNCIAS BIBLIOGRÁFICAS ALIANO FILHO, A. Metaheuŕısticas em um problema de rotação de culturas. Botucatu-SP, 2012. 136p. Dissertação (Mestrado) - UNESP. ARENALES, S.; DAREZZO, A. Cálculo Numérico: Aprendizagem com Apoio de Software. 1. ed. São Paulo: CENGAGE Learning, 2010. 364p. BARSANTE, L. S.; CARDOSO, R. T. N.; ACEBAL, J. L. Otimização do Controle de Gastos com Inseticidas e Machos Estéreis no Combate da Dengue através do Algoritmo Genético. SBPO, 2011. CASA DAS CIÊNCIAS. Aedes Aegypti. Dispońıvel em . Acesso em 08 de Agosto de 2013., 2013. COELLO, C. A. C. Fundamentals of Evolutionary Multi-Objective Optimization. Evolutionary Computation Group, p.1–23, 2010. CUNHA, M. C. C. Métodos Numéricos. 2. ed. Campinas: Editora da Unicamp, 2000. DEB, K.; KALYANMOY, D. Multi-Objective Optimization Using Evolutio- nary Algorithms. New York, NY, USA: John Wiley & Sons, Inc., 2001. DENGUE. Extensão geográfica da dengue. http : //www.dengue.org.br/dengue mapas.html. Acesso em 19 de Agosto de 2013, 2013. 52 DONALÍSIO, M. R.; GLASSER, C. M. Vigilância entomológica e controle de vetores do dengue. Revista Brasileira de Epidemiologia, v.5, n.3, 2002. EMBRAPA. Embrapa entra na guerra contra o mosquito da dengue. . Acesso em 19 de Setembro de 2013., 2010. ESTEVA, L.; YANG, H. M. Control of Dengue Vector by the Sterile Insect Technique Considering Logistic Recruitment. TEMA Tend. Mat. Apl. Comp., v.7, n.2, p.259–268, 2006. FRANCO, N. B. Cálculo Numérico. Campinas: Pearson Prentice Hall, 2006. FREITAS, F. G.; MAIA, C. L. B.; COUTINHO, D. P.; CAMPOS, G. A. L.; SOUZA, J. Aplicação de Metaheuŕısticas em Problemas de Engenharia de Software: Revisão de Literatura. II Congresso Tecnológico Infobrasil, 2009. FÍSICA COMPUTACIONAL. Integral. Dispońıvel em . Acesso em 19 de Setembro de 2013., 2013. GOLDBARG, M. C.; LUNA, H. P. L. Otimização combinatória e programação linear: modelos e algoritmos. 4. ed. Rio de Janeiro: Elsevier Editora Ltda., 2005. GOLDBERG, D. E. Genetic algorithms in search, optimization, and machine learning. The University of Alabama: Addison-Wesley, 1953. GUY, B.; SAVILLE, M.; LANG, J.; SIQUEIRA, J. B.; BRICKS, L. F. Desenvol- vimento de uma vacina tetravalente contra dengue. Revista Pan-Amaz Saude [online], v.2, n.2, p.51–64, 2011. INSTITUTO OSWALDO CRUZ. Dengue v́ırus e vetor. Dispońıvel em . Acesso em 09 de Agosto de 2013., 2013. 53 KEELING, M. J.; ROHANI, P. Modeling Infectious Diseases in Humans and Animals. New Jersey: Princeton University Press, 2008. 385p. MARINHO, L. A. C. Sintomas da Dengue. Dispońıvel em . Acesso em 09 de Agosto de 2013., 2013. MINISTÉRIO DA SAÚDE. Mosquito transgênico combaterá a dengue. . Acesso em 08 de Agosto de 2013., 2012. MINISTÉRIO DA SAÚDE. Diretrizes Nacionais para a Prevenção e Controle de Epidemias de Dengue. . Acesso em 19 de Setembro de 2013, 2013. ORGANIZAÇÃO MUNDIAL DA SAÚDE. DENGUE Guidelines for diagnosis, tre- atment, prevention and control., 2009. ORGANIZAÇÃO MUNDIAL DA SAÚDE. Dengue and dengue haemorrhagic fever. Dispońıvel em . Acesso em 09 de Agosto de 2013, 2013. REEVES, R. C.; BEASLEY, J. E. Modern Heuristic Techniques for Combinatorial Problems. London: Blackwell Scientific Publications, 1993. TEIXEIRA, M. G.; BARRETO, M. L.; GUERRA, Z. Epidemiologia e medidas de prevenção do Dengue. Inf. Epidemiol. Sus [online], v.8, n.4, p.5–33, 1999. TEIXEIRA, M. G.; COSTA, M. C.; COELHO, G.; BARRETO, M. L. Recent shift in age pattern of dengue hemorrhagic fever, Brazil. Emerging Infectious Diseases, v.14, n.10, 2008. 54 THOMÉ, R. C. A. Controle Ótimo Aplicado na Estratégia de Combate ao Aedes Aegypti Utilizando Inseticida e Mosquitos Estéreis. Campinas, 2007. 136p. Tese de doutorado - Universidade Estadual de Campinas. 8 APÊNDICE A MÉTODOS NUMÉRICOS Aqui serão apresentadas técnicas e conceitos matemáticos e computa- cionais necessários na resolução do trabalho. A.1 Método Numérico de Runge-Kutta O Runge-Kutta é o método mais utilizado dentre os apropriados para calcular a solução aproximada de problemas de valor inicial (PVI) devido sua sim- plicidade, precisão e versatilidade nas aplicações, (Cunha, 2000). Este método pode ser entendido como um aperfeiçoamento do método de Euler, com uma estimativa melhorada da derivada da função. Desenvolvido em 1901 por dois matemáticos alemães Carl David Tolmé Runge (1856-1927) e Wilhelm Kutta (1867-1944), o método apresenta precisão equi- valente ao método de Taylor, com a vantagem de evitar o cálculo de derivadas de ordem elevada evitando a complexidade anaĺıtica e um significativo esforço compu- tacional. De acordo com (Arenales & Darezzo, 2010), o método de Runge-Kutta é baseado na avaliação da função f(x, y) em alguns pontos. Definição Seja o PVI: 56 y′(x) = f(x, y) y(x0) = y0 O método geral de Runge-Kutta de R-estágios é definido por: yn+1 = yn + hϕR(xn, yn, h), em que ϕR(xn, yn, h) = c1k1 + c2k2 + . . .+ cRkR, c1 + c2 + . . .+ cR = 1, com k1 = f(xn, yn) k2 = f(xn + ha2, yn + h(b21k1)) k3 = f(xn + ha3, yn + h(b31k1 + b32k2)) ... kR = f(xn + haR, yn + h(bR1k1 + . . .+ bR,R−1kR−1)), sendo a2 = b21 a3 = b31 + b32 a4 = b41 + b42 + b43 ... aR = bR1 + bR2 + . . .+ bR,R−1. A aproximação de yn+1 é calculada a partir de yn e uma “média” de valores da função f(x, y) em vários pontos e os parâmetros cR, aR e brs podem ser escolhidos de modo que o método tenha a mesma ordem de um método de Taylor, definindo a ordem do método de Runge-Kutta. 57 A.1.1 Runge-Kutta de ordem 4 Dentre os métodos de Runge-Kutta, o de quatro estágios é o mais utilizado na solução numérica de problemas com equações diferenciais ordinárias, por ser uma combinação de simplicidade, alta precisão e economia. Na definição de Franco (2006), o método é dado por: yn+1 = yn + h 6 (k1 + 2(k2 + k3) + k4) onde k1 = f(xn, yn), k2 = f(xn + 1 2 h, yn + 1 2 hk1), k3 = f(xn + 1 2 h, yn + h(b31k1 + b32k2)), k4 = f(xn + h, yn + h(b41k1 + b42k2 + b43k3)). (7) conhecido como método de Runge-Kutta de ordem 4. Os parâmetros c1, c2, c3, c4, a1, a2, b21, b31, b41, b42 e b43 são determi- nados a partir de k1, k2, k3 e k4 por Taylor, em torno do ponto (xn, yn) até a ordem 4, expressando o método da seguinte forma: yn+1 = yn + (. . .)h+ (. . .)h2 + (. . .)h3 + (. . .)h4 +O(h5) e, então, iguala-se os coeficientes de h, h2, h3 e h 4 com o método de Taylor de ordem 4. Abaixo temos sua representação algoŕıtmica, Tabela A.1.1. 58 Tabela 5. Pseudocódigo do método RK de ordem 4. Dados y0, h e f(x, y) 1 Para k = 0, 1, 2, . . . , faça 2 m0 = hf(xk, yk) 3 m1 = hf ( xk + h 2 , yk + m0 2 ) 4 m2 = hf ( xk + h 2 , yk + m1 2 ) 5 m3 = hf(xk+1, yk +m1) 6 yk+1 = yk + 1 6 (m0 + 2m1 + 2m2 +m3) A.1.2 Método de Runge-Kutta de ordem 4 para Sistemas de Equações Diferenciais Para resolver um Sistema de Equações Diferenciais, Franco (2006) con- sidera um sistema de n equações diferenciais de primeira ordem:  y′1 = f1(x, y1, y2, . . . , yn) y′2 = f2(x, y1, y2, . . . , yn) ... y′n = fn(x, y1, y2, . . . , yn) (8) podendo ser escrito da forma: y′ = f(x, y) onde y, y′ e f são vetores cujas componentes são yi, y ′ i e fi, com i = 1, 2, . . . , n. Da mesma forma que nas equações diferenciais foi necessário impor uma condição inicial, para que o sistema apresente uma única solução, temos então: y(x0) = y0 59 Para resolver o sistema (8) pelo método de Runge-Kutta de ordem 4, devemos calcular a cada passo: yn+1 = yn + hØR(xn, y n 1 , y n 2 , . . . , y n m, h) em que yn+1 =  y1(n+1) y2(n+1) ... ym(n+1) , yn =  y1(n) y2(n) ... ym(n) , ØR =  ØR1 ØR2 ... ØRm , h ∈ ℜ sendo ØR,i(xn, y n i , ..., y n m) = c1k1(i) + c2k2(i) + c3k3(i) + c4k4(i) i = 1, . . . ,m com c1 + c2 + c3 + c4 = 1 e k1(i) = fi(xn, y n 1 , y n 2 , . . . , y n m) k2(i) = fi ( xn + 1 2 h, yn1 + 1 2 hk1(i), y n 2 + 1 2 hk1(i), . . . , y n m + 1 2 hk1(i) ) k3(i) = fi ( xn + 1 2 h, yn1 + 1 2 hk2(i), y n 2 + 1 2 hk2(i), . . . , y n m + 1 2 hk2(i) ) k4(i) = fi(xn + h, yn1 + hk3(i), y n 2 + hk3(i), . . . , y n m + hk3(i)) i = 1, . . . ,m A.2 Integração Numérica Outro conceito matemático relevante para este estudo é a integração numérica. Segundo Franco (2006), integrar numericamente uma função y = f(x) num dado intervalo [a, b] consiste, em geral, em integrar um polinômio Pn(x) que aproxime f(x) no dado intervalo. De maneira geral, temos: I = ∫ b a f(x)dx, 60 onde f(x) é cont́ınua com derivadas cont́ınuas em [a, b]. Quando posśıvel, determina- se a primitiva F (x) de f(x), tal que F ′(x) = f(x) e o valor da integral será: I = ∫ b a f(x)dx = F (b)− F (a). Graficamente, pode-se relacionar a integral I = ∫ b a f(x)dx com a área S representada na Figura 3, considerando a função f(x) ≥ 0, ∀x ∈ [a, b], entre a curva e o eixo das abcissas. Fonte: F́ısica Computacional (2013) Figura 26 - Área S relacionada à integração I A.2.1 Regra 1/3 de Simpson Em Franco (2006), para estabelecer a fórmula de Simpson, interpola-se f(x) usando um polinômio de 2o grau (polinômio de Lagrange) que coincide com essa função nos pontos x0, x1 e x2. Em seguida, toma-se n = 2, x0 = a, x1 = (a+ b)/2 e x2 = b no polinômio considerado, obtendo:∫ x2 xo f(x) = Σ2 k=0fkhC 2 k (9) de forma que, integrando os polinômios de 2o grau, estabelece-se os pesos da fórmula de Simpson: 61 w0 = ∫ x2 x0 (x− x1)(x− x2) (x0 − x1)(x0 − x2) dx = h 3 w1 = ∫ x2 x0 (x− x0)(x− x2) (x1 − x0)(x1 − x2) dx = 4h 3 w2 = ∫ x2 x0 (x− x0)(x− x1) (x2 − x0)(x2 − x1) dx = h 3 O cálculo das integrais pode ser simplificado lembrando que x1−x0 = h, x2−x1 = h e x2−x0 = 2h. Substituindo os valores de C2 k , k = 0, 1, 2 em (9), obtemos: ∫ x2 x0 f(x)dx ≃ h [ 1 3 f(x0) + 4 3 f(x1) + 1 3 f(x2) ] ou seja, ∫ x2 x0 f(x)dx ≃ h 3 [f(x0) + 4f(x1) + f(x2)] fórmula conhecida como Regra de 1/3 de Simpson. A.2.2 Regra 1/3 de Simpson Generalizada A generalização da Regra de 1/3 de Simpson para integração em um intervalo [a, b] é feita dividindo este intervalo em um número par (2n) de subintervalos de tamanho h = b− a 2n em que a = x0 e b = x2n em cada 2 subintervalos consecutivos. Assim: ∫ x2n x0 f(x)dx ≈ h 3 [ f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + 2f(x4) + . . .+ 2f(x2n−2) + 4f(x2n−1) + f(x2n) ] sendo esta a Regra de 1/3 de Simpson Generalizada.