Métodos de Calibração de Modelos hidrológicos Carlos Ruberto Fragoso Júnior Marllus Gustavo Ferreira Passos das Neves Sumário Conceito básicos O que é calibração? Problemas comuns na calibração de modelos hidrológicos Ciclo da calibração Métodos de calibração Função objetivo Técnicas numéricas Busca aleatória Técnicas iterativas; Busca direta; Técnicas de otimização global; Algoritmos genéticos Critérios de parada O que é calibração? Procura de valores dos parâmetros de um modelo matemático que resultem em uma boa concordância entre dados observados e calculados O erro é minimizado! Calibração - Otimização Otimização encontrar o mínimo ou o máximo de uma função 160 140 120 100 80 60 40 20 0 0 10 20 30 Calibração - Otimização Otimização calibração: quando a função a ser minimizada representa a discrepância entre os dados calculados por um modelo e os dados medidos Calibração o modelo “quer ser” o sistema que ele vai simular Problemas comuns em modelos hidrológicos Encontrar um conjunto ótimo de parâmetros que ajusta um evento de cheia ou uma série de vazões ou vazões mínimas Exemplo encontrar o coeficiente do reservatório linear simples que ajusta adequadamente uma recessão de um hidrograma Q(t+Dt) = Q(t) . exp(-Dt/k) Exemplo Exemplo 3:50 Exemplo Q(t+Dt) = Q(t) . exp(-Dt/k) Primeiro teste: k = 20 Problemas na calibração de modelos hidrológicos Modelos hidrológicos geralmente tem muitos parâmetros Não lineares Técnicas de otimização automáticas Usar Funções Objetivo Ciclo da calibração Critérios para mudança dos parâmetros Critérios para um “bom ajuste” (Função objetivo) Rodar o modelo Novos valores de parâmetros Critérios de parada Verificar o erro ou FO Métodos de calibração Métodos de calibração Tentativa e erro (Manual) Técnicas numéricas Ajusta os parâmetros manualmente baseado nos resultados Usa algoritmos numéricos para encontrar um conjunto de parâmetros ótimo Aleatório Assume faixa de probabilidade para cada parâmetro Funções Objetivo (FO) Medida do erro – objetivo é minimizar a FO Diferentes funções objetivo Somatório dos erros: compensação de erros Somatório do módulo dos erros Somatório dos erros ao quadrado Somatório de erros relativos Somatório dos desvios dos inversos da vazão Erro de volume (bias) Coeficiente de eficiência de Nash-Sutcliffe 70 Observada 60 simulada 50 vazões 40 30 20 10 0 0 20 40 60 80 tempo 100 120 140 Funções Objetivo N F1 (QOi QCi ) 2 i 1 N F2 QO i QC i i 1 N 1 1 2 F3 ( ) QO QC i i i 1 N F4 ( i 1 N QOi QCi 2 ) QOi Função quadrática dá um peso maior para as vazões maiores a tendência será ajustar melhor as máximas (série contínua) Função módulo também dá maior peso às vazões maiores e é pouco utilizada (superfície não comportada) Função para mínimos dá um peso maior para as vazões mínimas inverte a prioridade de F1 Função relativa procura retirar o peso excessivo dado aos valores absolutos de F1 2 F5 (logQOi logQCi ) i1 Exemplo N F1 (QOi QCi ) 2 i 1 N F2 QO i QC i i 1 N F3 ( i 1 3:50 1 1 2 ) QOi QCi Funções Objetivo Raiz do Erro Médio Quadrado (RSME) RMSE n i 1 ( X obs,i X mo del ,i ) 2 n Raiz do Erro Médio Quadrado Normalizado (NRSME) RMSE NRMSE X obs,max X obs,min Quanto mais perto de zero melhor RMSE NRMSE X obs Funções Objetivo Coeficiente de correlação de Pearson r n i 1 n i 1 ( xi x) ( yi y ) ( xi x) 2 n i 1 ( yi y ) 2 Coeficiente de eficiência de Nash-Sutcliffe E 1 n i 1 n ( X obs ,i X mo del ) 2 i 1 ( X obs ,i X obs ) 2 Quanto mais perto de 1 melhor Funções Objetivo - comparação Caso Ideal E=1 RMSE = 0 Funções Objetivo - comparação Caso 1 Máximas próximas, mínimas dispersas, tenta acompanhar o formato 3,20 3,00 2,80 2,60 2,40 2,20 2,00 1,80 1,60 1,40 1,20 1,00 0,80 0,60 0,40 0,20 0,00 ymodel yobs Média obs 3,50 R² = 0,8668 0 0,53,00 1 1,5 2 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 7,5 8 8,5 9 9,5 1010,51111,5 2,50 E = 0,791 RMSE = 0,318 2,00 1,50 1,00 0,50 0,00 0,00 0,50 1,00 1,50 2,00 2,50 3,00 3,50 Funções Objetivo - comparação Caso 2 Mínimas próximas, máximas dispersas, tenta acompanhar o formato 3,20 3,00 2,80 2,60 2,40 2,20 2,00 1,80 1,60 1,40 1,20 1,00 0,80 0,60 0,40 0,20 0,00 ymodel yobs Média obs 3,50 0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 7,5 8 8,5 9 9,5 1010,51111,5 3,00 R² = 0,8633 2,50 E = 0,843 RMSE = 0,276 2,00 1,50 1,00 0,50 0,00 0,00 0,50 1,00 1,50 2,00 2,50 3,00 3,50 Funções Objetivo - comparação Caso 3 Máximas próximas, mínimas dispersas, não acompanha o formato. Observar os pontos em torno da média E = 0,263 RMSE = 0,598 3,20 3,00 2,80 2,60 2,40 2,20 2,00 1,80 1,60 1,40 1,20 1,00 0,80 0,60 0,40 0,20 0,00 ymodel yobs Média obs 3,50 3,00 R² = 0,4031 0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 7,5 8 8,5 9 9,5 1010,51111,5 2,50 2,00 1,50 1,00 0,50 0,00 0,00 0,50 1,00 1,50 2,00 2,50 3,00 3,50 Funções Objetivo - comparação Caso 4 Mínimas próximas, máximas dispersas, não acompanha o formato. Observar os pontos em torno da média E = 0,541 RMSE = 0,471 3,20 3,00 2,80 2,60 2,40 2,20 2,00 1,80 1,60 1,40 1,20 1,00 0,80 0,60 0,40 0,20 0,00 ymodel yobs Média obs 3,50 0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 7,5 8 8,5 9 9,5 1010,51111,5 3,00 2,50 2,00 1,50 1,00 R² = 0,6427 0,50 0,00 0,00 0,50 1,00 1,50 2,00 2,50 3,00 3,50 Funções Objetivo - comparação Caso 5 Longe dos máximos e dos mínimos 3,20 3,00 2,80 2,60 2,40 2,20 2,00 1,80 1,60 1,40 1,20 1,00 0,80 0,60 0,40 0,20 0,00 ymodel yobs Média obs 3,50 3,00 0 0,5 1 1,5 2 2,5 3 3,5 4 4,5 5 5,5 6 6,5 7 7,5 8 8,5 9 9,5 1010,51111,5 2,50 E = 0,573 RMSE = 0,455 2,00 1,50 1,00 R² = 0,5997 0,50 0,00 0,00 0,50 1,00 1,50 2,00 2,50 3,00 3,50 Técnicas de otimização Cálculo analítico Técnicas numéricas Busca aleatória Busca direta Algoritmos genéticos Cálculo analítico Encontrar pontos da função em que a derivada é zero pode ser rápido, mais elegante, mas depende de alguns requisitos via de regra incomuns nas F.O. de calibração de um modelo chuva-vazão (funções de picos múltiplos, funções descontínuas, ausência da forma analítica da função) 160 140 F(x) a b.x c.x 2 120 100 80 60 40 dF 0 dx 20 0 0 10 20 30 Cálculo analítico 2 variáveis: F.O. analítica F: D R, D C R2 fechado e limitado, f diferenciável e de classe C2 f assume um valor máximo e um mínimo em D (teorema de Weierstrass) Condição necessária existam pontos críticos F 0 2F Condição suficiente incluir o estudo x2 da função Hessiana num ponto interior H(x, y) 2F (xo,yo) yx (a) (b) (c) (d) Se H(xo,yo) Se H(xo,yo) Se H(xo,yo) Se H(xo,yo) 2F xy 2F y2 > 0 e Fxx(xo,yo) > 0 (xo,yo) é ponto de mínimo local > 0 e Fxx(xo,yo) < 0 (xo,yo) é ponto de máximo local < 0 ponto de sela = 0 nada se afirma Exemplo Determine o mínimo da função F 3x2 3 0 x F 3y2 3 0 y F(x, y) x3 y3 3x 3y 4 Pontos críticos: (1,1), (1,-1), (-1,1), (-1,-1) Os pontos (1,1) e (-1,-1) são candidatos a mínimo e máximo, respectivamente: pelo cálculo das derivadas segundas 2F 2F (-1,-1) - 6 2 (-1,-1) 2 x y 2F 2F (1,1) 6 2 (1,1) 2 x y 6x 0 H(x, y) 36xy As candidaturas se confirmam 0 6y Características das Técnicas Numéricas Técnicas automáticas (Nash e Sutcliffe, 1970): Envolve sucessivas variações dos valores dos parâmetros de acordo com alguma regra ou padrão pré-concebidos de incrementos que levam em conta os resultados dos passos anteriores e em particular se a variação melhorou ou não o ajuste. Características das Técnicas Numéricas Definição do ponto de partida: critério para inicializar o processo de tentativa em geral depende mais do problema em questão do que do método Direção de pesquisa: identifica o vetor no qual serão realizadas as alterações das variáveis Espaçamento de cada tentativa: variação que ocorrerá na direção de pesquisa a cada tentativa Critérios de parada: definição dos critérios para aceitar uma determinada solução como o ótimo de uma função Técnicas numéricas - Busca Aleatória Conceito de probabilidade Mais simples Geram-se N valores dos parâmetros, calculam-se as F.O. e escolhe-se o menor valor de F.O Vantagens: funções descontínuas; picos múltiplos Desvantagens: demorado; não existe garantia de atingir o ponto ótimo global 160 “Ótimo” 140 120 100 80 60 40 20 0 0 10 20 30 Técnicas numéricas - Busca direta 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Estratégia de caminhar “morro acima” Função objetivo: F(x1,x2) 5 4.5 4 3.5 x2 3 2.5 2 1.5 Máximo local 1 Máximo global 0.5 0 0 0.5 1 1.5 2 2.5 x1 3 3.5 4 4.5 5 Início: ponto coordenadas (parâmetros) aleatórias 5 4.5 X2 = valor aleatório entre 0e5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 X1 = valor aleatório entre 0 e 5 Determina direção de busca: exemplo x2 = x2 + 0,3; x1 = x1 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Função objetivo melhorou? Não, então tenta no outro sentido 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 F.O. melhorou? Sim, então continua no mesmo sentido 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 F.O. melhorou? Sim, então continua no mesmo sentido 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 F.O. melhorou? Sim, então continua no mesmo sentido 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 F.O. melhorou? Não, então volta para o ponto anterior... ...e muda a direção de busca. 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 F.O. melhorou? Sim, então continua no mesmo sentido E assim segue até encontrar um ponto em que não existe direção de busca que melhore o valor da F.O. 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Método unidirecional ou univariacional É o método anterior 1. Direção de pesquisa paralela aos eixos; 2. Pesquisa em cada direção: espaçamento constante ou variável 3. Critério de parada desvantagens: (ao lado) Método da rotação das coordenadas (Rosenbrook) 1º ciclo igual ao univariacional 2º ciclo com rotação 2 alternativas para pesquisa em cada direção: método original que alterna a pesquisa de cada direção em cada tentativa Primeiro ciclo direção x1 x12 x10 S1 14; y2 = 88 x13 x11 S1 14 1,5x1 12,5; y3 = 63,25 x14 12,5 1,5x1,5 10,25 ; y4 = 34,56 x15 x16 Primeiro ciclo direção x2 x12 4 1 3; y8 = -12,01 x 22 3 1,5.1,0 1,5; y9 = -33,02 10,25 1,5x 2,25 6,88; y5 = 10,52 x 32 1,5 1,5.1,5 0,75; y10 = -47,67 6,88 1,5x3,375 1,82; y6 = 17,16 x 42 0,75 1,5.2,25 4,13; y11= -31,67 Rosenbrock: Método um pouco mais eficiente 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Direção de busca é a que potencialmente dará maior incremento da FO Limitação da busca direta: Ótimos locais 5 4.5 4 Região que atrai solução para o ótimo local 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Tentativa de contornar problema: Busca direta com inicialização múltipla 5 4.5 Várias tentativas; espera se que o ótimo global seja a melhor solução testada. 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Problema: Ineficiente e ineficaz quando a FO tem muitos ótimos locais Técnicas numéricas – Busca direta Busca direta (Rosenbrock e cia.) vantagens: funções descontínuas; otimização por simulação (funções que não podem ser expressas analiticamente - calibração de modelos) desvantagens: funções com picos múltiplos Técnicas numéricas – Busca direta Técnicas numéricas – Algoritmos genéticos (AGs) Pais mais adaptados têm maior probabilidade de gerar filhos (sobrevivência do mais apto = seleção natural) Darwin Técnicas numéricas – Algoritmos genéticos (AGs) • Conceitos da biologia (população, reprodução e gerações) • Filhos são semelhantes aos pais (não completamente iguais) • Pais mais “adaptados” maior probabilidade de reproduzir Natureza indivíduos mais adaptados (com mais “aptidão”) maior probabilidade de sobreviver até chegar à fase reprodutiva e participar do processo de reprodução. Algoritmo: pontos com maior F.O. (aptidão) têm maior probabilidade de serem escolhidos para participar dos complexos Técnicas numéricas – Algoritmos genéticos (AGs) Início Inicialização da população Cálculo da aptidão Solução encontrada? Nova população Seleção Reprodução Mutação Fim Algoritmos genéticos (AGs) Gera-se uma população inicial (possíveis soluções) População avaliada recebe uma nota ou aptidão (qualidade da solução) Mais aptos selecionados e os menos aptos descartados Os selecionados podem sofrer operações de crossover e mutação geram filhos 1ª geração Nova avaliação processo repetido até se tornar satisfatório Algoritmo genético “puro” 1 - gera população (pontos aleatórios) 160 140 120 100 80 60 40 20 0 0 10 20 30 Algoritmo genético “puro” 2 - escolhe pontos para participar do processo de “reprodução” (aqueles com melhor aptidão tem maior probabilidade de escolha 160 140 120 100 80 60 40 20 0 0 10 20 30 Aptidão pode ser igualada à F.O., mas via de regra, é função da F.O. (evitar valores negativos ou convergência prematura extremos locais) Algoritmo genético “puro” 3 - Exemplo de reprodução: escolhidos dois pontos PAIS 160 140 120 100 80 60 40 20 0 0 10 Xa=8 20 Xb=19 30 Binário: Xa=01000 Xb=10011 Selecionados com probabilidade “proporcional à aptidão” mais aptidão, mais chance de serem escolhidos Algoritmo genético “puro” Genética: filhos “recebem” cromossomos dos pais 01000 10011 01011 10000 Filhos: É determinado um (ou mais) ponto de “corte” (aleatório) Filho 1: parte dos “cromossomos” do pai e parte da mãe Filho 2: outra parte dos “cromossomos” do pai e parte da mãe Xa=01011 = 11 Xb=10000 = 16 Algoritmo genético “puro” Na operação de crossover, os pais transmitem aos filhos algumas características 160 140 120 pais 100 80 60 40 filhos 20 0 0 10 20 30 Na mutação (a seguir), surgem novas características próprias dos filhos transmiti-las para outra geração Algoritmo genético “puro” Mutação: evento de baixa probabilidade Genética: filhos “recebem” cromossomos dos pais 01000 10011 É determinado um (ou mais) ponto de “corte” (aleatório) 01011 10100 Filho 1: sem mutação Filhos: Filho 2: mutação Xa=01011 = 11 Xb=10100 = 20 Algoritmo genético “puro” Reprodução de todos os pontos escolhidos resulta na nova geração 160 140 120 100 80 60 40 20 0 0 10 20 30 Mutação melhora a diversidade na população, mas pode não preservar a informação original deve haver uma taxa pequena (0,1 a 5%), mas suficiente Algoritmo genético “puro” Depois de algumas gerações 160 140 120 100 80 60 40 20 0 0 10 20 30 Algoritmo genético “puro” Algumas desvantagens do algoritmo genético puro • Números binários • Transformação de variáveis de base decimal para binária -0,05 +180,3 Variável Y 0000000000 1111111111 Usando 10 bits; Resolução = 0,176 decimal Algoritmo genético “puro” Algumas vantagens do algoritmo genético puro Otimização com números inteiros Diâmetros comerciais Algoritmo genético “puro” Representação binária Historicamente importante tradicional Fácil de utilizar e manipular e simples de analisar teoricamente Problema com parâmetros contínuos e o usuário desejando trabalhar com boa precisão numérica muito armazenamento convergência lenta do algoritmo Representação real melhor compreendida Há diversos operadores para n° real (crossover e mutação) Evolução de complexos misturados (Shuffled complex evolution) SCE - UA Usa técnicas de busca aleatória algoritmos genéticos simplex (Nelder e Mead) Proposto por Duan, Gupta e Sorooshian (U. Arizona) Descrito no livro Sistemas Inteligentes da ABRH SCE-UA Espaço de busca Cria-se uma população de s pontos de forma aleatória População dividida em algumas comunidades (complexos ou conjuntos) Evoluem de forma independente o espaço de busca é explorado em todas as direções Após algumas gerações comunidades envolvidas entre si MISTURA origem a outras comunidades SCE-UA Espaço de busca Subcomplexo ou subconjunto extraído aleatoriamente de um complexo representa os PAIS Maior chance (distribuição triangular de probabilidades) aos melhores pais (menor FO) geração de descendentes Eles geram outros complexos Por vezes, descendentes são substituídos MUTAÇÃO SCE-UA Geram-se novos pontos pra recomeçar o processo até chegar a um nível satisfatório SCE-UA Geram-se novos pontos pra recomeçar o processo até chegar a um nível satisfatório SCE-UA Geram-se novos pontos pra recomeçar o processo até chegar a um nível satisfatório 5 4.5 4 3.5 3 Passo 1 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 2 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 4 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 5 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 6 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 7 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 8 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 9 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 10 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 4.5 4 3.5 3 Passo 20 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 1 - Geração aleatória de pontos 5 Cada ponto conj. de parâmetros 4.5 4 3.5 No caso ao lado cada ponto é um par de parâmetros x1, x2 3 2.5 2 1.5 Geração aleatória com uma distribuição uniforme 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 2 – Hierarquização e formação complexos 5 Hierarquização segundo as FOs de cada ponto 4.5 4 Divisão em complexos “comunidades” 3.5 3 2.5 cada um tem pontos bons e pontos ruins 2 1.5 Exemplo: complexos de 4 pontos 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos Cada complexo tem a oportunidade de evoluir de forma Independente as “comunidades” evoluem de forma independente 5 4.5 4 3.5 3 2.5 Nas “comunidades” surgem os “casais” subcomplexos 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos Seleção de “pais” ou “reprodutores” um subcomplexo é retirado aleatoriamente de cada complexo. 5 4.5 4 3.5 3 A probabilidade de um ponto do complexo participar do sub-complexo é proporcional à FO (aptidão para ser “pai”) 2.5 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 Exemplo: sub-complexo de 3 pontos extraído de um complexo de 4 pontos 4.5 4 3.5 3 A probabilidade de um ponto do complexo participar do sub-complexo é proporcional à FO (aptidão para ser “pai”) 2.5 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 Procedimento de reprodução semelhante ao método Simplex de Nelder e Mead 4.5 4 3.5 3 2.5 Define pior ponto do sub-complexo 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 Procedimento de reprodução semelhante ao método Simplex de Nelder e Mead 4.5 4 3.5 3 2.5 Define centróide dos melhores pontos 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 4.5 Passo de reflexão: distância a = distância b 4 3.5 a 3 Verifica valor da FO no novo ponto, se é melhor do que pior ponto, novo ponto é aceito, se não, vai para o passo de contração. 2.5 b 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 4.5 Passo de contração: distância a = distância b 4 a 3.5 b 3 Verifica valor da FO no novo ponto, se é melhor do que pior ponto, novo ponto é aceito, se não, cria ponto aleatório 2.5 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 4.5 Um novo ponto é gerado no espaço definido pelos limites mínimo e máximo de cada um dos parâmetros no complexo 4 3.5 3 2.5 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 3 – Evolução dos complexos sub-complexos 5 Um novo ponto é gerado no espaço definido pelos limites mínimo e máximo de cada um dos parâmetros no complexo 4.5 4 3.5 Após um n° de gerações prédeterminado (β) comunidades se encontram (mistura de complexos) nova divisão é feita 3 2.5 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 4 – Mistura de complexos 5 Mistura de “comunidades” permite o compartilhamento de informações que cada comunidade adquiriu de forma independente 4.5 4 3.5 3 Complexos novamente agrupados novo conjunto de pontos nova reorganização... 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 4 – Teste de convergência Pontos da amostra testados caso satisfaçam um critério de convergência previamente definido algoritmo se encerra. Caso contrário recomeça no passo 2. 5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 3:500 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 SCE-UA Pais mais adaptados têm maior probabilidade de gerar filhos (sobrevivência do mais apto = seleção natural) Probabilidade de escolha 1 0 Valor da FO Posição no ranking 1) Classificar os pontos do complexo em ordem de FO (ranking) 2) Atribuir probabilidade de escolha para participar do subcomplexo segundo a função do desenho Exemplo Sub-Complexo Complexo 5 5 4.5 4.5 4 4 3.5 3.5 3 3 2.5 2.5 2 2 1.5 1.5 1 1 0.5 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Dois pontos do complexo ficaram fora do sub-complexo. Não necessariamente os piores pontos ficam fora Filhos são semelhantes aos pais Genética: filhos “recebem” cromossomos dos pais 5 Algoritmo SCE-UA: No lugar dos “casais” estão os “sub-complexos”, que são “casais” de q pontos 4.5 4 3.5 a 3 2.5 b 2 1.5 1 0.5 0 3:50 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Calibração automática com SCE-UA Cada ponto valores dos parâmetros escolhidos. FO = coeficiente de Nash-Sutcliffe. Para ser avaliada realizar uma simulação completa (por exemplo, 10 anos de dados diários). 5 4.5 4 3.5 3 700 2.5 600 calculada observada 2 Vazão (m3 /s) 500 400 1.5 300 1 200 0.5 100 0 0 0 01/jun/72 01/jul/72 31/jul/72 30/ago/72 29/set/72 29/out/72 28/nov/72 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Exemplo de aplicação do SCE-UA COLLISCHONN, W. ; TUCCI, C. E. M. Calibração automática global do modelo IPH-2. In: Simpósio Brasileiro de Recursos Hídricos, 2001, Aracaju SE. Anais do Simpósio Brasileiro de Recursos Hídricos, 2001. Objetivos: testar a robustez, a eficiência e o uso de 2 FOs diferentes com série contínua de vazões diárias Robustez capacidade de encontrar sempre o mesmo ponto ótimo em uma série de tentativas independentes indicativo forte de que o ótimo global foi encontrado Eficiência capacidade de encontrar o ponto ótimo em um mínimo número de avaliações da função (iterações do modelo hidrológico) Exemplo: teste 1 Como ter certeza de que um método de calibração pode atingir o ótimo global? Testá-lo, usando uma série de vazões geradas pelo próprio modelo adotando valores de parâmetros arbitrariemente O método de calibração tem que atingir estes valores I0 Ib h Ks Kbas Rmax Alf 50,0 1,0 0,8 5,0 100,0 4,0 2,0 Resultados teste 1 10 aplicações sucessivas atingido sempre o ótimo global (conjunto de parâmetros que gerou a série sintética), em menos do que 10.000 avaliações da FO 140 Valor do parâmetro ao longo do processo 120 Valor do parâmetro I0 Literatura mostra testes semelhantes com métodos Rosenbrock e outros, que não conseguem superar este teste 100 I0 = 50 80 60 40 20 0 0 2000 4000 6000 8000 Número de avaliações da função 10000 12000 Teste 2: calibração dados observados Calibração do modelo IPH-2 (10 vezes) Calibração I0 Ib h Ks Kbas Rmax Alf R2 1 36,04 0,46 0,93 7,52 11,11 2,80 19,99 0,85915 2 36,04 0,46 0,93 7,52 11,16 2,80 19,99 0,85915 3 36.03 0.46 0.93 7,52 11,05 2,80 19,99 0,85915 4 35.91 0.46 0,93 7.56 11.95 2.81 19.69 0,85915 5 36.02 0.46 0,93 7.52 11.09 2.80 19.98 0,85915 6 36.04 0.46 0.93 7.52 11.14 2.80 19.99 0,85915 7 36.04 0.46 0,93 7.52 11.12 2.80 19.99 0,85915 8 36.05 0.46 0,93 7.52 11.13 2.80 19.99 0,85915 9 36.03 0.46 0.93 7.52 11.11 2.79 19.99 0,85915 10 36.04 0.46 0.93 7.52 11.16 2.80 19.99 0,85915 Eficaz indicativo forte de que foi atingido o ótimo global (apesar de não se conseguir provar) Otimização multi-objetivo Objetivo otimizar, ao mesmo tempo, várias FOs vários aspectos da qualidade do ajuste entre hidrogramas calculados e observados Praticamente impossível encontrar um conjunto de parâmetros que produza, ao mesmo tempo, valores ótimos de várias FOs diferentes (imperfeições dos modelos e dos dados de entrada) Otimização multi-objetivo Quase sempre, a opção pela otimização de uma FO implicará na não otimização de outra situação normalmente enfrentada na calibração de modelos chuva vazão bom ajuste dos picos de vazão implica em maus resultados durante os períodos de estiagem A solução não é um ponto, mas uma região de Pareto, ou região de soluções não inferiores ou não dominadas Otimização multi-objetivo 12 11 10 9 F1 F2 Região de Pareto 8 7 F(x) 6 5 4 3 2 1 0 -1 Ótimo F2 Ótimo F1 -2 0 1 2 3 4 5 6 x Região de Pareto 7 8 9 10 Otimização multi-objetivo 12 11 10 9 F1 F2 Região de Pareto 8 7 F(x) 6 5 Piorou em F1 4 3 2 1 Mínimo em F2 0 -1 -2 0 1 2 3 4 5 x x 6= 6 7 8 9 10 Otimização multi-objetivo A otimização mono-objetivo pode ser utilizada de forma iterativa para resolver o problema multi-objetivo: uso de coeficientes, ou pesos, atribuídos a cada uma das FOs a cada calibração diversos pontos que definem, aproximadamente, a região de Pareto (RP) do problema de calibração multiobjetivo “hierarquização de Pareto” ou Pareto ranking encontrar, em apenas um procedimento de otimização, vários pontos na RP MOCOM-UA MOCOM-UA (Multiple-Objective Complex Evolution Universidade do Arizona) Variação do algoritmo SCE-UA (Duan et al., 1992), desenvolvido para problemas monoobjetivo Inicia com a definição dos valores limite dos n parâmetros Adota-se o n° de pontos ns (ns conjuntos de parâmetros) tamanho da população São avaliadas nf FOs matriz de resultados [ns, nf] MOCOM-UA Exemplo de matriz de resultados: n = 1 (1 parâmetro, ns = 3 (três pontos) e nf = 2 (2 FOs) 12 N° do ponto 1 2 3 11 10 9 F1 F2 Região de Pareto 8 7 F(x) 6 5 4 3 x F1 F2 3 4 7 2 1 0 -1 -2 0 1 2 3 4 5 x 6 7 8 9 10 1 0 9 9 4 1 MOCOM-UA Exemplo de 2 parâmetros: n = 2 (dois parâmetros, ns = 5 (cinco pontos) e nf = 2 (duas FOs) Ponto Qualquer (X,Y) MOCOM-UA Com os valores das FOs conjuntos de parâmetros avaliados e hierarquizados pelos critérios de dominância e não dominância Passos: 1. Toma-se os ns pontos ou conjuntos e identifica-se aqueles que são dominados e os que são não dominados 2. Aos pontos que são não dominados atribui-se o índice 1 3. Os pontos com o índice 1 são retirados e os pontos restantes são novamente analisados 4. Aos pontos que são não dominados nesta segunda análise, atribui-se o índice 2 MOCOM-UA Com os valores das FOs conjuntos de parâmetros avaliados e hierarquizados pelos critérios de dominância e não dominância Passos: 5. Os pontos com o índice 2 também são retirados e os pontos restantes são analisados 6. Os passos se repetem até que se encontre um grupo de pontos em que não podem ser definidos dominados e não dominados. Estes pontos recebem o índice MAXR, onde MAXR é o número de passos necessários para não existir mais dominância entre os pontos Quanto menor o índice mais perto da região de Pareto MOCOM-UA: dominância e não dominância 12 11 x = 3 vs x = 7 quem domina quem? 10 9 F1 F2 Região de Pareto 8 7 F(x) 6 5 4 3 2 1 0 -1 -2 0 1 2 3 3 x= 4 5 x 6 x 7= 7 8 9 10 MOCOM-UA: dominância e não dominância 12 11 x = 6 vs x = 7 quem domina quem? 10 9 F1 F2 Região de Pareto 8 7 F(x) 6 5 4 3 2 1 0 -1 -2 0 1 2 3 4 5 x x =6 6 x 7= 7 8 9 10 MOCOM-UA: evolução dos pontos Exemplo de 2 parâmetros: geram-se NMAXR complexos (grupos de pontos). NMAXR n° de pontos que recebem o pior índice na hierarquização Cada complexo possui n+1 (3 pontos), e a cada ponto dele é associada uma probabilidade de que esse ponto participe do processo de evolução Cada complexo evolui uma única vez, gerando NMAXR novos pontos que substituem os NMAXR piores pontos anteriores, mantendo a população estável MOCOM-UA: geração de pontos Gerados ns = 5 pontos (conjuntos de pares X,Y) MOCOM-UA: hierarquização Pontos hierarquizados (dominância e não dominância) MAXR = 3 e há somente 1 ponto que recebe este índice NMAXR = 1 B e C não dominados B C MOCOM-UA: criação de complexos Gerado somente 1 complexo (NMAXR = 1) São escolhidos para participar dele os pontos B e C (além do ponto D) MOCOM-UA: criação de complexos Gerado somente 1 complexo (NMAXR = 1) MOCOM-UA: evolução dos complexos Centróide dos melhores pontos Método simplex de Nelder e Mead (semelhante ao algoritmo SCEUA) MOCOM-UA reflexão Ponto R é aceito se estiver na região válida e se for não dominado em relação aos outros 2 pontos MOCOM-UA contração Caso contrário, o ponto de contração K é aceito imediatamente MOCOM-UA Cada complexo evolui uma única vez, gerando NMAXR novos pontos que substituem os NMAX piores pontos anteriores, mantendo a população estável. A seguir, pontos novamente analisados e hierarquizados processo se repete, até que, na etapa de hierarquização, todos os pontos recebam o mesmo índice não há pontos melhores ou piores. Nesta situação, normalmente, deverá ter sido encontrada uma boa amostra de pontos sobre a região de Pareto MOCOM-UA Região de Pareto Aproximação da região de Pareto Exemplo IPH 2 Faixa válida dos parâmetros Parâmetr Unidade o Io mm.Dt-1 Ib mm.Dt-1 H Ks Dt Ksub Dt Rmáx mm Alf - Valor mínimo 10 0,1 0,0 0,01 30,0 0,0 0,01 Valor máximo 300 10 1,0 10,0 40,0 9,0 20,0 2 FO: Erro volume e RMSE Geração 1 8 7 Erro no volume 6 5 4 3 2 1 0 15000 17000 19000 21000 Soma desvios quadrados 23000 25000 Geração 10 8 7 Erro no volume 6 5 4 3 2 1 0 15000 17000 19000 21000 Soma desvios quadrados 23000 25000 Geração 20 8 7 Erro no volume 6 5 4 3 2 1 0 15000 17000 19000 21000 Soma desvios quadrados 23000 25000 Geração 50 8 7 Erro no volume 6 5 4 3 2 1 0 15000 17000 19000 21000 Soma desvios quadrados 23000 25000 Geração 138 8 7 Erro no volume 6 5 4 3 2 1 0 15000 17000 19000 21000 Soma desvios quadrados 23000 25000 Avaliação da incerteza: usar todos os conjuntos e gerar vários hidrogramas Propagação da incerteza: Q90 calculada, por exemplo, vai de 8,9 a 10,5 m3.s-1, sendo que a Q90 observada é de 9,1 m3.s-1 Quando se otimiza, aceita-se que todos os dados são bons, ou seja, a otimização não faz distinção entre dados bons e dados ruins Sugestões de leitura TUCCI, C. E. M.; COLLISCHONN, W. Ajuste Multiobjetivo dos Parâmetros de um Modelo Hidrológico. Revista Brasileira de Recursos Hídrcos, Vol. 8, n. 3, Jul/Set 2003 p. 27-39 BRAVO, J. M.; COLLISCHONN, W; TUCCI, C. E. M. Verificação da Eficiência e Eficácia de um Algoritmo Evolucionário Multi-objetivo na Calibração Automática do Modelo Hidrológico IPH II. Revista Brasileira de Recursos Hídrcos, Vol. 14, n. 3, Jul/Set 2009 p. 37- 50 Yapo, P. O.; Gupta, H. V.; Sorooshian, S. 1998 Multi-objective global optimization for hydrologic models. Journal of Hydrology, Vol. 204 pp. 83-97. Sorooshian, S.; Gupta, V. K. 1995 Model calibration In: Singh, V. J. (editor) Computer models of watershed hydrology. Water Resources Publications, Highlands Ranch. 1130 p. Sugestões de leitura Duan, Q.; Sorooshian, S.; Gupta, V. 1992 Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resources Research Vol. 28 No. 4. pp. 1015-1031. Duan, Q.; Sorooshian, S.; Gupta, V. 1994 Optimal use of the SCE – UA global optimization method for calibrating watershed models. Journal of Hydrology, Vol 158 pp. 265-284. Sugestões de leitura Bonabeau, E.; Dorigo, M.; Theraulaz, G. 2000 Inspiration for optimization from social insect behaviour. Nature Vol. 406 July pp.39-42. Goldberg, D. 1989 Genetic Algorithms in Search, Optimization and Machine Learning Addison-Wesley, 412 pp. Klemes, V. 1986 Operational testing of hydrological simulation models. Hydrological Sciences Journal V. 31 No. 1 pp. 13-24. Nash e Sutcliffe, 1970 (Journal of Hydrology) Particle Swarm Optimization