Viajar é a opção de lazer de muitas pessoas, e hoje os aplicativos que calculam rotas e nos ajudam a chegar nos lugares utilizando GPS são bastante úteis. Eles resolvem um dos problemas clássicos da otimização: o problema do caixeiro viajante (Traveling Salesman Problem). É bem provável que este seja o problema mais estudado na área de otimização discreta.

Índice

Por que ele é importante?

Imagine que você seja um representante comercial e precise visitar clientes em 5 cidades diferentes. Cada cidade deve ser visitada apenas uma vez e, ao final, você volta para a primeira. Você precisa escolher a menor rota possível entre elas. Quantas opções de rotas existem?

Você precisa escolher uma das cinco cidades para começar, depois uma das quatro para ser a próxima, aí ficam mais três… e assim vai. Matematicamente, você tem 5 x 4 x 3 x 2 = 120 opções. Neste caso pequeno seria possível pedir que o computador calculasse a distância percorrida em todas as combinações e retornasse a menor.

Mas e quando temos mais cidades? Se falarmos apenas das capitais brasileiras, temos 27 cidades (26 estados e um Distrito Federal). A quantidade de combinações é um número de 28 zeros.

Isso significa que precisamos de métodos que calculem estas rotas de maneira mais inteligente, sem precisarmos verificar todas as combinações.

Eu andei pesquisando na internet, mas não encontrei algum lugar em que houvesse o cálculo de uma rota entre todos os municípios brasileiros. Então a tarefa aqui será aproximar a menor rota entre mais de 5.000 municípios brasileiros.

O número de municípios brasileiros não é fixo, mas a melhor lista que consegui, com as posições geográficas de cada um, contém 5.509 municípios. Ela foi retirada originalmente deste link.

Existem alguns métodos para calcular esta rota, sem eles seria impossível calcular, por força bruta, qual é a melhor, mesmo que tivéssemos todos os computadores do mundo à nossa disposição.

Os arquivos utilizados neste artigo estão disponíveis neste repositório GitHub. O script usa Python puro, então recomendo rodar o script usando PyPy para ir mais rápido.

Dados

Para cada município temos o nome, estado, e a latitude e longitude. A distância entre eles será calculada baseada nessas informações. Apesar de haver diferença com a realidade, já que as estradas de uma cidade a outra não formam, necessariamente, uma linha reta, esta será uma boa aproximação.

Dificilmente chegaremos ao mínimo absoluto com o método demonstrado aqui, mas ele é uma maneira de se aproximar bastante deste número. Para efeitos práticos, aproximar-se do mínimo absoluto sem precisar de um grande poder computacional, e em pouco tempo, é o bastante.

Preparando os dados

Com o XLS transformado em CSV, vamos criar funções dentro do script Python para processar os dados, colocando-os em estruturas adequadas para rodarmos os algoritmos.

Este script é baseado no código disponibilizado no curso Discrete Optimization, e na implementação que utilizei para resolver uma das tarefas, que são alguns TSPs.

Primeiro criamos uma estrutura de dados para armazenar as informações de cada cidade, usaremos a named tuple. Esta estrutura cria uma subclasse da tuple original cujos dados podem ser acessados usando nomes, em vez dos índices da tuple padrão.

Agora criamos uma função para carregar o arquivo CSV e armazenar cada linha numa named tuple contendo o ID, latitude, longitude, nome e estado onde fica a cidade. Estas serão armazenadas numa lista. Retornamos a lista das cidades, e a quantidade.

town = namedtuple("town", ['id','name','lat', 'long'])

def parseInput():
    
    with open('MunicipiosBrasil.csv','r') as input_data_file:
        town_list = [l.split(';')[:4] for l in input_data_file.read().split('\r')[1:]]
        
        town_list = [town(int(t[0])-1,t[3] ,math.radians(float(t[1].replace(',','.'))), math.radians(float(t[2].replace(',','.')))) for t in town_list]
        
        townCount = int(len(town_list))

        return town_list, townCount

Note que, decidi subtrair 1 de cada ID para termos um índice começando em 0, e além de formatar a latitude e longitude de maneira que o Python pudesse transformá-las em floats, já aproveitei para transformá-las em radianos. Originalmente elas estão em graus decimais. Isso facilitará os cálculos de distância.

Resolvendo o problema

Mixed Integer Program

Uma maneira de obter uma solução bastante próxima da melhor é usar um software que resolva MIPs. Este softwares resolvem problemas de otimização discreta aproximando-os com Programação Linear, de maneira contínua, e depois explorando árvores binárias cujos nós são soluções com algumas variáveis discretas fixas.

Existem várias formulações para resolver este problema, e é possível saber quão próximos estamos do ponto mínimo (o intervalo entre o ponto mínimo e a nossa solução) durante o processo de solução.

O melhor software gratuito para resolver MIPs que encontrei foi o SCIP, que pode ser utilizado em Python com o módulo Numberjack. Mas fique atento, dependendo da formulação, um problema deste tamanho exigirá bastante memória e tempo para chegar à solução.

Cálculo da distância

Para calcularmos a distância entre dois pontos originalmente identificados pela latitude e longitude decimais, vamos transformá-los em radianos e usar a Aproximação Equiretangular. Como as distâncias entre os pontos, em geral, são pequenas, não teremos problemas com erros muito grandes, e esta acaba sendo uma fórmula bastante simples.

def length(town1, town2):
    R = 6371.
    x = (town2.long - town1.long) * math.cos((town2.lat + town1.lat)/2.)
    y = town2.lat - town1.lat
    return math.sqrt(x**2 + y**2) * R

Você pode obter mais informações sobre esta fórmula, e outras maneiras de calcular esta distância neste link, que é a fonte que utilizei para decidir qual fórmula utilizar.

(Meta) Heurísticas

Greedy

def greedy(sol, townCount,town_list):

    for i in xrange(townCount-1):
        cost = length(town_list[sol[i]], town_list[sol[i+1]])
        for j in xrange(i+2,townCount):
            cost2 = length(town_list[sol[i]], town_list[sol[j]])

            if cost2 < cost:
                old = sol[i+1]
                sol[i+1] = sol[j]
                sol[j] = old
                cost = cost2

    return sol

A solução que normalmente nos vêm à cabeça primeiro é simplesmente selecionar a cidade com a menor distância como a próxima a ser visitada no percurso. Esta é uma solução rápida e em boa parte dos casos vai atingir um resultado satisfatório. Além disso, podemos utilizá-la como solução inicial para ajudar algoritmos mais complexos a chegar mais rápido a uma solução melhor.

No nosso caso, esta solução nos dá uma distância total de 156.760 km. Esta será a base de comparação das outras soluções.

Iterated Local Search 2-Opt

Outra maneira de aproximar a distância mínima é usar uma busca local. Este método, diferente do MIP, funciona melhor com problemas grandes, mas normalmente não nos dá a menor rota.

Basicamente consiste em encontrar soluções vizinhas à nossa atual. Por exemplo, trocando duas cidades de posição. Dentre as maneiras de encontrar soluções vizinhas, uma que funciona muito bem, e é bastante utilizada é a K-OPT.

Se visualizarmos o problema como uma rede (ou um grafo): cada cidade sendo um nó, conectado à cidade anterior e à próxima, e o valor desta conexão é a distância entre elas. Neste caso, K é o número de nós que vamos trocar de posição para ver se melhoramos a solução.

Imagem retirada da Wikipedia

Neste exemplo será utilizado o 2-OPT, então buscaremos vizinhos trocando 2 nós. Em código isso é feito simplesmente selecionado duas cidades, e invertendo a ordem das cidades entre elas.

def selectNextTour2Opt(solution, townCount):
    first = randint(0,townCount)
    last = randint(first,townCount)
    
    slice = solution[first:last+1]
    slice.reverse()
    new_solution = solution[:first] + slice + solution[last+1:]

    return new_solution

Num exemplo simples: imagine que as cidades A – B – C – D estão conectadas, e resolvemos aplicar o 2-OPT, selecionando as cidades B e C como início e fim do trajeto. Vamos inverter o trajeto B – C, que se tornará C – B, e teremos uma nova solução A – C – B – D. Se a distância do trajeto desta solução for menor do que a anterior, aceitamos, caso contrário, descartamos.

Na prática este procedimento desfaz trajetos que contenham rotas que cruzem sobre as outras, e por isso, são ineficientes. Em nosso exemplo vamos utilizar um 2-OPT aleatório devido ao tamanho do problema, mas uma outra maneira de decidir os elos que serão trocados é selecionar o primeiro elo que diminua a distância, verificando cada um deles, ou verificar exaustivamente e encontrar o que apresente a maior redução entre todos eles.

Simulated Annealing

Um dos problemas do método acima é que ele pode ficar parado em pontos mínimos locais. Existem vários métodos para “escapar” desses pontos, e buscar uma mínima melhor, mesmo sem garantias que ela seja global.

Uma destas maneiras é o Simulated Annealing. Baseado no processo de aquecimento e resfriamento da indústria de metalurgia, este procedimento possui dois parâmetros: temperatura e coeficiente de resfriamento. Basicamente, em vez de aceitar apenas trajetos que possuam uma distância menor, nós temos a probabilidade de aceitar um trajeto que piore, temporariamente, a distância, na busca por uma solução melhor no futuro.

    for i in xrange(1000000):

        new_solution = selectNextTour2Opt(solution,townCount)
        
        newCost = costSolution(new_solution,town_list,townCount)
        if newCost < obj:
            obj = newCost
            solution = new_solution[:]
        elif expCoinFlip(newCost,obj,T):
            solution = new_solution[:]
            obj = newCost


        if obj < min_obj:
            
            min_obj = obj
            print 'Epoch',i,'Objective', min_obj

        T = 0.999 * T



        if T <= final_T:
            T = init_T

Esta probabilidade é proporcional à temperatura. Pense da seguinte maneira: no início, com a temperatura alta, nós estamos abertos a explorar várias soluções, mesmo que elas não melhorem nosso objetivo, conforme o tempo vai passando, decidimos caminhar em direção à mínima local de uma destas soluções.

No nosso caso, ainda teremos o reaquecimento, que é simplesmente mandar o algoritmo voltar à temperatura inicial após atingirmos uma temperatura muito baixa.

Resultados

Este algoritmo requer o ajuste de alguns parâmetros para que possamos atingir o melhor resultado possível. O valor ideal para estes parâmetros depende da estrutura do problema, e uma das maneiras mais simples de se encontrá-los é testando algumas combinações.

Não quero complicar muito, então os parâmetros ajustados, e os respectivos valores testados serão:

Temperatura inicial: 1, 10
Coeficiente de resfriamento: 0,99; 0,999; 0,9999

Vamos ver quais são os valores que nos ajudam a alcançar a distância mínima em 1 milhão de iterações. Assim como este problema, encontrar os parâmetros é, normalmente, um problema com pontos ótimos locais. Se tivéssemos muitos valores para visitar, seria recomendado fazer uma busca aleatória, mas como estamos restritos a 6 combinações, e a avaliação leva pouco tempo, vou testar todas elas.

No fim, a combinação com a menor distância em 1 milhão de iterações foi a temperatura inicial 1, e o coeficiente de resfriamento 0,999. Mas como a diferença entre as opções foi pequena, acredito que tenha sido uma questão de sorte. A distância final foi 151.125 km, uma redução aproximada de 5.635 km, ou 3,59%.

Decidi rodar o mesmo script, com os parâmetros anteriores, por 10 milhões de iterações. Após cerca de 40 minutos a distância atingida foi 142.794 km, uma redução aproximada de 13.966 km, ou 8,91%. Por se tratar de um problema com muitas cidades, e cada uma delas bem próximas de outras, acredito que rodar por mais tempo ainda possa diminuir bastante esse número.

Seja o primeiro a saber das novidades em Machine Learning. Me siga no LinkedIn.