CÁSSIO AUGUSTO USSI MONTI OTIMIZAÇÃO APLICADA À ENGENHARIA FLORESTAL LAVRAS - MG 2018 CÁSSIO AUGUSTO USSI MONTI OTIMIZAÇÃO APLICADA À ENGENHARIA FLORESTAL Dissertação apresentada à Universidade Federal de Lavras, como parte das exigências do Programa de Pós- Graduação em Engenharia Florestal, área de concentração em Manejo Florestal, para obtenção do título de mestre. LAVRAS - MG 2018 Ficha catalográfica elaborada pelo Sistema de Geração de Ficha Catalográfica da Biblioteca Universitária da UFLA, com dados informados pelo(a) próprio(a) autor(a). Monti, Cássio Augusto Ussi. Otimização Aplicada à Engenharia Florestal / Cássio Augusto Ussi Monti. - 2018. 104 p. : il. Orientador(a): Lucas Rezende Gomide. . Dissertação (mestrado acadêmico) - Universidade Federal de Lavras, 2018. Bibliografia. 1. Regressão Não Linear. 2. Algoritmo Genético. 3. Planejamento do Transporte Florestal. I. Gomide, Lucas Rezende. . II. Título. O conteúdo desta obra é de responsabilidade do(a) autor(a) e de seu orientador(a). CÁSSIO AUGUSTO USSI MONTI OTIMIZAÇÃO APLICADA À ENGENHARIA FLORESTAL Dissertação apresentada à Universidade Federal de Lavras, como parte das exigências do Programa de Pós- Graduação em Engenharia Florestal, área de concentração em Manejo Florestal, para obtenção do título de mestre. APROVADA em 07 de março de 2018. Prof. Dr. Lucas Rezende Gomide - UFLA Prof. Dr. Fausto Weimar Acerbi Júnior – UFLA Prof. Dr. Carlos Alberto Araújo Júnior - UFMG Orientador Prof. Dr. Lucas Rezende Gomide LAVRAS - MG 2018 À minha família. Dedico. AGRADECIMENTOS À Universidade Federal de Lavras, especialmente ao Departamento de Ciências Florestais, pela oportunidade. À CAPES, pela concessão da bolsa de mestrado. Ao professor Lucas Rezende Gomide, pela orientação, paciência e disposição para ajudar. Aos companheiros do Laboratório 49, tornando este período de mestrado mais agradável e divertido. Aos companheiros pós-graduandos do LEMAF/UFLA. À todos os funcionários do DCF/UFLA e LEMAF/UFLA. À minha família, pelos conselhos e apoio total em todas as minhas decisões. Aos grandes amigos que fiz em Lavras, especialmente aos Camisa 10. MUITO OBRIGADO! RESUMO A modelagem de sistemas é aplicada em vários ramos da engenharia, para representar qualquer tipo de situação operacional através de técnicas estocásticas, como regressão linear e não linear, ou determinística, como modelo de roteamento veicular. Ultimamente é comum a aplicação de inteligência computacional como alternativa a estas técnicas, como heurísticas e meta- heurísticas. Com base nestas informações, a gestão planeja a colheita e o transporte florestal, com objetivo de abastecer a fábrica. Tal planejamento utiliza a modelagem determinística de sistemas através de modelos de programação matemática oriundos da disciplina de pesquisa operacional. Dessa forma, a estrutura da dissertação contemplou três capítulos, sendo a primeira parte composta por levantamento bibliográfico do estado da arte sobre o ajuste de modelos estatísticos e teoria de grafos aplicados ao transporte e otimização de problemas. Na segunda parte buscou-se abordar métodos de ajuste de modelos não lineares, propondo o desenvolvimento de um método alternativo aos algoritmos usualmente aplicados. O mesmo prevê a estimação dos parâmetros iniciais de maneira automática utilizando algoritmo genético como opção. Assim, foram testados três métodos clássicos de ajustes, Gauss-Newton, Levenberg-Marquardt e NL2SOL. A relação hipsométrica foi adotada como base para os testes computacionais, sendo testada uma combinação de 12 bases de dados florestais distintas e 12 modelos não lineares, em 100 repetições para garantir uma confiabilidade da informação. Dentre os métodos testados, Levenberg-Marquardt apresentou maior frequência de sucesso dos ajustes, seguido por NL2SOL e Gauss-Newton. O algoritmo genético não produziu diferença significativa na predição dos parâmetros, podendo ser utilizado na parametrização de modelos estatísticos de regressão. Ao ultima parte da dissertação compôs um estudo envolvendo uma modelagem determinística para o problema de roteamento de veícular aplicado ao transporte florestal (PRVTF). O modelo envolvendo variáveis inteiras e contínuas, considerou em sua função multiobjetivo a redução no número de caminhões, distância percorrida, horas- extra, eficiência de carga e movimentação de gruas como processo integrado. Os resultados indicaram ausência de horas-extra e redução de 72,92% no uso dos caminhões, consequentemente promovendo ainda um aumento no número de viagens/veículo. Em apenas 3,17% das viagens à eficiência de carga não foi obtida, sendo um fato previsto na remoção final da madeira nos talhões. Identificaram-se ainda pontos de melhoria na movimentação das gruas entre talhões devido ao ócio operacional. O modelo PRVTF atua como sistema de gestão integrada eficiente e eficaz no transporte florestal. Palavras-chave: Algoritmo genético. Programação linear inteira mista. Sistema integrado. Parametrização não linear. ABSTRACT System modeling is applied in various branches of engineering to represent any type of operational situation through stochastic techniques, such as linear and non-linear regression, or deterministic, as a model of vehicle routing. Lately it is common application of computational intelligence as an alternative to these techniques, such as heuristics and metaheuristics. Based on this information, planning harvest and transport forestry were design, in order to supply the factory. Such planning uses the deterministic modeling of systems through mathematical programming models from the operational research discipline. Thus, the structure of the dissertation contemplated three chapters, the first part being a bibliographical survey of the state of the art on the fitting of statistical models and theory of graphs applied to the transport and optimization of problems. In the second part, we tried to approach methods of fitting nonlinear regression models, proposing the development of an alternative method to the algorithms usually applied. The same predicts the estimation of the initial parameters automatically using genetic algorithm as an option. Thus, three classical fitting methods, Gauss-Newton, Levenberg-Marquardt and NL2SOL, were tested. The hypsometric relation was adopted as the basis for the computational tests, and a combination of 12 distinct forest databases and 12 nonlinear models were tested in 100 replicates to guarantee information reliability. Among the tested methods, Levenberg-Marquardt showed a greater frequency of fitting success, followed by NL2SOL and Gauss-Newton. The genetic algorithm did not produce significant difference in the prediction of the parameters, and can be used in the parameterization of statistical regression models. The last part of the dissertation composed a study involving a deterministic modeling for the vehicle routing problem applied to forest transportation (VRPFT). The model involving integer and continuous variables considered in its multiobjective function the reduction in the number of trucks, distance traveled, overtime, load efficiency and crane movement as an integrated process. The results indicated absence of overtime and a reduction of 72.92% in the use of trucks, consequently promoting an increase in the number of trips / vehicles. In only 3.17% of the trips to the load efficiency was not obtained, being a fact predicted in the final removal of the wood in the plots. It was also identified points of improvement in the movement of cranes between blocks due to operational idle. The VRPFT model acts as an efficient and effective integrated management system applied to forest transportation. Key-words: Genetic Algorithm. Mixed integer linear programming. Integrated system. Non-linear parameterization. SUMÁRIO PRIMEIRA PARTE ....................................................................... 10 1 INTRODUÇÃO GERAL ............................................................... 11 1.1 OBJETIVO ..................................................................................... 13 2 REFERENCIAL TEÓRICO........................................................... 13 2.1 Regressão não linear ....................................................................... 13 2.2 Algoritmo genético ......................................................................... 16 2.3 Problema de roteamento de veículos .............................................. 20 2.3.1 Teoria de grafos .............................................................................. 20 2.3.2 Problema de roteamento de veículos clássico ................................ 22 2.4 Otimização Multiobjetivo ............................................................... 26 3 CONSIDERAÇÕES GERAIS ........................................................ 28 REFERÊNCIAS ............................................................................. 28 SEGUNDA PARTE - ARTIGOS ................................................... 39 ARTIGO 1 - A HYBRID METHOD FOR FITTING NONLINEAR REGRESSION MODELS .............................................................. 40 1 INTRODUCTION .......................................................................... 42 2 MATERIAL AND METHODS ..................................................... 45 2.1 Data base structure ......................................................................... 45 2.2 Fitting approach .............................................................................. 47 2.2.1 Nonlinear regression models ......................................................... 47 2.2.2 Hybrid method............................................................................... 49 2.2.2.1 Stage 1: Genetic dynamic approach .............................................. 50 2.2.2.2 Stage 2: Statistical approach.......................................................... 55 2.3 Experiment analysis ....................................................................... 60 3 RESULTS....................................................................................... 61 3.1 Statistical regression models and database pattern ......................... 61 3.2 Hybrid method................................................................................ 62 3.2.1 Stage 1: Genetic algorithm ............................................................ 62 3.2.2 Stage 2: Statistical fitting methods ................................................ 64 3.3 Genetic algorithm versus Hybrid method ...................................... 66 4 DISCUSSION AND CONCLUSION ............................................ 69 REFERENCES ............................................................................... 73 ARTIGO 2 – OTIMIZAÇÃO MULTIOBJETIVO NO TRANSPORTE DE MADEIRA .................................................... 82 1 INTRODUÇÃO ............................................................................. 84 2 MATERIAL E MÉTODOS ........................................................... 86 2.1 Descrição da rede de transporte de madeira ................................... 86 2.2 Modelo de programação linear inteira mista multiobjetivo ........... 88 2.2.2 Restrições técnicas de operacionalização ....................................... 89 2.2.3 Características computacionais e avaliação dos resultados ............ 93 3 RESULTADOS .............................................................................. 93 4 DISCUSSÃO E CONCLUSÕES ................................................... 98 REFERÊNCIAS ........................................................................... 103 10 PRIMEIRA PARTE 11 1 INTRODUÇÃO GERAL O setor florestal brasileiro, associado às florestas plantadas, apresenta grande diversidade de produtos abastecendo indústrias de celulose e papel, de painéis de madeira reconstituída, siderúrgicas movidas a carvão e biomassa, entre outras. Segundo a Indústria Brasileira de Árvores (IBÁ), em 2014, o PIB do setor cresceu 1,7% movido primordialmente pela expansão da exportação de celulose. O setor de árvores plantadas apresenta maior crescimento do PIB, cerca de 12,6%, do que o setor agropecuário (0,4%), indústria (-1,2%) e setor de serviços (0,7%) juntos para o ano de 2014 e é 17 vezes maior que o crescimento do PIB brasileiro no mesmo ano. A participação do setor de árvores plantadas no PIB do Brasil é de 1,1% e (IBÁ 2014) de toda riqueza gerada no país, 5,5% é do PIB florestal, o que representa que a cada hectare plantado, o setor foi responsável por adicionar BRL 7,8 mil ao PIB nacional em 2014. Comparando com a soja, por exemplo, cuja produção brasileira é referencia mundial, os valores marginais foram de BRL 4,9 mil/ano por hectare plantado, 1,6 vezes menor que a contribuição do setor florestal de árvores plantadas, demonstrando a sua força mesmo em tempos de queda na economia. Ainda segundo a Ibá (2014), as exportações do setor atingiram receita de USD 8,49 bilhões, aumento de 2,5% em comparação ao recorde obtido em 2013 (USD 8,28 bilhões), reforçando a influência significativa do setor de florestas plantadas na economia do Brasil. As principais empresas do setor de florestas plantadas constantemente investem na busca por redução dos custos e melhoria nas estimativas dentro dos povoamentos florestais, proporcionando maior precisão na simulação dos orçamentos. Para isto, a modelagem de dados é essencial na definição da meta volumétrica periódica a ser colhida com finalidade principal de abastecimento 12 da fábrica. Portanto, comumente são aplicados modelos estatísticos de regressão para estimar o montante de madeira estocado nos plantios florestais. Em alguns casos tais modelos são não lineares e apresentam algumas peculiaridades, embora alguns tipos de modelos representem o comportamento biológico do crescimento das árvores e, portanto, são mais utilizados. Ainda possibilitam melhorar as estatísticas de resíduos trazendo maior acurácia às estimativas e consequentemente às etapas subsequentes do planejamento florestal, como a colheira, transporte de madeira e abastecimento da fábrica. Um forte contribuinte para o aumento no ônus da atividade florestal é o transporte da madeira até a fábrica. O transporte florestal, no Brasil, é realizado pelo modal rodoviário, principalmente, devido a grande malha viária e diferentes tipos de veículos, com várias capacidades de carga (MACHADO, LOPES e BIRRO, 2009). Este tipo de transporte corresponde a 59% da carga transportada, enquanto que o ferroviário representa 24%, no geral. No Canadá e Rússia, o transporte ferroviário apresenta mais força que o rodoviário, 46 e 81%, respectivamente (COSTA, 2006, citado por MACHADO, LOPES e BIRRO, 2009). Segundo Kariniemi (2005), 80% do transporte de madeira da Finlândia, país de economia fortemente ligada ao setor florestal, é realizado pelo modal rodoviário, 16% por ferroviário e os outros 4% do transporte é via aquático. Contudo, Tahvanainen e Anttila (2011) afirmam que o transporte ferroviário e aquaviário dependem do transporte rodoviário, da floresta até o terminal de recarga, mostrando a importância do modal rodoviário para outros modais neste país. Segundo Machado, Lopes e Birro (2000), o transporte ferroviário é o que possibilita maior deslocamento de carga, tornando esta modalidade mais eficiente e barata que as demais, entretanto o rodoviário tem mais incentivos por parte do governo, devido a políticas visando trazer multinacionais do ramo, como motivo principal. Segundo Epstein, Rönnqvist e Weintraub (2007), o 13 transporte corresponde a uma grande porção dos custos operacionais e, portanto, é necessário que seja bem gerenciado. Ainda segundo os autores, minimizar o custo do transporte é equivalente a minimizar o trabalho que ele demanda para a companhia, reduzindo, além disso, o impacto ambiental devido a emissões de gases do efeito estufa, ruído, poeira de estradas não pavimentadas, além do risco de acidente e redução de eficiência por congestionamentos. Portanto, a próxima parte desse trabalho aborda o tema geral de otimização aplicada à duas vertentes estudadas na engenharia florestal. Dentre elas, o Artigo 1 trata de ajustes de modelos não lineares aplicados á bases de dados florestais. Para proporcionar a parametrização dos modelos selecionados foi aplicado o algoritmo genético com intuito de fornecer os parâmetros iniciais para os métodos clássicos de ajuste (Gauss-Newton, Levenberg-Marquardt e NL2SOL). O Artigo 2 aborda a otimização multiobjetivo da atividade de transporte de madeira em situação real de empresa do setor florestal. Foi desenvolvido modelo matemático que simula o comportamento dos caminhões e gruas para o transporte de madeira dos talhões selecionados até a fábrica. 1.1 OBJETIVO O objetivo geral deste trabalho foi aplicar técnicas de otimização ao meio florestal, abordando problemas cotidianos do setor tanto por métodos aproximativos quanto por métodos exatos de resolução de problemas. 2 REFERENCIAL TEÓRICO 2.1 Regressão não linear 14 Análise de regressão estatística é frequentemente aplicada em muitos campos da ciência como biologia, pesquisas médicas e agricultura (PAN e FANG, 2002). Na área florestal, vários modelos clássicos de regressão linear têm sido aplicados, como encontrado em Wang et al. (2008), Wagner (2013), Spake e Doncaster (2017) e Groenendijk et al. (2017), abordando os modelos de Curtis, Clutter e o polinômio do quinto grau, por exemplo. Este tipo de análise linear tem objetivo de ajustar um modelo a uma amostra de dados derivada de uma população, utilizando princípios matemáticos, proporcionando um espectro generalista para a variável resposta. Modelos de regressão linear seguem a estrutura funcional entre o vetor dependente Y e o vetor independente X (Equação 1). E(Y) = f (X, β) (1) A estimativa de E(Y) é possível através da medição da tendência central para Y em relação a X por meio do vetor de parâmetros β na sua forma linear (LePAGE, 2011). O autor ainda afirma que a estimativa de Y é especificada pela função f que depende de parâmetros desconhecidos e o componente aleatório ε = Y - f (X, β), sobre o qual algumas pressuposições estatísticas são consideradas. Para diferentes instâncias observadas, assume-se que os resíduos são independentes e devem satisfazer Eε = 0, Eε = σ² > 0, isto é, os resíduos devem ter média zero e variância positiva e constante. Com objetivo de se obter a estimativa da variável de interesse, alguns métodos de parametrização são considerados clássicos na literatura, pelo seu frequente uso, o mais figurativo é o método dos mínimos quadrados ordinais (MQO). Segundo Pan e Fang (2002), no MQO o vetor de parâmetros β de um modelo multivariado linear é dado pelo estimador de mínimos quadrados , conforme (2). 15 = YX’(X’X)-1 (2) deve obter o quadrado mímino dos resíduos, traço, determinante e maior valor de colinearidade da soma de quadrados multivariada (PAN e FANG, 2002). Assim sendo, se MQO é considerado ótimo no que se refere à minimizar a covariância do traço [min tr(Cov( ))], então ele é chamado de melhor estimador linear inviezado (MELI) ou estimador de Gauss-Markov, produzindo os valores ótimos para cada parâmetro. Semelhante aos modelos de regressão linear, os não lineares apresentam a variação principal de um comportamento não linear de pelo menos um dos parâmetros em função de variável resposta (BATES and WATTS, 1988). A grande vantagem desta abordagem advem da possibilidade de melhorias na análise estatística com relação às estatísticas de resíduos e ajuste da função aos dados. Outra vantagem importante no meio florestal é o comportamento biológico que alguns modelos não lineares apresentam, importante nas estimativas de crescimento e produção de maciços florestais. Em ciências florestais, vários trabalhos são publicados anualmente abordando algum tipo de análise estatística não linear (RETSLAFF et al., 2015; MEHTÄTALO et al., 2015; MISIK et al., 2016; FISCHER e SCHÖNFELDER, 2017). Diferentemente da análise linear, a regressão não linear não garante a solução ótima na estimativa do vetor de parâmetros, portanto métodos aproximativos são aplicados com objetivo de obter soluções próximas ao valor ótimo em menor tempo possível. Pázman (2011) afirma que o método iterativo de Gauss-Newton é o mais utilizado para estimar o vetor de parâmetros em análises de regressão não linear. Contudo, o autor afirma que este método resulta em não convergência para alguns casos, demandando a escolha dos parâmetros iniciais com maior parcimônia, sendo esta uma desvantagem deste tipo de 16 análise comparada com a regressão linear. Jiang e Liu (2011) afirmam que o método de Levenberg-Marquardt proporciona maior flexibilidade do que Gaus- Newton na estimativa, resultando em maior frequência de convergências com a mesma confiança estatística. Como alternativa a estes métodos de parametrização clássicos, emergem alguns algoritmos de busca e refinamento, também chamados evolucionários. O maior expoente desta classe de algoritmos é o algoritmo genético, desenvolvido por Holland (1973) e massivamente aplicado em trabalhos de várias áreas, inclusive na estatística (Yao e Sethares, 1994; Chatterjee et al., 1996; Karr et al., 1998; Kapanoglu et al., 2007; Chen, 2015; Jin et al., 2016). 2.2 Algoritmo genético Os primeiros estudos abordando os conceitos de algoritmos adaptativos e evolutivos foram introduzidos pelos trabalhos independentes de Holland (1962) e Bremermann (1962), mas já estudados por Fraser (1957). Paralelamente, Fogel (1962) trouxe o tema de programação evolutiva em seu trabalho Autonomous Automata, bem como Rechenberg (1975) abordou o tema estratégias evolutivas, reforçando o conceito de algoritmos naturais, compilados por Victorino (2005). O algoritmo genético (AG), principal expoente da classe, foi abordado primeiramente por Holland (1973), em seu trabalho Genetic algorithms and the optimal allocation of trials, ao qual afirma a existência de uma classe de problemas de otimização caracterizados pela (1) complexidade e incerteza inicial; (2) necessidade de adquirir novas informações rapidamente para reduzir a incerteza e (3) requisitos das novas informações a serem exploradas e adquiridas com performances médias crescentes e taxa de aquisição de informações constantes. Com esta explanação, o autor descreve a problemática da execução do AG, abordando ainda o conceito de estruturas 17 passíveis de busca, isto é, os cromossomos, testados por uma função de aptidão, ou função fitness, utilizando os termos atuais (Figura 1). O algoritmo genético é formado genericamente pelas seguintes partes de seu processo: a) Codificação: É a forma cujos dados brutos serão avaliados. Segundo Gomide (2009), um cromossomo formado por um número de genes, onde seus alelos assumem valores binários, podem produzir uma estrutura que pode assumir 2n diferentes tipos de esquemas, sendo que n é o número alelos deste gene. Mitchell (1996), afirma que a codificação binária é a mais utilizada neste algoritmo, contudo, há outras opções como com números reais e codificação Gray (VICTORINO, 2005). b) Função de Aptidão (fitness): É a conexão do algoritmo com o problema real a ser resolvido (SAKAWA, 2002). Segundo Gomide (2009), o fitness é a expressão fenotípica de um indivíduo, calculado a partir da combinação dos genes presentes no cromossomo com objetivo de medir a qualidade do indivíduo. c) População: Divide o problema em um conjunto de soluções factíveis contidas compondo o cromossomo e definidos numa população com dimensão pré-definida. Esta estratégia permite maior probabilidade de o algoritmo sair de ótimos locais (GOMIDE, 2009). Segundo Oliveira (2004), a definição do tamanho da população inicial é fundamental para boa convergência do algoritmo, sendo que os fatores que atuam na população são regras baseadas em seleção, reprodução e mutação (HERTZ e WIDMER, 2003). 18 d) Operador de Seleção: É o processo de selecionar os melhores indivíduos para a próxima geração. O processo de seleção deve ser balanceado com a mutação e reprodução, no processo de parametrização, para proporcionar melhoria no fitness médio a cada geração (MITCHELL, 1996). Existem vários tipos de métodos de seleção, segundo Gomide (2009), são eles: Elitista, proporção, roleta viciada, dentre outros. Contudo, o autor comenta que o método mais aplicado é o da roleta viciada. e) Operador de Reprodução: Também chamado de crossover, o método consiste na produção de dois novos indivíduos obtidos pela troca gênica entre seus pais (ROTHLAUF, 2006). Segundo Reeves (2003), o mecanismo é aplicado a partir de escolhas compostas por regras aleatórias que são responsáveis pela recombinação de genes entre os cromossomos. f) Operador de Mutação: Segundo Victorino (2005), este operador é responsável pela introdução e manutenção da diversidade genética da população, proporcionando maior área de busca no espaço de soluções. A diferença entre mutação e reprodução é a característica aleatória da mutação tentando a melhoria de soluções individuais, enquanto a reprodução percorre o espaço factível também de maneira aleatória, mas direcionada (MITSUO e CHENG, 2000). Contudo, segundo Mitsuo e Cheng (2000), a mutação permite maior eficiência na exploração da vizinhança de um cromossomo devido sua característica de realizar alterações mínimas na estrutura cromossômica. 19 Figura 1: Exemplificação do fluxo de execução simplificada do algoritmo genético Fonte: Razali e Geraghty (2011). Ao longo das últimas décadas, diversos trabalhos abordaram para as mais diversas finalidades no ramo de otimização, como por exemplo, uso do AG multiobjectivo para otimização de sistema de distribuição de água (TANYIMBOH e SEYOUM, 2016) e solução de problemas de agendamento (KIM, WALEWSKI e CHO, 2016), uso do AG e SIG para otimização de alocação de terrenos urbanos (LI e YEH, 2007), otimização de portfólio de investimentos ajudando na tomada de decisão (ZAINASHEV, 2011) e até aplicado para otimizar etapas de outro mecanismo de inteligência artificial, como é o caso do trabalho de Aydogan, Karaoglan e Pardalos (2012), cujos autores utilizaram o AG para otimizar o algoritmo de classificação baseado em regras fuzzy para problemas de grande dimensão. Amplamente utilizado em 20 problemas de roteamento de veículos multiobjetivo (CHAND et al., 2010; CHAND e MOHANTY, 2013), roteamenteo com múltiplo depósitos (LAU et al., 2010; WANG et al., 2016), abordando problemas online de roteamento de veículos (ABDALLAH, ESSAM e SARKER, 2017) e testando variações nos operadores genéticos para o problema de roteamento de veículos (PULJIC´ e MANGER, 2013; VAIRA e KURASOVA, 2013; PIERRE e ZAKARIA, 2017). 2.3 Problema de roteamento de veículos 2.3.1 Teoria de grafos No meio científico inúmeras situações reais podem convenientemente ser descritas por diagramas, formados por um conjunto de pontos unidos por linhas, agrupando determinados pares destes pontos. A abstração matemática desta situação remete ao conceito de grafos. Bondy e Murty (1976) descrevem matematicamente um grafo G como constituído pela ordem tripla (V(G), E(G), ψG), sendo V(G) o conjunto não vazio de vértices, E(G) o conjunto de arestas e ψG o funcional que associa cada aresta de G à um par ordenado de vértices do grafo G. A Figura 2 ilustra um grafo contendo as estruturas apresentadas. 21 Figura 2. Representação de estrutura em grafo. Fonte: Bondy and Murty (1976). Operacionalmente a teoria dos grafos é amplamente utilizada em rotas terrestres, aplicados à logística do transporte de mercadorias, visando reduzir a distância percorrida e, consequentemente, o custo da atividade. Um exemplo aplicado é o trabalho de Pereira (2008), que utilizou o conceito de grafos e redes de transporte para executar a gestão da distribuição de bebidas em perímetro urbano. Outro exemplo da aplicação é o trabalho de Araújo Júnior (2016), ao qual propõe um sistema multiagente para otimizar o transporte florestal embasado em redes de transporte, as quais deveriam ter suas distâncias minimizadas para maior eficiência da operação. Reduzir o percurso em uma rede de pontos é uma tarefa complexa, um problema clássico de redes interligadas é o do caixeiro viajante, cujo objetivo é percorrer todos os pontos do grafo produzindo o menor caminho possível. O problema do caixeiro viajante (PCV) clássico apresenta, segundo Dantzig, Fulkerson e Johnson (1954), formulação matemática composta por 3 restrições e uma função objetivo, que minimiza a soma das distâncias percorridas. Utilizado por Almoustafa, Hanafi e Mladenovi´c (2013), o PCV é alicerce para inferir acerca do problema de roteamento veicular, comumente utilizado por empresas em redes logísticas de transporte de bens. 22 Inúmeros trabalhos têm sido publicados buscando otimizar o caminhamento em redes de pontos, como é o caso de Moon et al. (2002), sendo aplicado o algoritmo genético com topologia sorteada, definida pela ordenação dos vértices em grafo orientado, para minimizar a soma das distâncias percorridas no problema do caixeiro viajante com restrições precedentes, sendo um dos problemas combinatoriais mais complexos. O trabalho consistiu em comparar o algoritmo genético proposto, com o algoritmo branch-and-bound e algoritmo genético clássico, formado por estratégias de penalidades, o que mostrou que a nova estratégia foi superior aos algoritmos tradicionais. Outro exemplo de trabalho aplicado à teoria dos grafos foi observado em Nazif e Lee (2015), em qual os autores otimizaram os operadores do algoritmo genético aplicado ao roteamento de veículos capacitado, cujos resultados apontaram que a modificação proposta foi competitiva, quando comparada com heurísticas e meta-heurísticas aplicadas ao mesmo problema, como o algoritmo de busca tabu proposto por Toth e Vigo (2003), o algoritmo savings ants (SA) proposto por Reimann, Doerner e Hartl (2004), memória genética adaptativa (SEPAS), proposto por Tarantilis (2005) e variações do algoritmo genético propostas por Prins (2004). 2.3.2 Problema de roteamento de veículos clássico Originalmente proposto por Dantzig e Ramser (1959), com o nome de Truck Dispatching Problem, o problema de roteamento veicular (PRV) generaliza o problema do caixeiro viajante (PCV), como afirmam os autores, os quais foram os primeiros a solucionar este tipo de problema. Segundo seus os autores, o problema do caixeiro viajante consiste em determinar o caminho mais curto que passa por pontos determinados uma única vez. Já o PRV estende o PCV de maneira a especificar que devem ser realizadas determinadas entregas 23 em cada ponto considerado no grafo, exceto o ponto final. Sendo que a capacidade de carga do veículo é menor do que a demanda de produtos a serem entregues, fazendo com que o veículo deva executar n vezes o deslocamento característico do caixeiro viajante. Os autores ainda afirmam que este tipo de problema apresenta complexidade natural de resolução pelo número elevado de combinações a serem testadas, sendo que um número relativamente pequeno de pontos, 15, produza um número elevadíssimo de rotas a serem testadas, 653.837.184.000, sendo que a representação matemática da permutação das rotas combinadas é n!/2, onde n é o número de vértices de um grafo. Clarke e Wright (1964) encontraram melhor solução para este problema frente aos seus idealizadores, confirmando a expansão do desenvolvimento desses algoritmos de resolução. Este tipo de problema é classificado por Lenstra e Rinnooy Kan (1981) como NP-hard, sendo de difícil resolução e, dependendo de sua complexidade, não são solucionados pelos métodos exatos em tempo operacional para a tomada de decisão. Entretanto, com o surgimento de cenários mais complexos, há busca por outros métodos de resolução mais eficientes, como nos trabalhos de Stepanov e Smith (2009) e Subramanian et al. (2012), cujos autores utilizam algoritmo híbrido e roteamento multi-objetivo para definir as melhores rotas para veículos em uma rede de transporte, evidenciando a necessidade de algoritmos mais robustos para determinadas situações. A Figura 3 é uma representação do problema de roteamento de veículos com 1 depósito e 5 consumidores, cujo um objetivo possível é percorrer, uma única vez, o menor caminho entre os pontos, realizando coleta/entrega de material nestes pontos. Este tipo de problema é formado por múltiplos caixeiros viajantes, contudo, no setor florestal, há uma adaptação no conceito do problema, sendo que o veículo não executa o percurso apenas uma vez, pelo fato do volume a ser retirado de cada talhão é maior do que a capacidade de carga do veículo. 24 Figura 3: Representação esquemática do roteamento de veículos. Fonte: Surekha e Sumathi (2011). Algumas variações do problema clássico de roteamento de veículos, como em problemas com entrega e coleta simultânea, vêm sendo utilizadas para testar novos algoritmos. O problema de entrega e coleta simultâneas estende o problema de roteamento de veículos transportando bens de um depósito até um cliente, mas também de clientes para o depósito (WASSAN e NAGY, 2014). Trabalhos como o de Subramanian et al. (2010) implementaram heurística paralela para o problema de roteamento com entrega e coleta simultânea. A mesma variação do problema de roteamento de veículos foi utilizada por Jin Ai e Kachitvichyanukul (2009) e Gajpal e Abad (2009), na qual os autores aplicaram particle swarm optimization (PSO) e um sistema de colônia de formigas, respectivamente. Segundo Salhi e Nagy (1999), no problema de roteamento de veículos com backhauling, ou coleta e entrega simultânea, não é somente requerida a entrega da mercadoria aos clientes, mas também o carregamento dos veículos para retorno ao depósito. Ainda segundo os autores, este tipo de problema, na prática, deve ser realizado após todas as entregas serem efetuadas, pois pode dificultar o rearranjo de mercadoria no veículo e produzir um percurso 25 infactível do ponto de vista econômico. Este é mais uma variação do problema clássico de roteamento de veículos, cujo alguns trabalhos se destacam, como é o caso de Surekha e Sumathi (2011) e Ho et al. (2008), cujos autores trabalharam com algoritmo genético na solução de problemas de roteamento com múltiplos depósitos. A combinação de variações do problema de roteamento também é comum no meio científico, tentando trazer o máximo de realidade nas análises, como no trabalho de Nagy e Salhi (2005), que propuseram um algoritmo heurístico que assume entrega e coleta de mercadorias de maneira simultânea, diferentemente da literatura comum, considerando ainda um único depósito e múltiplos depósitos no experimento. Uma situação comum em organizações que se utilizam de transporte de mercadorias é a heterogeneidade da frota de veículos, no que se refere à capacidade de carga de cada tipo de veículo empregado, caracterizando mais uma variação no problema clássico de roteamento. Segundo Wei, Zhang e Lim (2014), a operação de carregamento e entrega de mercadorias em clientes dispersos geograficamente, a partir de um ponto central de distribuição, representa alta porção dos custos correntes de uma empresa, considerando ainda a disposição da carga no interior do veículo e frota mista de veículos, que proporciona maior flexibilidade no planejamento do custo efetivo do transporte. Alguns algoritmos são aplicados para resolver este tipo de problema com tempo hábil de processamento, tendo em vista que o método exato não produz este tipo de resultado para problemas de grande dimensão, como exemplo podem ser citados os trabalhos de Yousefikhoshbakht, Didehvar e Rahmati (2015), Low et al. (2014), Jiang et al. (2014) e Salhi, Imran e Wassan (2014), nos quais os autores aplicam uma série de algoritmos heurísticos e meta-heurísticos para resolução de problemas de roteamento de veículos com frota mista. Na área florestal alguns trabalhos vêm sendo publicados nos últimos anos utilizando o modelo de roteamento de veículos como objeto de estudo, 26 como nos casos Carlsson e Rönnqvist (2007), cujos autores abordaram o planejamento tático de transporte de madeira utilizando programação linear com fluxo constante entre pontos de demanda e entrega (backhauling). Outro exemplo é o trabalho de Audy, D’Amours e Rönnqvist (2012) abordando métodos de planejamento e sistemas de suporte de decisão para o roteamento de veículos aplicado ao transporte de madeira, em que os autores realizam uma revisão sobre o tema apresentando abordagens para o meio florestal. Outros exemplos na área florestal são encontrados em Weintraub et al. (1996), Doerner et al. (2007) e Flisberg, Lidén e Rönnqvist (2007). 2.4 Otimização Multiobjetivo A complexidade dos problemas operacionais observados em empresas do setor florestal são de difícil representação matemática e, por isso, geralmente são executados de maneira empírica. Um dos fatores que mais dificulta sua representação é sua característica multiobjetivo pela interação de padrões difusos aos quais a função objetivo do problema pode seguir. O termo otimização multiobjetivo é sinônimo de termos comuns na otimização, como vetor de otimização, otimização multicritério ou otimização multiperformance (COELHO, 2000). Segundo Osyczka (1985), citado por Fotakis et al. (2012), otimização multiobjetivo é um problema de busca na qual um vetor de variáveis de decisão é considerado para satisfazer as restrições e otimiza o vetor de sua função, sendo o objetivo de cada vetor a função objetivo do problema. Uma série de trabalhos foram realizados ao longo dos anos aplicando ao setor florestal o conceito de otimização multiobjetivo (MARINGANTI et al. 2011; FOTAKIS et al. 2012; GROOT et al. 2012; TÓTH et al. 2013; CAMBERO e SOWLATI 2016). 27 Alguns autores apresentam como principais abordagens para a definição da função multiobjetivo a aplicação de pesos nas sub-funções (TOTH et al. 2013; CAMBERO e SOWLATI 2016) e aplicação do Ranqueamento de Pareto (MARINGANTI et al. 2011; FOTAKIS et al. 2012; GROOT et al. 2012). Segundo Fotakis et al. (2012), o melhor método multiobjetivo de otimização é a aplicação do Ranqueamento de Pareto, afirmando como sendo sua maior vantagem o uso das unidades reais de cada objetivo, evitando sua conversão e consequentemente aplicação de pesos, que muitas vezes derregulam o funcionamento do modelo. A Figura 4 esquematiza o comportamento das soluções dentro do sistema de Ranqueamento de Pareto Figura 4: Comportamento das soluções de um problema multiobjetivo de acordo com o sistema de Ranqueamento de Pareto. Fonte: Monfared et al. 2012. Contudo, a conversão dos objetivos em unidades monetárias, quando possível, justifica sua aplicação no que se refere a análise de decisão sobre a simulação produzida pelo modelo. Segundo Cambero e Sowlati (2016), quando se aplica pesos aos sub-objetivos o fator de preferência é inserido ao sistema de 28 tomada de decisão, sendo que cada objetivo pode eventualmente ser mais almejado que outro em determinadas situações operacionais. 3 CONSIDERAÇÕES GERAIS Aplicação de modelos, estocásticos ou determinísticos, representa um avanço no campo científico e prático pela possibilidade de representar situações reais de maneira matemática. Modelos são representações da realidade e, por muitas vezes, não são capazes de forcener o nível de detalhe requerido ou a precisão desejada em alguma análise. Contudo, emerge como um bom balizador na tomada de decisão principalmente na gestão de empresas e em todos os níveis de planejamento. A grande vantagem dessa aplicação é analisar sistemas reais de maneira a proporcionar vantagens econômicas, sociais e ambientais, através da estimativa de dados, redução de custos, emissão de gases CFC, impactos físicos, químicos, biológicos e sociais, maximização de receitas, da produção etc. REFERÊNCIAS ABDALLAH, A.; M.; F.; M.; ESSAM, D.; L.; SARKER, R.; A. On solving periodic re-optimization dynamics vehicle routing problems. Applied Soft Computing, Amsterdam, v.55, p.1-12, 2017. ALMOUSTAFA, S.; HANAFI, S.; MLADENOVI´C, N. New exact method for large asymmetric distance-constrained vehicle routing problem. European Journal of Operational Research, Kidlington, v.226, p.386-389, 2013. ARAÚJO JÚNIOR, C.; A. Um sistema multiagente para otimização do transporte florestal. Tese de doutorado. Universidade Federal de Viçosa, Viçosa, Brasil, 2016. AUDY, J.;F.; D’AMOURS, S.; RÖNNQVIST, M. Planning methods and decision support systems in vehicle routing problems for timber transportation: 29 A review. Centre Interuniversitaire de Recherche sur les Rèseaux d’entreprise, la Logistique et le Transport (CIRRELT), Québec, v.38, 2012. AYDOGAN, E.; K.; KARAOGLAN, I.; PARDALOS, P.; M. hGA: Hybrid genetic algorithm in fuzzy rule-based classification systems for high- dimensional problems. Applied Soft Computing, Amstedam, v.12, p.800-806, 2012. BATES, D.M., WATTS, D.G. Nonlinear regression analysis and its applications. Wiley, Amsterdam, pp.1-90, 1988. BONDY, J. A.; MURTY, U. S. R. Graph theory with applications. Departamento de Combinação e Otimização. Universidade de Waterloo, Ontario, Canada, editora North Holland. Capítulo 1, p.264, 1976. BREMERMANN, H.; J. Optimization through evolution and recombination. Editor(s): M.C. Yovits, G.T. Jacobi & G.D. Goldstine, Spartan Book, Self- Organizing Systems, Washington, pp. 93-106, 1962. CAMBERO, C.; SOWLATI, T. Incorporating social benefits in multi-objective optimization of forest-based bioenergy and biofuel supply chain. Applied Energy, Amsterdam, v.178, p.721-735, 2016. CARLSSON, D.; RÖNNQVIST, M. Backhauling in forest transportation: Models, Methods, and practical usage. Canadian Journal Forest Research, Ottawa, v.37, n.12, p.2612-2623, 2007. CARVALHO, M. B. de. Aplicações de meta-heurística genética e fuzzy no sistema de colônia de formigas para o problema do caixeiro viajante. Dissertação de Mestrado. Universidade Estadual de Campinas, Campinas, Brasil, 2007. CATTARUZZA, D.; ABSI, N.; FEILLET, D.; VIDAL, T. Memetic algoritmo for multi trip vehicle routin problem. European Journal of Operational Research, Kidlington, v.236, p.833-848, 2014. CHAND, P.; MISHRA, B.; S.; P.; DEHURI, S. A multi objective genetic algorithm for solving vehicle routing problem. Internatinal Journal of Information Technology and Knowledge Management, Netherland, v.2, n.2, p.503-506, 2010. 30 CHAND, P.; MOHANTY, J.; R. Multi objective genetic approach for solving vehicle routing problem. International Journal of Compurter Theory and Engineering, Ottawa, v.5 n.6, p.846-849, 2013. CHATTERJEE, S., LAUDATO, M., LYNCH, L.A. Genetic algorithm and their statistical applications: an introdution. Computational Statistics & Data Analysis, v.22, p.633-651, 1996. CHEN, J. A new hybrid genetic algorithm for parameter estimation of nonlinear regression modeling. Long, S., Dhillon, B.S. (eds), Proceedings of the 15th International Conference on Man-Machine-Enviroment System Engeneering, Lecture Notes in Electrical Engineering, 2015 CLARKE, G.; WRIGHT, J.W. Scheduling of vehicle from a central depot to a number of delivey points. Operations Research, Catonsville, v.12, p.568-581, 1964. COELLO, C.A. An updated survey of GA-basedmultiobjective optimization techniques. ACM Computing Surveys (CSUR), Boston, v.32, n.2, p.109–143, 2000. DANTZIG, G. B.; FULKERSON, D. R.; JOHNSON, S. M. Solution of a large- scale traveling salesman problem. Operations Research, Catonsville, v.2, p.393-410, 1954. DANTZIG, G. B; RAMSER, J. H. The truck dispatching problem. Management Science, New York, v.6, n.1, pp. 80-91, 1959. EPSTEIN, R.; RÖNNQVIST, M.; WEINTRAUB, A.. Forest transportation (Chapter 20). Editor(s): WEINTRAUB, A.; ROMERO, C.; BJØRNDAL, T.; EPSTEIN, R; MIRANDA, J. Handbook of operations research in nature resourcers. Springer, 1 ed, p.611, 2007. FISCHER, C.; SCHÖNFELDER, E. A modefied growth function with interpretable parameters applied to the age-height relationship of individual trees. Canadia Journal of Forest Research, v.47, n.2, p.166-173, 2017. FLISBERG, P.; LIDÉN, B.; RÖNNQVIST, M. A hybrid method based on linear programming and tabu search for routing of logging trucks. Computer & Operational Research, Netherland, v.36, n.4, p.1122-1144, 2007. FOGEL, L.; J. Autonomous automata. Industrial Research, v.4, p.14-19, 1962. 31 FOTAKIS, D. G.; SIDIROPOULOS, E.; MYRONIDIS, D.; IOANNOU, K. Spatial genetic algoritmo for multi-objective forest planning. Forest Policy and Economics, Amstedam, v.21, p.12-19, 2012. FRASER, A.; S. Simulation of genetic systems by automatic digital computers: I. Introduction. Australian Journal Biololy Science, Clayton, v.10, p.484-491, 1957. GAJPAL, Y.; ABAD, P. An ant colony system (ACS) for vehicle routing problem with simultaneous delivery and pick-up. Computers and Operational Research, London, v.36, p.3215-3223, 2009. GOMIDE, L.; R. Planejamento espacial florestal. Teses de Doutorado. Universidade Federal do Paraná, Curitiba, 2009. GROENENDIJK, P.; BONGERS, F.; ZUIDEMA, P.; A. Using tree-ring data to improve timber-yield projections for African wet tropical forest tree species. Forest Ecology and Management, v.400, p.396-407, 2017. GROOT, J.C.J.; OOMEN, G.J.M.; ROSSING, W.A.H. Multi-objective optimization and design of farming systems. Agricultural Systems, Amsterdam, v.110, p.63-77, 2012. HERTZ, A.; WIDMER, M. Guidelines for the use of meta-heuristcs in combinatorial optimization. European Journal of Operational Research, Kidlington, v.151, n.2, p.247-252, 2003. HO, W.; HO, G.; T.; S.; JI, P.; LAU, H.; C.; W. Hybrid genetic algorithm for the multiple-depot vehicle routing problem. Engineering Applications of Artificial Intelligence, Amsterdam, v.21, p.548-557, 2008. HOLLAND, J.; H. Genetic algorithms and the optimal allocation of trials. Society for Industrial and Applied Mathematics Journal on Computing, Philadelphia, v.2, n.2, p.88-105, 1973. HOLLAND, J.; H. Outline for a logical theory of adaptive systems. Journal of Association for Computing Machinery, New York, v. 3, p. 297-314, 1962. HOLLAND, J.; H. Outline for a logical theory of adaptive systems. Journal of Association for Computing Machinery, New York, v. 3, p. 297-314, 1962. 32 IBÁ. Indústria Brasileira de Árvores, p.77, 2015. JIANG, J.; NG, K.; M.; POH, K.; L.; TEO, K.; M. Vehicle routing problem with heterogeneous fleet and time windows. Expert Systems with Applications, Amsterdam, v.41, p.3748-3760, 2014. JIANG, L.;C.; LIU, R.;L. Segmented taper equations with crown ratio and stand density for Da-hurian Larch (Larix gmelinii) in Northeastern China. Journal of Forestry Research, Ottawa, v.22, n.3, p.347-352, 2011. JIN AI, T.; KACHITVICHYANUKUL, V. A particle swarm optimization fot the vehicle routing problem with simultaneous puckup and delivery. Computers and Operational Research, London, v.36, p.1693-1702, 2009. JIN, Y.F., YIN, Z.Y., SHEN, S.L., HICHER, P.Y. Selection of sand models and identification of parameters using an enhanced genetic algorithm. International Journal for Numerical and Analytical Methods in Geomechanics, v.40, p.1219-1240, 2016. KAPANOGLU, M., KOC, I.O., ERDOGMUS, S. Genetic algorithms in parameter estimation for nonlinear regression models: an experimental approach. Journal of Statistical Computation and Simulation, v.77, n.10, p.851-867, 2007. KARINIEMI, A. Puunkorjuu ja kaukokuljetus vuonna. Metsätehon Katsaus, v.19, n.4, 2005. KARR, C.L., WECK, B., FREEMAN, L.M. Solutions to systems of nonlinear equations via a genetic algorithm. Engineering Applications of Artificial Intelligence, v.11, p.369-375, 1998. KIM, K.; WALEWSKI, J.; CHO, Y.; K. Multiobjective construction schedule optimization using modified niched pareto genetic algorithm. Journal of Management in Engineering, Reston, v.32, n.2, p.1-12, 2016. LAU, H.; C.; W.; CHAN, T.; M.; TSUI, W.; T.; PANG, W.; K. Application of genetic algorithms to solve the multidepot vehicle routing problem. IEEE Transections Automation Science and Engineering, v.7, n.2, p.383-392, 2010. LENSTRA, J.; RINNOOY KAN, A. Complexity of vehicle routing and scheduling problems. Networks, Malden, v.11, p.221-227, 1981. 33 LePAGE, R. Linear regressin models. LOVRIC, M. (Ed.) International encycopedia of statistical science. Springer, p.1673, 2011. LI, X.; YEH, A.; G.O. Integration of genetic algorithms and GIS for optimal location search. International Journal of Geographical Information Science, London, v. 19, n.5, p.581-601, 2007. LOW, C.; CHANG, C.;M.; LI, R.; K.; HUANG, C.; L. Coordination of production scheduling and delivery problems with heterogeneous fleet. International Journal of Production Economics, Amsterdam, v.153, p.139- 148, 2014. MACHADO, C.; C.; LOPES, E.; S.; BIRRO, M.; H.; B. Elementos básicos do transporte florestal rodoviário. Editora UFV. Universidade Federal de Viçosa, Viçosa. p.167, 2 ed., 2009. MACHADO, C.; C.; LOPES, E.; S.; BIRRO, M.; H.; B. Elementos básicos do transporte florestal rodoviário. Editora UFV. Universidade Federal de Viçosa, Viçosa. p.167, 1 ed., 2000. MARINGANTI, C.; CHAUBEY, I.; ARABI, M.; ENGEL, B. Application of a multi-objective optimization method to provide least cost alternatives for NPS pollution control. Environmental Management, Gewerbestrasse, v.48, p.448- 461, 2011. METAHTÄTALO, L.; DE-MIGUEL, S.; GREGOIRE, T. Modeling height- diameter curves for prediction. Canadian Journal of Forest Research, Ottawa, v.45, p.826-837, 2015. MISIK, T.; ANTAL, K.; KÁRÁSZ, I.; TÓTHMÉRÉSZ, B. Non-linear height- diameter models for three woody, understory species in a temperate oak forest in Hungary. Canadian Journal of Forest Research, v.46, n.11, p.1337-1342, 2016. MITCHELL, W. An introduction to genetic algorithms. Bradford Book, p.158, 1996. MITSUO, G.; CHENG, R. Genetic algorithms and engineering optimization. Wiley, p.512, 2000. 34 MONFARED, M.D.; MOHADES, A.; REZAEI, J. Convex hull ranking algorithm for multi-objective evolutionary algorithms. Scientia Iranica, Sharif, v.18, n.6, p.1435-1442, 2011. MOON, C.; KIM, J.; CHOI, G.; SEO, Y. An efficiente genetic algorithm for the traveling salesman problem with precedent constraints. European Journal of Operational Research, Kidlington, v.140, p.606-617, 2002. NAGY, G.; SALHI, S. Heuristic algorithms for single and multiple depot vehicle routing problem with pickup and deliveries. European Journal of Operational Research, Kidlington, v.162, p.126-141, 2005. NAZIF, H.; LEE, L..; S. Optimised crossover genetic algorithm for capacitated vehicle routing problem. Applied Mathematical Modelling, Amsterdam, v.36, p.2110-2117, 2012. OLIVEIRA, S.; A. Metaheurísticas aplicadas ao planejamento da expansão da transmissão de energia elétrica em ambiente de processamento distribuído. Tese de Doutorado. Universidade Estadual de Campinas, Campinas, 2004. OSYCZKA, A. Multicriteria optimization for engineering design. In: Gero, J.S. (Ed.), Design Optimization. Academic Press, Cambridge, MA, p.193–227, 1985. PAN, J.;X.; FANG, K.;T. Growth curve models and statistical diagnostics. Springer, Amsterdam, p.387, 2002. PÁZMAN, A. Nonlinear regression. LOVRIC, M. (Ed.) International encycopedia of statistical science. Springer, Amsterdam, p.1673, 2011. PEREIRA, G.; R. Aplicação da gestão baseada em atividades à distribuição urbana de bebidas. Tese de doutorado. Universidade Federal do Rio de Janeiro, Rio de Janeiro, Brasil, 2008. PIERRE, D.; M.; ZAKARIA, N. Stochastic partially optimized cyclic shift crossover for multi-objective genetic algorithms for the vehicle routing problem with time-windows. Applied Soft Computing, Amsterdam, v.52, p.863-876, 2017. 35 PRINS, C. A simple and effective evolutionary algorithm for the vehicle routing problem, Computers and Operations Research, Amsterdam, v.31, p.1985– 2002, 2004. PULJIC´, K.; MANGER, R. Comparison of eight evolutionary crossover operators for the vehicle routing problem. Mathematical Communications, Osijek, v.18, p.359-375, 2013. RASHIDI, H.; FARAHANI, R.Z A hybrid ant colony system for partially dynamics vehicle routing problem. American Journal of Operational Research, Rosemead, v.2, n.3, p.31-44, 2012. RAZALI, N.; M.; GERAGHTY, J. Genetic algorithm performance with different selection strategies in solving TSP. Proceedings of the World Congresso on Engineering, London, vol II, July 6-8, 2011. RECHENBERG, I. Evolutionsstrategies – Optimierung technischer systeme nach prinzipien der biologischen evolution. Feddes Repertorium, Stuttgart, v.86, n.5, p.337-340, 1975. REEVES, C.; R. Handbook of metaheuristics. Editor(s): GLOVER, F.; KOCHENBERGER, G.; A. Kluwer Academic Publicer, p.557, 2003. REIMANN, M. DOERNER, R.; HARTL, R. D-ants: saving based ants divided and conquer the vehicle routing problem. Compute and Operational Research, Amsterdam,v.31, p.563-591, 2004. REINA, C. D. Roteirização de veículos com janelas de tempo utilizando algoritmo genético. Dissertação de Mestrado. Escola Politécnica da Universidade de São Paulo, São Paulo, Brasil, 2012. RETSLAFF, F.;A.;S.; FILHO, A.;F.; DIAS, A.;N.; BERNETT, L.;G.; FIGURA, M.;A. Curvas de sítio e relações hpsométricas para Eucalyptus grandis na região dos Campos Gerais, Paraná. Cerne, v.21, n.2, p.219-225, 2015. ROTHLAUF, F. Representation fo genetic and evolucionary algorithms. Springer, p.325, 2006. SAKAWA, M. Genetic algorithm and fuzzy multiobjetive optimization. Kluwer Academic Publicers, Springer US, Book, v.14, n.1, p.288, 2002. 36 SALHI, S.; IMRAN, A.; WASSAN, N.; A. The multi-depot vehicle routing problem with heterogeneous vehicle fleet: Formulation and a varible neighborhood search implementaion. Computers and Operational Research, London, v.52, p.315-325, 2014. SALHI, S.; NAGY, G. A cluster insertion heuristic for single and multiple depot vehiclo routing problem with backhauling. Journal of Operational Research Society, London, v.50, p.1034-1042, 1999. SOLOMON, M. M. Algoritms for vehicle routing and scheduling problems with time window constraints. Operations Research, Catonsville, v.35, n.2, p.254- 266, 1987. SPAKE, R.; DONCASTER, C.; P. Use of meta-analysis in forest biodiversity research: key challenges and considerations. Forest Ecology and Management, v.400, p.429-437, 2017. STEPANOV, A.; SMITH, J. M. Multi-objetive evacuation routing in transportation network. European Journal of Operational Research, Kidlington, v.198, p.435-446, 2009. SUBRAMANIAN, A.; DRUMMOND, L.; M.; A.; BENTES, C.; OCHI, L.; S.; FARIAS, R. A parallel heuristic for the vehicle routing problem with simultaneous puckup and delivery. Computers and Operational Research, London, v.37, p.1899-1911, 2010. SUBRAMANIAN, A.; PENNA, P.H.V.; UCHOA, E.; OCHI, L.S. A hybrid algoritmo for heterogeneous fleet vehicle routing problem. European Journal of Operational Research, Kidlington, v.221, p.285-295, 2012. SUREKHA, P.; SUMATHI, S. Solution to multi-depot vehicle routing problem using genetic algorithms. World Applied Programming, Dadaab, v.1, n.3, p.118-131, 2011. TAHVANAINEN, T.; ANTTILA, P. Supply chain cost analysis of long-distance transportation of energy wood in Finland. Biomass and Bioenergy, v.35, p.3360-3375, 2011. TANYIMBOH, T.; T.; SEYOUM, A.; G. Multiobjctive evolutionary optimization of water distribution systems: Exploiting diversity with infeasible solutions. Journal of Environmental Management, Amsterdam, v.183, p.133- 141, 2016. 37 TARANTILIS, C.; D. Solving the vehicle routing problem with adaptive memory programming methodology, Computers and Operations Research, Amsterdam, v.32, p. 2309–2327, 2005. TOTH, P.; VIGO, D. The granular tabu search and its application to the vehicle routing problem. Informs Journal on Computing, Catonsville, v.15, p.333- 348, 2003. TÓTH, S.F.; ETTL, G.J.; KÖNNYU, N.; RABOTYAGOV, S.S.; ROGERS, L.W.; COMNICK, J.M. ECOSEL: Multi-objective optimization to sell forest ecosystem services. Forest Policy and Economics, Amsterdam, v.35, p.73-82, 2013. VAIRA, G.; KURASOVA, O. Genetic algorithm and VRP: the behaviour of a crossover operator. Baltic Journal Modern Computing, Latvia, v.1, n.3-4, p.161-185, 2013. VICTORINO, I.; R.; S. Otimização de um reator industrial de produção de álcool cíclico utilizando algoritmos genéticos. Tese de Doutorado. Universidade Estadual de Campinas, Campinas, 2005. WAGNER, H.; H. Rethinking the linear regression model for spatial ecological data. Ecology, Amsterdam, v.94, n.11, p.2381-2391, 2013. WANG, S.; LU, Z.; WEI, L.; JI, G.; YANG, J. Fitness-scaling adaptative genetic algorithm with local search for solving vehicle routing problem. Simulations of Urban Transportation Systems, London, v.92, n.7, p.601-616, 2016. WANG, Y.; BAUERLE, W.; L.; REYNOLDS, R.; F. Predicting the growth of deciduous tree species in response to water stress: FVS-BGC model parametrization, application, and avaliation. Ecological Modelling, v.217, p.139-147, 2008. WASSAN, N.; A.; NAGY, G. Vehicle routing with deliveries and pick-up: Modelling issues and meta-heuristics solution approaches. International Journal of Transportation, Sandy Bay, v.2, n.1, p.95-110, 2014. WEI, L.; ZHANG, Z.; LIM, A. An adaptative variable neighborhood search for a heterogeneous fleet vehicle routing problem with three-dimensional loading 38 constraints. IEEE Computational Intellingence Society, Osaka, v.9, n. 4, p.18- 30, 2014. WEINTRAUB, A.; EPSTEIN, R.; MORALES, R.; SERON, J.; TRAVERSO, P. Truck scheduling system improves efficiency in the forest industries. Institute of Operational Research and the Management Science, Maryland, v.26, n.4, p.1-12, 1996. YAO, L., SETHARES, W.A. Nonlinear parameter estimation via the genetic algorithm. IEEE Transaction on Signal Processing, v.42, n.4, p.927-935, 1994. YOUSEFIKHOSHBAKHT, M.; DIDEHVAR, F.; RAHMATI, F.Solving the heterogeneous fixed fleet open vehicle routing problem by a combined metaheuristic algorithm. International Journal of Production Research, London, v.52, n.9, p.2565-2575, 2014. ZAINASHEV, T. Using a genetic algorithm in the optimal investment portfolio decision-making process. Proceedings of the International Multiconference of Engineers and Computer Scientists, Hong Kong, v.1, 2011. 39 SEGUNDA PARTE - ARTIGOS 40 ARTIGO 1 - A hybrid method for fitting nonlinear regression models Cássi Augusto Ussi Monti, Lucas Rezende Gomide, Fausto Weimar Acerbi Júnior e Antônio Carlos Ferraz Filho Artigo redigido conforme normas da Revista Forest Ecology and Management (Submissão em 22/02/2018), sendo versão sujeita a alterações. 41 ABSTRACT Regression analysis is widely applied in many fields of science and engineering to predict any hard to determine variable. In general, nonlinear regression models are more complex than linear regression. The former models request initial parameters to be fitted and the procedure to estimate these values is not simple. In this work, we propose a new hybrid method to fit nonlinear regression models using height diameter models as an example. The proposed hybrid method with dynamic range strategy has two stages; in the first stage we adopted the genetic algorithm (GA) to automated prediction of the initial parameters. The second stage consisted of applying the Gauss-Newton (GAGN), NL2SOL (GAPORT) or Levenberg- Marquardt (GALM) algorithms to get the fitting. We tested these three algorithms combined with the genetic algorithm. Twelve well-known models were chosen from literature and 12 databases were selected to test the hybrid method performance. According to the results, the GA was able to predict the initial parameters, helping the traditional algorithms to converge perfectly in most cases. However, the success of fitting rate was reduced for models with four parameters. We found that the Levenberg- Marquardt algorithm was superior to Gauss-Newton and NL2SOL for all database and models. The tests showed that using only GA to predict the final parameters was also robust and similar to the hybrid method. The classical 42 methods (stage 2) had no statistical difference to GA (stage 1). We conclude that hybrid method is effective to estimate the initial parameters. Keywords: nonlinear parameterization; dynamic range; genetic algorithm; fitting curves; hypsometric relation. 1 INTRODUCTION In nonlinear regression models, the relation between the predictor and response variables follows the general form Y=f(t,β), where f is a nonlinear function with at least one parameter (Fischer and Schönfelder, 2017). The main advantages of these models are associated to the biological interpretation identified in some examples of nonlinear models. In biological behavior models, generally β0 is defined as the asymptote (maximum potential of response variable), β1 is the biological constant, and the growth rate in the maximum biological potential is β2 (Fekedulegn et al., 1999). The nonlinear regression models disadvantage is associated to requirement of initial parameters for the algorithm’s solution to converge. In Forestry, some nonlinear models is often applied to predict the stand variable patterns, e.g. total height (Huang et al., 1992; Fang and Bailey, 1998; Mehtätalo et al. 2015), tree crown (Fu et al., 2013; Fu et al., 2017), dominant height (Rayner, 1992; Fontes et al., 2003; Lekwardi et al., 2012), growth and yield (Shvets and Zeide, 1996; Mønness, 2011; Fischer and Schönfelder, 2017). 43 However, these functional relations only work with high correlation between variables. In general, linear model fitting is obtained by optimal methods, i.e. ordinary least squares regression (OLS), which is used to predict the model’s parameters. On the other hand, nonlinear regression fitting method is always approximated (Chumney and Simpson, 2006) and not optimal. The OLS was firstly exposed by Gauss (1809), which launched its principles to minimize the residual sum of squares. The residuals are derived from the differences between observed (yi) and estimated ( ) values. The estimated parameter vector is the least squares estimator of the response variable (Sprott, 1978). Moreover, errors must be random, independent and must follow normal distribution with zero mean and constant σ². These assumptions are the beginning of classical statistic’s theory proposed by Fisher (1925) and Snedecor (1937). Under these circumstances, the OLS has changed according to Gauss- Markov theorem, where the least squares estimator from some general regression model could provide a minimum variance and an unbiased linear estimator (Lovric, 2011). It is the ordinary least squares estimator of this model, i.e. provides an optimum error minimization. Advantages and disadvantages are observed in those procedures due to loss of information (Fekedulegn et al., 1999; Özçelik et al., 2013; Fu et al., 2013; Fu et al., 2017a; Fu et al., 2017b). 44 Algorithms such as Gauss-Newton and Levenberg-Marquardt were developed to predict the nonlinear regression`s parameters (Draper and Smith, 1981). The first algorithm assumes a modification of Newton’s method incorporating linear search. This strategy provides savings in function and derivative evaluation time, which is the least squares problem resulting on small residuals at the solution and fast local convergence (Nocedal and Wright, 1999). The authors also point out Levenberg-Marquardt method as promissory due to avoidance of the Gauss-Newton method weaknesses, the Jacobian’s rank- deficient shape. However, when the Hessian matrix components are ignored in the local convergence phases, these two methods are similar. The initial parameters values are needed to initiate the algorithm convergence and provide an accurate model fit. An inconvenience about using nonlinear regression models (Tomioka et al., 2007). This disadvantage may increase if the initial values provided do not make sense or a biological pattern is not achieved. Commonly, the fitting process returns non-convergence outcomes. However, the partial derivative is another flexible option to obtain the initial parameters and this approach usually results in excellent outputs (Fekedulegn et al., 1999). Additionally the parameter estimation is not user friendly as in linear regression and requires specialists to guarantee the convergence (Ratkowsky, 1983). 45 Recent research has focused on solving the non-convergence situation by applying artificial intelligence based on robust search algorithms or recognition pattern approaches. This has been attributed to a range of factors, including complexity of variables and infeasible solution in operational time. In fact, optimization problems need nonconventional methods to be solved, as genetic algorithm for instance (Yao and Sethares, 1994, Chatterjee et al., 1996; Karr et al., 1998; Kapanoglu et al., 2007; Chen, 2015; and Jin et al., 2016). However, problems arise when we desire to model some functional and biological pattern using nonlinear regression model without knowing the models initial parameters values. This problem may not be solved in a straightforward manner by using traditional algorithms or guessing them in some cases. Considering these issues, this research relies on proposing a hybrid algorithm method to predict the initial parameters for nonlinear regression models applied in forestry. 2 MATERIAL AND METHODS 2.1 Data base structure The database used to test the proposed model fitting the algorithm was designed to encompass several forest types. The idea is to perform a complete and robust test using the hybrid algorithm with dynamic range. We used 5,834 46 individuals trees observations from 12 different locations, where the diameter at breast height (d) and total height (h) were measured (Table 1). The areas do not have any geographical correlation or similar structure, increasing the sites diversities and patterns. We also took care to represent a great level of landscapes, soil, vegetation types, hydrologic regimes and climate. In order to test the fitting conditions, we selected some types of vegetation classified by Veloso et al. (1991) as Brazilian Savannah (Cerradão), Rainforest (Mata Atlântica) and Semi-Deciduous Seasonal Forest (Floresta Estacional Semi- decidual). We also used the database from the commercial species (Hardwood and Softwood) most used in Brazilian plantations, represented by many growth stages and silvicultural treatments, as follows: Eucalyptus spp (2x2m; 2, 6, 8 and 40 years) and Pinus spp (2.2x2.2m; 5, 7 and 8 years – 1 selective thinning at 8 years). Finally, we selected two Brazilian’s native species which has a monodominance behavior (Araucaria augustifolia (Bertol.) Kuntze and Xylopia brasiliensis Spreng.), important for ecological, cultural and social aspects. 47 Database N Average Max Min SD d h d h d h d H SA* 100 19.7 12.3 43.0 17.1 7.0 7.1 80.5 22.7 NMT** 1,975 9.6 6.5 99.0 30.0 2.9 1.4 54.8 30.2 NX** 73 24.3 15.1 50.6 21.5 6.1 7.5 126.4 33.8 NSD** 1,975 11.7 9.7 77.4 25.0 3.0 2.1 82.5 36.0 NC** 609 11.1 7.1 56.8 20.6 3.1 1.8 71.6 28.4 HE2* 229 11.7 15.1 14.9 40.3 4.8 10.0 15.6 22.7 HE6* 357 16.9 24.1 24.3 27.9 6.3 7.0 30.6 29.0 HE8* 28 19.1 24.4 29.1 32.5 6.5 9.8 63.6 56.0 HE40* 188 17.1 25.1 38.2 50.0 2.7 4.0 84.2 121.6 SP5* 100 11.7 9.6 20.4 12.7 4.4 5.1 31.0 15.5 SP7* 100 14.7 13.7 24.1 17.8 6.3 7.9 35.1 16.9 SP8* 100 14.4 14.8 23.1 18.8 7.1 10.4 30.3 15.9 Table 1 Summary of descriptive statistics from the study’s database involving even-aged (*) and uneven-aged (**) forest types. (Notes: N – number of trees measured; d – diameter at breast height (cm); h – total tree height (m); SA – Araucaria augustifolia; NMT – Atlantic forest trees (Mata Atlântica); HE8 – 8 years-old Eucalyptus spp; HE40 – 40 years-old Eucalyptus spp; NX – Xylopia brasiliensis trees; SP7 – 7 years-old7 Eucalyptus spp; SP5 – 5 years-old 5 Eucalyptus spp; SP81S – 8 years-old after thinning Pinus spp; HE6 - 6 years-old Eucalyptus spp; HE2 – 2 years-old Eucalyptus spp; NSD – Deciduous forest trees (Semi-Decidual); NC –savanna (Cerradão; Average – average of observed h-d data; Max – maximum value of d-h data; Min – minimum value of d-h data; SD – standard deviation of d-h data). 2.2 Fitting approach 2.2.1 Nonlinear regression models Twelve nonlinear biological growth models were chosen aiming to assess the hybrid method robustness. These type of models have been frequently applied to predict variables in forest science, e.g. total tree height (Temesgen et al., 2008; Özçelik et al., 2013; Misik et al. 2016), dominant tree height (Wang et 48 al., 2007), individual tree volume (Huang and Titus, 1999) and biomass (Parresol, 2001). Prediction of the total height (h) of a plot is widely applied in forestry studies, usually presenting weak correlation with diameter at breast height (d), mostly measured variable. Tree height estimation or hypsometric relation (h-d) has financial importance to reduce the inventory costs (Ribeiro et al., 2016), since fewer trees have to be measured in the field. The hypsometric patterns tend to have nonlinear behavior, while the linear h-d equations may only be applicable for the tree sizes tested and standard conditions (Bi et al., 2012). Due to the large applicability in forest management, we used the h-d relation to test our algorithm. Twelve nonlinear regression models were selected from the forest science literature (Table 2). 49 Parameters Model Formula ID Reference 2 1 Negative Exponential1 ( ) 0 (1 exp( ))βd dh      2 Burkhart2   1 0 β β ex= p ε+ d h d        3 3 Monomolecular1     0 1 2β 1 - β exp( -β ε= ) + d h d 4 Mitcherlich1 ( ) 0 1 2β β β d dh     5 Gompertz1 ( ) 0 1 2 exp( exp( ))β β βd dh        6 Logistic3 ( ) 0 1 2 /1 exp(( ) / )β β βd dh     7 Chapman- Richards4 2β ( ) 11 exp( )(1 )ββd dh     8 Bailey2 2β ( ) 0 1 (1 exp(- ))β βdh d     4 9 Von Bertalanffy1 3 3 1 1 β 1 β 2( ) 10 ( exp( β ))ββd dh        10 Bailey2 3β ( ) 0 1 2β (1 β exp( β ))dh d     11 Zeide2 3β ( ) 0 2β exp( β1 exp( β ))dh d       12 Richards2 3β ( ) 0 1 2β (1 β exp( β )dh d     Table 2 Nonlinear regression models used to test the hybrid algorithm with dynamic range (Note: h(d) – total tree height as function of diameter at breast height (d); β0 – asymptote parameter; β1 – the curve slope parameter; β2 – the rate to reach the asymptote parameter; β3 – allometric constant parameter; ε – random error; e – mathematical exponential expression.1 – Fekedulegn et al. (1999); 2 – Fang and Bailey (1998); 3 – Calegario (2002); 4 – Retslaff et al. (2015)). 2.2.2 Hybrid method In this section we describe the hybrid method to fit any nonlinear regression models without knowledge of the initial parameter values. The method has two stages: a) Genetic-dynamic approach: the genetic algorithm coded to predict the initial parameters values, and b) Statistical approach: the 50 model fitting applying the classical statistical nonlinear regression algorithms (Gauss-Newton, Levenberg-Marquardt or NL2SOL). We have been renamed it as GAGN (Genetic algorithm + Gauss-Newton), GALM (Genetic algorithm + Levenberg-Marquardt) and GAPORT (Genetic algorithm + NL2SOL). 2.2.2.1 Stage 1: Genetic dynamic approach The genetic algorithm (GA) implementation was set up after initial parameterization tests. In fact, every model tested estimates the total tree height using the predicted parameter of each individual from the population. According to this procedure we defined the population size, fitness function and genetic operators (selection, crossover and mutation). The population size is represented by a group of individuals to be evaluated for every generation. In this work, we have found a maximum population size of 300 individuals and 100 generations as stop criteria. Steady state strategy was adopted to define the population from the generation n+1. The evaluation function or fitness f strategy adopted was based on the residual sum of squares in Equation 1, where hi is the observed total tree height and i ∈ {1, …, N}; and is the estimated total tree height value from the tested nonlinear regression model.   2N i – i i hf   h (1) 51 Directly connected to the fitness function, the parameters range was defined previously as initial genetic algorithm set. The parameters assume a random generation structure within the symmetric range interval [-r, +r] for all parameters. Initially, this strategy assigns r with high values (500) to define the parameter’s interval search by the GA. The arbitrary number was chosen considering the database used and the model’s complexities. Therefore, the dynamic range adaptation was requested to accelerate the algorithm convergence to better local solution bounds. Thus, the range interval changes dynamically and the algorithm explores promissory solution areas. The r values change according to the rule r = min {βp n+1}, where βn+1 ∈ {1,…, P} and to the best fitness reached in any generation. Previously tested, the dynamic range function (dyrange) substantially improved the success of algorithm convergence rate (Fig. 1). 52 1: Begin function: 2: dyrange ( f n , f n+1 , βn+1){ 3: if ( f n > f n+1): 4: r = min {βp n+1}  βn+1 ∈ {1,…, P} 5: else: r = r 6: return range [+r, -r] 7: } 8: End function Fig. 1. Procedure of the dynamic range used to guide the random process of genetic algorithm. Complementary, the genetic operators imitates the theory of evolution by means of natural selection with traits that enhance survival and reproduction. This procedure is the key of the algorithm and randomly simulated n evaluative generation by genetic operators. 1) Selection Operator: this operator plays a most adapted individual selector, similar to the Darwin’s Theory of biological selection that states that the most adapted individual has the greater probability to prevail in nature. The computational approach is implemented by a random search from the population, creating a new set and selecting the parents according to their fitness values (Kapanoglu et al. 2007). There are several ways to implement the selection operator, for more 53 information see Tomioka et al. (2007), Hadi and Gonzalez-Andujar (2009), Mitra et al. (2011) and Karami and Wiens (2014). The tournament selection (Equation 2) was chosen to control the diversity losses (Jin et al. 2016), selecting two individuals (f(x1) and f(x2)) from the parents pool (population). In order to obtain the most adapted parents, this operator chooses the best value of fitness in every pairwise defined states. Finally, the best individuals are selected to crossover. Selection (2) 2) Crossover Operator: is the main operator of the genetic algorithm, which provides exploitation phases of the solution search (VanderNoot and Abrahams, 1998). The authors also explain about the effect of crossover in the solution qualities by swapping genes between parents. Considering these aspects, we previously tested the exchange rates of the crossover operator. Due to the reduced number of model parameters (2, 3 or 4), one gene swapping from the parent’s chromosomes and creation of two offspring was adopted (Figure 2), as observed in VanderNoot and Abrahams (1998), Yoshimoto et al. (2003), Sevaux and Mineur (2007) and Tomioka et al. (2007). 54 Fig 2. Crossover operator structure applied to change a random locus at a random individual scheme. 3) Mutation Operator: according to VanderNoot and Abrahams (1998), mutation performs a solution exploration increasing the diversity of the population. This operator imitates a biological mutation. Darwins’s Theory says there are some “random changes” in an individual and, whether these changes may provide skills increasing, than it will be kept from parents to children, maintaining differences from the others individuals (diversity). This statement is about chromosome mutation, providing some gene changes into living being’s characteristics. Artificially, the randomly selected gene from the chromosome is replaced within another parameter range. The mutation maintains genetic population diversity and provides an escape solution from some local optima (Yoshimoto et al., 2003; Roush and Branton, 2005; Sevaux and Mineur, 2007; Tomioka et al., 2007). We applied one locus change per generation and 0.6 random mutation rates into the population and k offspring set k ∈ {1, …, π}, and π works as an offspring size vector. 55 2.2.2.2 Stage 2: Statistical approach The second stage of the hybrid algorithm is associated with the nonlinear regression model fitting. The initial parameters derived from the best fitness output from the GA convergence are used in this stage. We tested three algorithms largely applied in statistical software with iterative approach for nonlinear least square regression: Gauss-Newton (GN), Levenberg-Marquardt (LM) and NL2SOL (PORT). Consider the data set (xi, yi) to be fitted, where xi is the independent and yi is the dependent variable, i ∈ {1, …, N}, βj the parameters, j ∈ {1, …, M}, f is a nonlinear function of the predictor variables and the parameters and ε is the related random errors for each observation i in Equation 3. yi = f(xi , βj) + εi (3) The GN method uses a linear approximation to iteratively improve an initial parameters vector guess β0 for β, where β is a vector of parameters that minimize the residual sum of squares, and keep improving the estimates until there is no change (Bates and Watts, 1988). The authors stated that this approximation follows the first order Taylor series (Equation 4), which introduces the first order derivative matrix, i.e. Jacobian matrix J, in order of 56 each parameter [Jij = ∂f (xi, β)/∂βj]; and the component that gives the improvement iteratively, the Gauss increment [δ = (β – β0)]. f (xi , β) ≈ f (xi , β0 j)+Jn1(β1-β0 1)+Jn2(β2-β0 2)+…+Jnm(βj-β0 j) (4) The approximation can be re-written as in Equation (5). ɳ(β) ≈ ɳ(β0)+Jnm δ (5) Where ɳ(β) is the expectation surface into the response space. The general residuals approximation Z(β) can be expressed as a function of observation vector Y in the first iteration as in Equation (6). Z0(β) ≈ Y–[ɳ(β0)+J0 nm δ0] = Z0-J0 nmδ0 (6) According to Bates and Watts (1988), the Gauss increment may increase the sum of squares when the δn overflows to the region ɳ(β), so a step factor λ in the direction δ0 may provide a decrease in the sum of squares by a parameter modification in each iteration [β1=β0+λδ0], where λ is chosen to ensure that a solution n+1 is better than n [S(β1) < S(β0)]. The GN method then calculates a new Z1, new J1 and new δ1, and so on, until the increments do not change the current parameters of vector β. For more details about Gauss-Newton algorithm implementation see Bates and Watts (1988). The Gauss-Newton method was largely used as base to implementing others algorithms for nonlinear least squares problems, LM and PORT, for 57 example. The LM was idealized by Levenberg (1944) considering the Newton’s method into its update function and improved by Marquardt (1963), who incorporate the estimated local curvature information into the update function. The Levenberg-Marquardt algorithm was implemented as a robust technique by Moré (1978). The author approach the step bound Δ updating by its choice depending on the ratio ρ(β) between the actual reduction and the prediction reduction. This ratio is obtained by the decomposition on the linear system through measuring the agreement between the linear and the nonlinear function (Moré, 1983). This step bound concept is a prior attempt of trust-region approach, using ρ(β) to choose Δ given by the expression as in Equation 7. ρ(β)= (7) To update Δ, Equation 7 must be kept at a reasonable value, such that if ρ(β) ≥ 3/4, than Δ increases, but if ρ(β) ≤ 1/4, than Δ must be decreased. For more specific rules for updating Δ, see Moré (1983). The author still stated that α > 0 is the Levenberg-Marquardt parameter if | f (α) | ≤ σΔ, where σ ∈ (0,1) specifies the relative error in || Dβ(α) ||, where D is a diagonal matrix which takes into account the scaling of problem, i.e. the data base module (Moré, 1983). 58 The last tested algorithm by hybridization within GA was NL2SOL (PORT). According to Dennis (1977), the PORT algorithm is more reliable than the Gauss-Newton or Levenberg-Marquardt methods in large residual cases, even in the first algorithm versions. Dennis, Gay and Welsch (1981) gave the last NL2SOL (PORT) version used in recent packages and software. This algorithm is called by the authors as an adaptive nonlinear least squares because it uses some procedures found in GN and LM, such Gauss increment δ and trust- region concept that was truly implemented in LM by Moré and Sorensen (1983). This method has a strong similarity with GN if the Hessian approximation in PORT is not accounted for at the beginning, than the GN algorithm is clearly used (Dennis, Gay and Welsch, 1981). However, authors claims that PORT algorithm uses a double dogleg optimization (DDO) applied to find the Gauss increment. The DDO has a modification based on trust-region concept by Dennis and Mei (1975). According to the PORT developers, the important thing is the idea of having at βn a local quadratic model q of f and an estimate of a region in which q is trusted to represent f. The next point βn+1 is chosen to approximately minimize the q in this region, in order to use the improved information about f (βn+1) to update the size and shape of the trust- region. 59 In order to demonstrate the main procedure of hybrid algorithms we show its pseudo-code. It covers both stages, the GA procedures and the classic regression methods (Fig. 3). Hybrid method with dynamic range algorithm 1: Begin 2: stage 1: GA; 3: Initial parameters step: 4: population size L = 300 + offspring vector size k = 140 5: range r n ∈ [500, -500] n ∈ {1,…,P} 6: initial solution - f ( r n, h(d)) 7: Initialize iterative procedure: 8: While generation ≤ 100, Do: 9: Selection function 10: Crossover function 11: Mutation function 12: Dynamic range function 13: Loop 14: return βn= r n, where r n from the best solution f* 15: stage 2: adjustment phase; 16: GAGN as a function of βn 17: GAPORT as a function of βn 18: GALM as a function of βn 19: End Fig. 3. Procedure of the hybrid algorithm with all three fitting approaches. The hybrid algorithm implementation was developed in R software (R DEVELOPMENT CORE TEAM, 2016). In the first stage we coded a real genetic algorithm (GA) and in the second stage we used two functions to fit the nonlinear regression models. These functions were nls() for the Gauss-Newton and NL2SOL algorithms, developed on Bates and Watts (1988) and Dennis et al. (1981), respectively; and nls.lm() for the Levenberg-Marquardt algorithm, 60 developed by Moré (1978). All the computation tests were performance using Desktop/Intel® Core™ i5-3570 (3.4GHz). 2.3 Experiment analysis We selected 12 most applied nonlinear biological models in literature with highly nonlinearity features and varying number of coefficients (2, 3 and 4 parameters). A combinatorial problem arises when every model is matched with 12 databases. Besides, this procedure runs 100 times to validate the results due to the stochastic algorithm behavior. At the end, the experiment repetitions reached 14,400 (12x12x100) runs. Finally, we processed the experiment 57,600 times, when combined the 3 nonlinear regression methods (GN, PORT and LM) and GA. Statistical criteria derived from the error were applied to compare the regression methods: a) MAPE - mean absolute percentage error (8); b) MSD - mean squared deviation (9), and c) MAD -mean absolute deviation (10). These statistics are often applied to compare the goodness-of-fit of different fitting methods. MAPE = x100 (8) 61 MAD = (9) MSD = (10) Complementary, the first stage results were also used to check the quality of the GA solutions. In fact, this strategy was designed to predict the GA’s robustness to estimate the parameters without the need of the regression algorithms. These tendencies and patterns were compared using the disagreement test (Borba and Nakano, 2016). This test fits a simple linear regression model with two variables, in this case the variables come from the GA and hybrid methods. If the intercept is 0 and angular coefficient statistically equal to 1, than both methods are equal. However, this test does not control type I errors, so we used the Kappa, Pearson and Kendall coefficient of concordance to control it and accurately compare the height estimates between the different outputs. 3 RESULTS 3.1 Statistical regression models and database pattern 62 The databases have two groups of stands (plantation and native) and the number of successful fittings varied according to the model. We identified the higher h-d correlation of the plantations from HE6 (78.09%), and the lower to HE2 (48.84%), and 74.64% in average. Considering the natural forests, the best database was NX (74.76%) and the worse was NC (73.76%), with an average of 72.68% of correlation coefficient. The negative effect of low h-d correlation was the main source of a non- convergence of the algorithms for some databases. The Burkhart (model 2), Logistic (model 6) and Bailey (model 10) had the highest successful fitting rates, in each respective parameter group (2pg, 3pg and 4pg). Otherwise, the worse fitting models were the Negative Exponential (model 1), Mitcherlich (model 4) and Richards (model 12). However, the hybrid method is capable of predicting the initial parameters without previous analysis about the models and the database. The nonlinearity feature of these models set up some intrinsic frames which determine the model’s adherence with the database and consequent fitting. 3.2 Hybrid method 3.2.1 Stage 1: Genetic algorithm 63 The coded GA was powerful to predict the initial parameters for all tested models and databases, where the fitness function decreased over the iterations, as expected (Fig. 4). The GA improves the initial population using a randomized procedure. This has been attributed to a range of factors, including mutation rate, population size and also the dynamic range strategy. However the 4pg group was the hardest to minimize the fitness function. Fig 4. Best fitness average (Left panel) and Dynamic range average (Right panel) of parameters range behavior, for each parameter class. The good local optimum is usually hard to detect and requires high efforts for any meta-heuristic. Due to this problem, the dynamic range strategy was essential for the GA to converge to better solutions. The values range 64 interval changes when the algorithm finds the better solutions and this helps the stochastic pattern to improve the population fitness average. The robust convergence algorithm combined with short run-time is desired to find the initial parameters. This problem can be solved in a straightforward manner using GA. 3.2.2 Stage 2: Statistical fitting methods The minimal attempt proposed to fit the nonlinear regression models, i.e. the minimal times that the hybrid algorithm must be processed until it gets a fitting, was calculated (Figure 5). Based on the results, when the model complexity increases, the hybrid algorithm (stage 2) fails to predict the parameters until 23 proceedings from 100 attempts for 4pg models. However, the convergence depends on the results of the GA (stage 1). Considering only GALM, the fitting success was observed for the first attempt. 65 Hybrid Methods GALM GAGN GAPORT Fig 5. Requested attempts until the first fitting for each hybrid method tested and applying the best identified model (Bailey) for 4pg. (Note: * No fitting was achieved). According to the results, the 4pg models requested more attempts to reach the first successful fit. Considering only the GALM, the combination of algorithms was efficient at the first attempt in almost all tests. However, in some situations the initial predicted parameters were questionable, considering estimation quality and biological patterns. Comparing GAGN and GAPORT within the same initial parameters, the last method was superior in success convergences, although requesting more efforts until the first fitting. The same tendency was observed for 2pg and 3pg. There were two databases that yielded with no fitting by GAPORT (NC; SA) and six databases (NC; SP7; SP5; NX; NMT; SA) for GAGN. The fitting success rate decreased with increasing model parameter number for all methods. Therefore, the hybridization strategy was 66 responsible to predict the initial seed quality, i.e. the initial parameters with quality without knowing previous information about the data base nor models (Table 3). Number of parameters Hybrid Methods Sucess Average Rate (%) 2 GALM 100.0 GAGN 97.9 GAPORT 98.4 3 GALM 100.0 GAGN 60.5 GAPORT 66.5 4 GALM 100.0 GAGN 3.7 GAPORT 3.8 Table 3 The fitting success average rate of nonlinear models applying the hybrid method for all databases tested. 3.3 Genetic algorithm versus Hybrid method In order to show the hybrid method ability in estimating the initial parameters of the models, we showed only the worse data base analyzed, with the lowest h-d correlation and model fitting success rate, which was NC (Native Cerradão) composed by wooded savanna forest type. In general, the GA was consistent and robust as observed in hybrid method results. The worse database (NC) did not affect the GA performance according to residual statistics (Table 4). However, the GALM was better than 67 GA only for 4pg. On the other hand, the GAGN and GAPORT failed to fit the models in this group. In fact, residual analysis shows a high similarity between all methods. Therefore, the solutions from GA were always inferior with high correlation of all methods tested as expected. Number of parameters Statistics Hybrid Methods GA GALM GAGN GAPORT 2 MAPE (%) 25.6 25.6 25.6 25.6 MAD (m) 1.5 1.5 1.5 1.5 MSD (m²) 4.1 4.1 4.1 4.1 3 MAPE (%) 24.6 24.7 24.66 26.3 MAD (m) 1.4 1.4 1.42 1.5 MSD (m²) 3.7 3.7 3.67 4.4 4 MAPE (%) 24.9 NA NA 25.1 MAD (m) 1.4 NA NA 1.4 MSD (m²) 3.7 NA NA 3.8 Table 4 Average of residual statistics for genetic algorithm and hybrid methods considering all tested nonlinear regression models (Note: MAPE – Mean Absolute Percentage Error; MAD – mean absolute deviation; MSD – mean of standard deviation; Mean – average of high estimated h value; NA – No Answer). We applied the disagreement test and coefficient correlation to evaluate the similarity between the hybrid methods and genetic algorithm. Due to the GALM superiority among the others hybrid methods, we ran the test only between GA and GALM. The results showed similar β1 values for 2pg and 3pg, which means GA and hybrid methods, represented by GALM, do not disagree with 2pg (GA = 0.9944 and GALM = 1.0055) and 3pg models (GA = 0.9937 and GALM = 1.0062). The opposite pattern of β1 was observed for 4pg models 68 (GA= 0.8692 and GALM = 1.1373), demonstrating more disparity between the fitting strategies. Figure 6 shows the behavior of β1 for all model groups. a b c Fig 6. Disagreement test between GA and GALM for the NC database. (Note: a- Burkhart (2pg); b- Logistic (3pg); c- Bailey (4pg)). An additional analysis was done to guarantee the control of type II statistical error. We applied the Kendall and Pearson correlation and Cohen’s Kappa tests to reinforce the reliability of the estimated values from the GA for fitting curves. After the comparisons, we affirmed that the GA-GALM values agree, according to the Kappa test (1.0), Pearson (1.0) and Kendall (0.99) correlation for the best models of each parameters group. Therefore, the fitting methods are considered similar to each other, thus, GA may be used to predict the values of non-observed data with the same reliability found by classical statistical methods (stage 2). 69 4 DISCUSSION AND CONCLUSION We highlight some issues that have been previously raised concerning nonlinear regression models. These issues are usually related to the database, mathematical model shape, initial parameters value and fitting method. Fitting curves with high quality is not an easy task. In forest science, h-d correlation seems decrease over time demonstrating a degenerated pattern of hypsometric relation. In our work were identified a low h-d correlation due to the forest management regimes, age and native forest types influencing this pattern. Therefore, the fitting curve process becomes more complex and sometimes infeasible. Our results may reinforce such statement, suggesting an increasing effort for models with high number of parameters to fit properly. In fact, modeling process involves many aspects to guarantee unbiased predictions. This leads to new methods of artificial intelligence to solve these problems. Traditional algorithms with low success fitting rate have been declined in many statistical packaged software and frame languages. Therefore, there are knowledge gaps to be explored in regression analysis. Many of unsuccessful fittings could be associated to the method and our hybrid algorithm could help them to get a success fitting on negative exponential function (Fekedulegn et al., 1999) and Weibull (Mahanta and Borah, 2014) for instance. In correlated works raised the same application of artificial intelligence 70 as fitting curves strategy (Křivý et al., 2000; Tvrdík et al., 2007), bringing other scenarios to apply our hybrid method in the future. The found standard statistical packages showed more fail rates than intelligence computational approach. In some cases of our study we observed a similar tendency, but more associated to the GAGN and GAPORT than GALM or GA. Regardless of the database tested, the negative exponential functions were fitted with GALM and also GA. A notable improvement in residual statistics was not detected between the last two methods, which indicates their equality. Altunkaynak and Esin (2004) previously defined the initial parameter range applying the binary codification for each individual of the GA population. The authors found a similarity between standard deviation for Gauss-Newton and Genetic Algorithm in several tested models (Gompertz, Richards and Logistic), as we observed in our study. These results do not surprise due to several studies showing the GA potential to predict complex functions (Altunkaynak and Esin, 2004; Roush and Branton, 2005; Sevaux and Mineur, 2007; Dai et al., 2009; Zúñiga et al., 2014). In order to illustrate the positive effect of the hybrid methods, we identified initial parameters for all databases and models tested. The stochastic outcomes under uncertainties from GA are robust and efficient to predict the parameters as observed in Kapanoglu et al. (2007), Chen (2015) and Jin et al. 71 (2016). Our hybrid method may provide a framework for adapting in any statistical packaged software, mainly if the dynamic range function is applied. Tomioka et al. (2007) proposed a method to change the parameters domain using several genetic algorithms, called ADM, similar to our approach. They changed the likelihood function to search a new parameters domain (range) from the GA. A comparable strategy was applied in our algorithm, where the advantage of employing dynamic range makes th