TSP_Matheuristics_Entrega.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "1fd06caf",
      "metadata": {
        "id": "1fd06caf"
      },
      "source": [
        "## Combinatorial Optimization\n",
        "\n",
        "### Solving The Symmetric Traveling Salesman Problem - TSP\n",
        "\n",
        "**Authors:** Guilherme C. Cadori; Rafael A. de Oliveira\n",
        "\n",
        "**Date:** 27/09/2023"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4348a811",
      "metadata": {
        "id": "4348a811"
      },
      "source": [
        "#### 1) Loading Instances"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "bd42090c",
      "metadata": {
        "id": "bd42090c"
      },
      "outputs": [],
      "source": [
        "# Importing libraries\n",
        "# !pip install gurobipy\n",
        "import numpy as np\n",
        "import networkx as nx\n",
        "import matplotlib.pyplot as plt\n",
        "import random as rdm\n",
        "import gurobipy as gp\n",
        "import time\n",
        "import requests\n",
        "import pandas as pd\n",
        "import itertools\n",
        "import numba as nb\n",
        "import scipy\n",
        "import operator\n",
        "import collections\n",
        "from io import StringIO\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "381efaf2",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "381efaf2",
        "outputId": "65c103ee-44f1-4e03-f2e6-c3cb851b01f1"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Loading instances:\n",
            "Loaded instance 1\n",
            "Loaded instance 2\n",
            "Loaded instance 3\n",
            "Loaded instance 4\n",
            "Loaded instance 5\n",
            "\n",
            "Successfully loaded 5 instances.\n",
            "\n",
            "The following instances were created:\n",
            "Instance29\n",
            "Instance48\n",
            "Instance52\n",
            "Instance100\n",
            "Instance130\n"
          ]
        }
      ],
      "source": [
        "# Define a list of URLs to fetch data from\n",
        "urls = [\n",
        "    'https://raw.githubusercontent.com/guilhermecadori/TSPData/main/bays29.tsp',\n",
        "    'https://raw.githubusercontent.com/guilhermecadori/TSPData/main/att48.tsp',\n",
        "    'https://raw.githubusercontent.com/guilhermecadori/TSPData/main/berlin52.tsp',\n",
        "    'https://raw.githubusercontent.com/guilhermecadori/TSPData/main/kroA100.tsp',\n",
        "    'https://raw.githubusercontent.com/guilhermecadori/TSPData/main/ch130.tsp'\n",
        "]\n",
        "\n",
        "# Create an empty dictionary to store datasets\n",
        "datasets = {}\n",
        "\n",
        "print('Loading instances:')\n",
        "\n",
        "# Loop through the URLs and fetch data for each one\n",
        "for url in urls:\n",
        "    # Send a GET request to the URL\n",
        "    response = requests.get(url)\n",
        "\n",
        "    # Check if the request was successful (status code 200)\n",
        "    if response.status_code == 200:\n",
        "        data = response.text\n",
        "\n",
        "        # Split the data into lines\n",
        "        lines = data.split('\\n')\n",
        "\n",
        "        # Get the dimension from the 4th row (assuming it follows the format \"DIMENSION: <value>\")\n",
        "        dimension_line = lines[3].strip()\n",
        "        dimension = int(dimension_line.split(\":\")[1].strip())\n",
        "\n",
        "        # Initialize an empty list to store the data\n",
        "        data_list = []\n",
        "\n",
        "        # Start reading data from the 7th line until \"EOF\" is encountered\n",
        "        for line in lines[6:]:\n",
        "            if line.strip() == 'EOF':\n",
        "                break\n",
        "            else:\n",
        "                columns = line.split()\n",
        "                data_list.append(columns)\n",
        "\n",
        "        # Create a Pandas DataFrame from the data list\n",
        "        df = pd.DataFrame(data_list, columns=['nodeID', 'coord1', 'coord3'])\n",
        "        df = df.astype({'nodeID': 'float64', 'coord1': 'float64', 'coord3': 'float64'})\n",
        "\n",
        "        # Create a dataset name that includes the dimension\n",
        "        dataset_name = f'Instance{dimension}'\n",
        "\n",
        "        # Store the dataset in the dictionary with the name as the key\n",
        "        datasets[dataset_name] = df\n",
        "\n",
        "        loadedInstanceCount = len(datasets.keys())\n",
        "\n",
        "        print(f'Loaded instance {loadedInstanceCount}')\n",
        "\n",
        "    else:\n",
        "        print(f'Failed to retrieve data from {url}. Status code: {response.status_code}')\n",
        "\n",
        "# Counting the all the loaded instances\n",
        "totalInstanceCount = len(datasets.keys())\n",
        "\n",
        "print(f'\\nSuccessfully loaded {totalInstanceCount} instances.')\n",
        "\n",
        "# Checking created datasets\n",
        "print('\\nThe following instances were created:')\n",
        "\n",
        "for dataset in datasets.keys():\n",
        "    print(dataset)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5f028bb8",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5f028bb8",
        "outputId": "069f4dc3-6723-478c-a292-faeec6571931"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "Calculations completed. Check if all expected matrices were produced\n",
            "\n",
            "\n",
            "The following distance matrices were created:\n",
            "Instance29\n",
            "Instance48\n",
            "Instance52\n",
            "Instance100\n",
            "Instance130\n"
          ]
        }
      ],
      "source": [
        "# Creating distance matrices\n",
        "distance_matrices = {}\n",
        "\n",
        "for dataset_name, dataset_df in datasets.items():\n",
        "    distance_matrices[dataset_name] = scipy.spatial.distance.cdist(dataset_df, dataset_df)\n",
        "\n",
        "print('\\nCalculations completed. Check if all expected matrices were produced\\n')\n",
        "\n",
        "# Now, 'distance_matrices' contains the distance matrices for each dataset\n",
        "# Checking matrices' names\n",
        "print('\\nThe following distance matrices were created:')\n",
        "\n",
        "for dataset in distance_matrices.keys():\n",
        "    print(dataset)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "eb407b38",
      "metadata": {
        "id": "eb407b38"
      },
      "outputs": [],
      "source": [
        "#### Instância 1\n",
        "# Creating pair-wise distance dictionaries\n",
        "# Number of nodes (assuming it's a square matrix)\n",
        "num_nodes1 = len(distance_matrices['Instance29'])\n",
        "\n",
        "# Initialize an empty dictionary to store the pair-wise distances\n",
        "pairwise_distances1 = {}\n",
        "\n",
        "# Iterate through each node\n",
        "for i in range(num_nodes1):\n",
        "    for j in range(num_nodes1):\n",
        "        # Skip calculating distance between the same node\n",
        "        if i == j:\n",
        "            continue\n",
        "        # Create a key for the pair (i, j) and store the distance\n",
        "        pairwise_distances1[(i, j)] = distance_matrices['Instance29'][i, j]\n",
        "\n",
        "# Now, pairwise_distances will contain the distances between each pair of nodes\n",
        "# pairwise_distances1\n",
        "\n",
        "\n",
        "#### Instância 2\n",
        "num_nodes2 = len(distance_matrices['Instance48'])\n",
        "\n",
        "# Initialize an empty dictionary to store the pair-wise distances\n",
        "pairwise_distances2 = {}\n",
        "\n",
        "# Iterate through each node\n",
        "for i in range(num_nodes2):\n",
        "    for j in range(num_nodes2):\n",
        "        # Skip calculating distance between the same node\n",
        "        if i == j:\n",
        "            continue\n",
        "        # Create a key for the pair (i, j) and store the distance\n",
        "        pairwise_distances2[(i, j)] = distance_matrices['Instance48'][i, j]\n",
        "\n",
        "\n",
        "#### Instância 3\n",
        "num_nodes3 = len(distance_matrices['Instance52'])\n",
        "\n",
        "# Initialize an empty dictionary to store the pair-wise distances\n",
        "pairwise_distances3 = {}\n",
        "\n",
        "# Iterate through each node\n",
        "for i in range(num_nodes3):\n",
        "    for j in range(num_nodes3):\n",
        "        # Skip calculating distance between the same node\n",
        "        if i == j:\n",
        "            continue\n",
        "        # Create a key for the pair (i, j) and store the distance\n",
        "        pairwise_distances3[(i, j)] = distance_matrices['Instance52'][i, j]\n",
        "\n",
        "\n",
        "#### Instância 4\n",
        "num_nodes4 = len(distance_matrices['Instance100'])\n",
        "\n",
        "# Initialize an empty dictionary to store the pair-wise distances\n",
        "pairwise_distances4 = {}\n",
        "\n",
        "# Iterate through each node\n",
        "for i in range(num_nodes4):\n",
        "    for j in range(num_nodes4):\n",
        "        # Skip calculating distance between the same node\n",
        "        if i == j:\n",
        "            continue\n",
        "        # Create a key for the pair (i, j) and store the distance\n",
        "        pairwise_distances4[(i, j)] = distance_matrices['Instance100'][i, j]\n",
        "\n",
        "\n",
        "#### Instância 5\n",
        "num_nodes5 = len(distance_matrices['Instance130'])\n",
        "\n",
        "# Initialize an empty dictionary to store the pair-wise distances\n",
        "pairwise_distances5 = {}\n",
        "\n",
        "# Iterate through each node\n",
        "for i in range(num_nodes5):\n",
        "    for j in range(num_nodes5):\n",
        "        # Skip calculating distance between the same node\n",
        "        if i == j:\n",
        "            continue\n",
        "        # Create a key for the pair (i, j) and store the distance\n",
        "        pairwise_distances5[(i, j)] = distance_matrices['Instance130'][i, j]\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "994a079e",
      "metadata": {
        "heading_collapsed": true,
        "id": "994a079e"
      },
      "source": [
        "#### 2) Solving Approach: Pure Heuristic"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8qPjid34ZD6M",
      "metadata": {
        "hidden": true,
        "id": "8qPjid34ZD6M"
      },
      "outputs": [],
      "source": [
        "def tour_length(instance, distance_matrices, tour):\n",
        "    total_distance = 0\n",
        "    for i in range(len(tour)):\n",
        "        total_distance += distance_matrices[instance][tour[i]][tour[(i+1) % len(tour)]]\n",
        "    return total_distance\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "-2UbTcSy64of",
      "metadata": {
        "hidden": true,
        "id": "-2UbTcSy64of"
      },
      "outputs": [],
      "source": [
        "def calculate_edge_changes(instance, distance_matrices, tour, i, j):\n",
        "    n = len(tour)\n",
        "    prev_i = (i - 1 + n) % n\n",
        "    next_i = (i + 1) % n\n",
        "    prev_j = (j - 1 + n) % n\n",
        "    next_j = (j + 1) % n\n",
        "\n",
        "    if abs(i-j) == 1 or abs(i-j) == n - 1:\n",
        "        old_length = (\n",
        "            distance_matrices[instance][tour[prev_i]][tour[i]] +\n",
        "            distance_matrices[instance][tour[i]][tour[next_i]] +\n",
        "            distance_matrices[instance][tour[prev_j]][tour[j]] +\n",
        "            distance_matrices[instance][tour[j]][tour[next_j]] -\n",
        "            distance_matrices[instance][tour[i]][tour[j]]\n",
        "        )\n",
        "        new_length = (\n",
        "            distance_matrices[instance][tour[prev_i]][tour[next_i]] +\n",
        "            distance_matrices[instance][tour[prev_j]][tour[next_j]] +\n",
        "            distance_matrices[instance][tour[i]][tour[j]]\n",
        "        )\n",
        "    else:\n",
        "        old_length = (\n",
        "            distance_matrices[instance][tour[prev_i]][tour[i]] +\n",
        "            distance_matrices[instance][tour[i]][tour[next_i]] +\n",
        "            distance_matrices[instance][tour[prev_j]][tour[j]] +\n",
        "            distance_matrices[instance][tour[j]][tour[next_j]]\n",
        "        )\n",
        "        new_length = (\n",
        "            distance_matrices[instance][tour[prev_j]][tour[i]] +\n",
        "            distance_matrices[instance][tour[i]][tour[next_j]] +\n",
        "            distance_matrices[instance][tour[prev_i]][tour[j]] +\n",
        "            distance_matrices[instance][tour[j]][tour[next_i]]\n",
        "        )\n",
        "\n",
        "    return new_length - old_length\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "UI2758NSZnsB",
      "metadata": {
        "hidden": true,
        "id": "UI2758NSZnsB"
      },
      "outputs": [],
      "source": [
        "def lin_kernighan(instance, distance_matrices, tour):\n",
        "    initial_time = time.time()\n",
        "    n = len(tour)\n",
        "    improvement = True\n",
        "\n",
        "    while improvement:\n",
        "        improvement = False\n",
        "\n",
        "        for i, j in itertools.combinations(range(n), 2):\n",
        "            delta_length = calculate_edge_changes(instance, distance_matrices, tour, i, j)\n",
        "\n",
        "            if delta_length < 0:\n",
        "                tour[i], tour[j] = tour[j], tour[i]\n",
        "                improvement = True\n",
        "                break\n",
        "\n",
        "    final_time = time.time()\n",
        "    elapsed_time = final_time - initial_time\n",
        "\n",
        "    return tour, elapsed_time\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "1knbmbpBElfQ",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "hidden": true,
        "id": "1knbmbpBElfQ",
        "outputId": "795c5aaf-4690-4b8b-ff8b-a8f493098a16"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Instance29 - Starting Solution: 22741.77 | Final Solution: 10877.69 | Elapsed Time: 0.03\n",
            "Instance48 - Starting Solution: 162121.67 | Final Solution: 42608.25 | Elapsed Time: 0.09\n",
            "Instance52 - Starting Solution: 29138.6 | Final Solution: 9945.52 | Elapsed Time: 0.14\n",
            "Instance100 - Starting Solution: 179500.45 | Final Solution: 42786.04 | Elapsed Time: 0.52\n",
            "Instance130 - Starting Solution: 49541.99 | Final Solution: 12686.44 | Elapsed Time: 2.82\n"
          ]
        }
      ],
      "source": [
        "instance_tours = []\n",
        "instance_tours_length = []\n",
        "for instance in datasets.keys():\n",
        "    tours = []\n",
        "    tours_length = []\n",
        "    elapsed_time = 0\n",
        "    best_tour_length = float('inf')\n",
        "    for i in range(20):\n",
        "        initial_tour = rdm.sample(list(range(len(datasets[instance]))), len(datasets[instance]))\n",
        "        initial_length = tour_length(instance, distance_matrices, initial_tour)\n",
        "        improved_tour, elapsed_time_one = lin_kernighan(instance, distance_matrices, initial_tour)\n",
        "        improved_length = tour_length(instance, distance_matrices, improved_tour)\n",
        "        tours.append(improved_tour)\n",
        "        tours_length.append(improved_length)\n",
        "        elapsed_time += elapsed_time_one\n",
        "        if improved_length < best_tour_length:\n",
        "            best_tour_length = improved_length\n",
        "    instance_tours.append(tours)\n",
        "    instance_tours_length.append(tours_length)\n",
        "    print(instance,' - ', 'Starting Solution: ', round(initial_length, 2),' | ',\n",
        "          'Final Solution: ', round(best_tour_length, 2), ' | ',\n",
        "          'Elapsed Time: ', round(elapsed_time, 2), sep='')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "t-tR5i-Gnu-4",
      "metadata": {
        "id": "t-tR5i-Gnu-4"
      },
      "outputs": [],
      "source": [
        "count_arcs_instances = {}\n",
        "# Itera as Instâncias\n",
        "for idx, instance in enumerate(datasets.keys()):\n",
        "    max_reps = {}\n",
        "    solutions = np.array(instance_tours[idx])\n",
        "    # Itera as Soluções das Instâncias\n",
        "    arcs = []\n",
        "    count_arcs = {}\n",
        "    for solution in solutions:\n",
        "        # Itera os Nós da Solução\n",
        "        for idx, node in enumerate(solution):\n",
        "            arcs.append((solution[idx], solution[(idx + 1) % len(solution)]))\n",
        "    for arc in arcs:\n",
        "        if arc in count_arcs.keys():\n",
        "            continue\n",
        "        else:\n",
        "            count_arcs[arc] = arcs.count(arc)\n",
        "    count_arcs_instances[instance] = dict(sorted(count_arcs.items(), key=lambda item: item[1], reverse=True))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c07496ae",
      "metadata": {
        "id": "c07496ae"
      },
      "source": [
        "#### 3) Solving Approach: Pure Mathematical Model"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "53a81619",
      "metadata": {
        "id": "53a81619"
      },
      "source": [
        "##### 3.1) Instância 1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e7d04a71",
      "metadata": {
        "id": "e7d04a71",
        "outputId": "d88efa23-61df-40a5-e760-01c770f6118d",
        "scrolled": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 29 rows, 812 columns and 812 nonzeros\n",
            "Model fingerprint: 0xc0a0dc88\n",
            "Variable types: 0 continuous, 812 integer (812 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [1e+02, 2e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 2e+00]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 29 rows, 812 columns, 812 nonzeros\n",
            "Variable types: 0 continuous, 812 integer (812 binary)\n",
            "\n",
            "Root relaxation: objective 8.655120e+03, 44 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 8655.12021    0    6          - 8655.12021      -     -    0s\n",
            "     0     0 8713.77221    0    6          - 8713.77221      -     -    0s\n",
            "     0     0 8721.60291    0   14          - 8721.60291      -     -    0s\n",
            "     0     0 8727.68287    0    -          - 8727.68287      -     -    0s\n",
            "     0     0 8967.46702    0   16          - 8967.46702      -     -    0s\n",
            "     0     0 8967.46702    0   16          - 8967.46702      -     -    0s\n",
            "     0     2 8967.46702    0   16          - 8967.46702      -     -    0s\n",
            "*   18    12               2    9131.0426923 9084.17020  0.51%   5.8    0s\n",
            "*   23     6               4    9119.6993441 9084.17020  0.39%   5.2    0s\n",
            "*   25     6               6    9105.5402591 9084.17020  0.23%   5.0    0s\n",
            "*   33     0               4    9085.8315732 9085.83157  0.00%   3.8    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Zero half: 5\n",
            "  Lazy constraints: 7\n",
            "\n",
            "Explored 36 nodes (193 simplex iterations) in 0.06 seconds (0.01 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 4: 9085.83 9105.54 9119.7 9131.04 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 9.085831573160e+03, best bound 9.085831573160e+03, gap 0.0000%\n",
            "\n",
            "User-callback calls 214, time in user-callback 0.03 sec\n",
            "____________________________________________________________________________________________________\n",
            "\n",
            "Instância 1 - Função Objetivo: R$ 9085.83\n",
            "\n",
            "Elapsed time: 0.077 seconds\n",
            "\n",
            "Rota: [0, 23, 12, 15, 26, 7, 22, 6, 24, 18, 10, 21, 16, 13, 17, 14, 3, 9, 19, 1, 20, 4, 28, 2, 25, 8, 11, 5, 27]\n",
            "____________________________________________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance29']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance29'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances1.keys(),\n",
        "                     obj=pairwise_distances1,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances1.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('_' * 100)\n",
        "#### Retornando solução\n",
        "print('\\nInstância 1 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('_' * 100)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f5a6b4f",
      "metadata": {
        "id": "4f5a6b4f"
      },
      "source": [
        "##### 3.2) Instância 2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "40a45133",
      "metadata": {
        "id": "40a45133",
        "outputId": "553cce27-2e07-4d18-86ad-93ec949baa81"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 48 rows, 2256 columns and 2256 nonzeros\n",
            "Model fingerprint: 0x7aeef8e3\n",
            "Variable types: 0 continuous, 2256 integer (2256 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [1e+02, 8e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 2e+00]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 48 rows, 2256 columns, 2256 nonzeros\n",
            "Variable types: 0 continuous, 2256 integer (2256 binary)\n",
            "\n",
            "Root relaxation: objective 3.169020e+04, 77 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 31690.2032    0   14          - 31690.2032      -     -    0s\n",
            "     0     0 31742.6322    0    6          - 31742.6322      -     -    0s\n",
            "     0     0 31776.6984    0   14          - 31776.6984      -     -    0s\n",
            "     0     0 31802.2087    0    8          - 31802.2087      -     -    0s\n",
            "     0     0 31802.2087    0    8          - 31802.2087      -     -    0s\n",
            "     0     2 31802.2087    0    8          - 31802.2087      -     -    0s\n",
            "*  164   100              11    33801.942660 32933.7176  2.57%   6.1    0s\n",
            "*  234    70               7    33567.062881 33156.1869  1.22%   6.5    0s\n",
            "*  305    14               8    33541.935605 33531.5069  0.03%   6.0    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Zero half: 7\n",
            "  Lazy constraints: 22\n",
            "\n",
            "Explored 333 nodes (1959 simplex iterations) in 0.15 seconds (0.05 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 3: 33541.9 33567.1 33801.9 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 3.354193560452e+04, best bound 3.354193560452e+04, gap 0.0000%\n",
            "\n",
            "User-callback calls 829, time in user-callback 0.04 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 2 - Função Objetivo: R$ 33541.94\n",
            "\n",
            "Elapsed time: 0.1514 seconds\n",
            "\n",
            "Rota: [0, 7, 37, 30, 43, 17, 6, 27, 5, 36, 18, 26, 16, 42, 29, 35, 45, 32, 19, 46, 20, 31, 38, 47, 4, 41, 23, 9, 44, 34, 3, 25, 1, 28, 33, 40, 15, 21, 2, 22, 13, 24, 12, 10, 11, 14, 39, 8]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance48']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance48'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances2.keys(),\n",
        "                     obj=pairwise_distances2,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances2.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 2 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1b0fe8c9",
      "metadata": {
        "id": "1b0fe8c9"
      },
      "source": [
        "##### 3.3) Instância 3"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "edddf512",
      "metadata": {
        "id": "edddf512",
        "outputId": "1eaebe56-ff93-49f9-e828-8053e92ea6d6"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 52 rows, 2652 columns and 2652 nonzeros\n",
            "Model fingerprint: 0xf218a4b2\n",
            "Variable types: 0 continuous, 2652 integer (2652 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [2e+01, 2e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 2e+00]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 52 rows, 2652 columns, 2652 nonzeros\n",
            "Variable types: 0 continuous, 2652 integer (2652 binary)\n",
            "\n",
            "Root relaxation: objective 7.281303e+03, 88 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 7281.30305    0    -          - 7281.30305      -     -    0s\n",
            "     0     0 7282.63738    0    6          - 7282.63738      -     -    0s\n",
            "     0     0 7283.97740    0    6          - 7283.97740      -     -    0s\n",
            "     0     0 7285.87733    0    -          - 7285.87733      -     -    0s\n",
            "     0     0 7287.47024    0    -          - 7287.47024      -     -    0s\n",
            "     0     0 7617.82768    0    -          - 7617.82768      -     -    0s\n",
            "*    0     0               0    7649.8377518 7649.83775  0.00%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Zero half: 4\n",
            "  Lazy constraints: 4\n",
            "\n",
            "Explored 1 nodes (103 simplex iterations) in 0.03 seconds (0.01 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 1: 7649.84 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 7.649837751833e+03, best bound 7.649837751833e+03, gap 0.0000%\n",
            "\n",
            "User-callback calls 120, time in user-callback 0.03 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 3 - Função Objetivo: R$ 7649.84\n",
            "\n",
            "Elapsed time: 0.04 seconds\n",
            "\n",
            "Rota: [0, 21, 30, 17, 2, 16, 20, 41, 6, 1, 29, 22, 19, 49, 28, 15, 45, 43, 33, 34, 35, 38, 39, 36, 37, 47, 23, 4, 14, 5, 3, 24, 11, 27, 26, 25, 46, 12, 13, 51, 10, 50, 32, 42, 9, 8, 7, 40, 18, 44, 31, 48]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance52']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance52'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances3.keys(),\n",
        "                     obj=pairwise_distances3,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances3.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 3 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "01a81c35",
      "metadata": {
        "id": "01a81c35"
      },
      "source": [
        "##### 3.4) Instância 4"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "37305b3b",
      "metadata": {
        "id": "37305b3b",
        "outputId": "903fb9bb-21a5-44fd-a58d-62d2c0990349"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 100 rows, 9900 columns and 9900 nonzeros\n",
            "Model fingerprint: 0xc1ce0230\n",
            "Variable types: 0 continuous, 9900 integer (9900 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [1e+01, 4e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 2e+00]\n",
            "Presolve time: 0.01s\n",
            "Presolved: 100 rows, 9900 columns, 9900 nonzeros\n",
            "Variable types: 0 continuous, 9900 integer (9900 binary)\n",
            "\n",
            "Root relaxation: objective 1.994385e+04, 155 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 19943.8534    0   18          - 19943.8534      -     -    0s\n",
            "     0     0 20015.7990    0    8          - 20015.7990      -     -    0s\n",
            "     0     0 20054.5297    0   17          - 20054.5297      -     -    0s\n",
            "     0     0 20076.4222    0   19          - 20076.4222      -     -    0s\n",
            "     0     0 20076.4222    0   19          - 20076.4222      -     -    0s\n",
            "     0     2 20076.4222    0   19          - 20076.4222      -     -    0s\n",
            "* 1780  1283              47    22252.926949 21405.9602  3.81%   5.0    1s\n",
            "* 1786  1202              49    22170.461164 21405.9602  3.45%   5.0    1s\n",
            "* 2148  1001              42    22019.555961 21481.3649  2.44%   4.8    1s\n",
            "* 2298   866              35    21949.970526 21482.9647  2.13%   4.7    1s\n",
            "* 3037   550              30    21883.268378 21705.8681  0.81%   4.8    2s\n",
            "* 3225   392              34    21864.367034 21751.3313  0.52%   4.8    2s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 3\n",
            "  Lift-and-project: 2\n",
            "  Inf proof: 2\n",
            "  Zero half: 8\n",
            "  Lazy constraints: 21\n",
            "\n",
            "Explored 3850 nodes (18080 simplex iterations) in 2.23 seconds (1.18 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 6: 21864.4 21883.3 21950 ... 22252.9\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 2.186436703366e+04, best bound 2.186436703366e+04, gap 0.0000%\n",
            "\n",
            "User-callback calls 8112, time in user-callback 0.44 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 4 - Função Objetivo: R$ 21864.37\n",
            "\n",
            "Elapsed time: 2.2474 seconds\n",
            "\n",
            "Rota: [0, 46, 92, 27, 66, 57, 60, 50, 86, 24, 80, 68, 63, 39, 53, 1, 43, 49, 72, 67, 84, 81, 94, 12, 75, 32, 36, 4, 51, 77, 95, 38, 29, 47, 99, 40, 70, 13, 2, 42, 45, 28, 33, 82, 54, 6, 8, 56, 19, 11, 26, 85, 34, 61, 59, 76, 22, 97, 90, 44, 31, 10, 14, 16, 58, 73, 20, 71, 9, 83, 35, 98, 37, 23, 17, 78, 52, 87, 93, 15, 21, 69, 65, 25, 64, 3, 96, 55, 79, 30, 88, 41, 7, 91, 74, 18, 89, 48, 5, 62]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance100']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance100'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances4.keys(),\n",
        "                     obj=pairwise_distances4,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances4.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 4 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "518dc49e",
      "metadata": {
        "id": "518dc49e"
      },
      "source": [
        "##### 3.5) Instância 5"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "98e740ca",
      "metadata": {
        "id": "98e740ca",
        "outputId": "8f861841-8e16-4ae1-8571-97130dfd26d8"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 130 rows, 16770 columns and 16770 nonzeros\n",
            "Model fingerprint: 0xb8fad695\n",
            "Variable types: 0 continuous, 16770 integer (16770 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 1e+00]\n",
            "  Objective range  [7e+00, 9e+02]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 2e+00]\n",
            "Presolve time: 0.01s\n",
            "Presolved: 130 rows, 16770 columns, 16770 nonzeros\n",
            "Variable types: 0 continuous, 16770 integer (16770 binary)\n",
            "\n",
            "Root relaxation: objective 7.429119e+03, 203 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 7429.11915    0   26          - 7429.11915      -     -    0s\n",
            "     0     0 7450.84097    0   24          - 7450.84097      -     -    0s\n",
            "     0     0 7459.32375    0   27          - 7459.32375      -     -    0s\n",
            "     0     0 7466.14773    0   19          - 7466.14773      -     -    0s\n",
            "     0     0 7466.60802    0   26          - 7466.60802      -     -    0s\n",
            "     0     0 7466.60802    0   26          - 7466.60802      -     -    0s\n",
            "     0     2 7466.60802    0   26          - 7466.60802      -     -    0s\n",
            "* 1034   906              10    7813.0549906 7813.05499  0.00%   6.1    1s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 5\n",
            "  Zero half: 7\n",
            "  Lazy constraints: 1\n",
            "\n",
            "Explored 1034 nodes (6806 simplex iterations) in 1.44 seconds (0.96 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 1: 7813.05 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 7.813054990600e+03, best bound 7.813054990600e+03, gap 0.0000%\n",
            "\n",
            "User-callback calls 2346, time in user-callback 0.45 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 5 - Função Objetivo: R$ 7813.05\n",
            "\n",
            "Elapsed time: 1.4503 seconds\n",
            "\n",
            "Rota: [0, 15, 44, 4, 10, 75, 108, 60, 128, 123, 63, 68, 85, 87, 25, 6, 96, 69, 106, 126, 103, 42, 33, 16, 30, 26, 18, 14, 28, 23, 99, 115, 94, 78, 37, 11, 86, 80, 102, 76, 93, 88, 109, 97, 67, 62, 47, 24, 31, 35, 83, 112, 118, 110, 122, 100, 81, 56, 55, 64, 8, 51, 74, 73, 98, 72, 91, 105, 119, 52, 57, 48, 71, 90, 101, 5, 9, 13, 12, 66, 95, 121, 54, 59, 50, 41, 43, 92, 46, 39, 36, 21, 22, 32, 20, 17, 7, 2, 29, 58, 89, 124, 84, 65, 77, 120, 125, 107, 113, 82, 114, 27, 61, 104, 127, 111, 116, 129, 117, 79, 45, 19, 3, 34, 53, 1, 49, 70, 38, 40]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance130']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance130'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances5.keys(),\n",
        "                     obj=pairwise_distances5,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances5.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 5 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8306934d",
      "metadata": {
        "id": "8306934d"
      },
      "source": [
        "#### 4) Solving Approach: Matheuristics"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e6fdaa2e",
      "metadata": {
        "id": "e6fdaa2e"
      },
      "source": [
        "##### Instance I"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "264ac9c3",
      "metadata": {
        "id": "264ac9c3",
        "outputId": "99e57026-7f06-46f8-c690-0fa21c5274c3"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 30 rows, 812 columns and 904 nonzeros\n",
            "Model fingerprint: 0xf6cebd38\n",
            "Variable types: 0 continuous, 812 integer (812 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 2e+00]\n",
            "  Objective range  [1e+02, 2e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 1e+01]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 30 rows, 812 columns, 904 nonzeros\n",
            "Variable types: 0 continuous, 812 integer (812 binary)\n",
            "\n",
            "Root relaxation: objective 8.655120e+03, 46 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 8655.12021    0    6          - 8655.12021      -     -    0s\n",
            "     0     0 8713.77221    0    6          - 8713.77221      -     -    0s\n",
            "     0     0 8721.60291    0   14          - 8721.60291      -     -    0s\n",
            "     0     0 8727.68287    0    -          - 8727.68287      -     -    0s\n",
            "     0     0 8967.46702    0   16          - 8967.46702      -     -    0s\n",
            "     0     0 8967.46702    0   16          - 8967.46702      -     -    0s\n",
            "     0     2 8967.46702    0   16          - 8967.46702      -     -    0s\n",
            "*   18    12               2    9131.0426923 9084.17020  0.51%   5.9    0s\n",
            "*   23     6               4    9119.6993441 9084.17020  0.39%   5.3    0s\n",
            "*   25     6               6    9105.5402591 9084.17020  0.23%   5.1    0s\n",
            "*   33     0               4    9085.8315732 9085.83157  0.00%   4.0    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Zero half: 5\n",
            "  Lazy constraints: 7\n",
            "\n",
            "Explored 36 nodes (206 simplex iterations) in 0.06 seconds (0.01 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 4: 9085.83 9105.54 9119.7 9131.04 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 9.085831573160e+03, best bound 9.085831573160e+03, gap 0.0000%\n",
            "\n",
            "User-callback calls 214, time in user-callback 0.00 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 1 - Função Objetivo: R$ 9085.83\n",
            "\n",
            "Elapsed time: 0.0626 seconds\n",
            "\n",
            "Rota: [0, 23, 12, 15, 26, 7, 22, 6, 24, 18, 10, 21, 16, 13, 17, 14, 3, 9, 19, 1, 20, 4, 28, 2, 25, 8, 11, 5, 27]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance29']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance29'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances1.keys(),\n",
        "                     obj=pairwise_distances1,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances1.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "# Creating Restricted Space constraints\n",
        "# Define the maximum number of arcs allowed as a percentage of the total number of nodes\n",
        "max_arcs = 10 #/ 100 * (2 * num_nodes - 1)\n",
        "\n",
        "# Dictionary to store constraints\n",
        "constraints = {}\n",
        "\n",
        "# Access nodes specific to 'Instance29'\n",
        "nodes_instance29 = count_arcs_instances.get('Instance29', {})\n",
        "\n",
        "# Add constraints for 'Instance29'\n",
        "instance29_variables = [variables[i, j] for (i, j) in nodes_instance29.keys()]\n",
        "if len(instance29_variables) >= max_arcs:\n",
        "    constraint = m1.addConstr(gp.quicksum(instance29_variables) >= max_arcs, name=f'Constraint_Instance29')\n",
        "    constraints['Instance29'] = constraint\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Running optimization\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 1 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2ac8f6c2",
      "metadata": {
        "id": "2ac8f6c2"
      },
      "source": [
        "##### Instance II"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0c7eb3b6",
      "metadata": {
        "id": "0c7eb3b6",
        "outputId": "b071cad4-c5c6-4e3b-fe5f-3327baa170be"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 49 rows, 2256 columns and 2414 nonzeros\n",
            "Model fingerprint: 0x198f5bbd\n",
            "Variable types: 0 continuous, 2256 integer (2256 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 2e+00]\n",
            "  Objective range  [1e+02, 8e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 1e+01]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 49 rows, 2256 columns, 2414 nonzeros\n",
            "Variable types: 0 continuous, 2256 integer (2256 binary)\n",
            "\n",
            "Root relaxation: objective 3.169020e+04, 78 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 31690.2032    0   14          - 31690.2032      -     -    0s\n",
            "     0     0 31742.6322    0    6          - 31742.6322      -     -    0s\n",
            "     0     0 31776.6984    0   14          - 31776.6984      -     -    0s\n",
            "     0     0 31802.2087    0    8          - 31802.2087      -     -    0s\n",
            "     0     0 31802.2087    0    8          - 31802.2087      -     -    0s\n",
            "     0     2 31802.2087    0    8          - 31802.2087      -     -    0s\n",
            "*  109    78               9    33541.935605 32537.2081  3.00%   6.3    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Zero half: 7\n",
            "  Lazy constraints: 19\n",
            "\n",
            "Explored 229 nodes (1438 simplex iterations) in 0.12 seconds (0.05 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 1: 33541.9 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 3.354193560452e+04, best bound 3.354193560452e+04, gap 0.0000%\n",
            "\n",
            "User-callback calls 609, time in user-callback 0.06 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 2 - Função Objetivo: R$ 33541.94\n",
            "\n",
            "Elapsed time: 0.1316 seconds\n",
            "\n",
            "Rota: [0, 7, 37, 30, 43, 17, 6, 27, 5, 36, 18, 26, 16, 42, 29, 35, 45, 32, 19, 46, 20, 31, 38, 47, 4, 41, 23, 9, 44, 34, 3, 25, 1, 28, 33, 40, 15, 21, 2, 22, 13, 24, 12, 10, 11, 14, 39, 8]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance48']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance48'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances2.keys(),\n",
        "                     obj=pairwise_distances2,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances2.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "# Creating Restricted Space constraints\n",
        "# Define the maximum number of arcs allowed as a percentage of the total number of nodes\n",
        "max_arcs = 10 #/ 100 * (2 * num_nodes - 1)\n",
        "\n",
        "# Dictionary to store constraints\n",
        "constraints = {}\n",
        "\n",
        "# Access nodes specific to 'Instance29'\n",
        "nodes_instance29 = count_arcs_instances.get('Instance48', {})\n",
        "\n",
        "# Add constraints for 'Instance29'\n",
        "instance29_variables = [variables[i, j] for (i, j) in nodes_instance29.keys()]\n",
        "if len(instance29_variables) >= max_arcs:\n",
        "    constraint = m1.addConstr(gp.quicksum(instance29_variables) >= max_arcs, name=f'Constraint_Instance48')\n",
        "    constraints['Instance48'] = constraint\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Running optimization\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 2 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "40222c82",
      "metadata": {
        "id": "40222c82"
      },
      "source": [
        "##### Instance III"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f922468b",
      "metadata": {
        "id": "f922468b",
        "outputId": "7427be56-e53f-461f-8d6d-b94d7945c9b9"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 53 rows, 2652 columns and 2831 nonzeros\n",
            "Model fingerprint: 0xf419209e\n",
            "Variable types: 0 continuous, 2652 integer (2652 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 2e+00]\n",
            "  Objective range  [2e+01, 2e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 1e+01]\n",
            "Presolve time: 0.00s\n",
            "Presolved: 53 rows, 2652 columns, 2831 nonzeros\n",
            "Variable types: 0 continuous, 2652 integer (2652 binary)\n",
            "\n",
            "Root relaxation: objective 7.281303e+03, 88 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 7281.30305    0    -          - 7281.30305      -     -    0s\n",
            "     0     0 7282.63738    0    6          - 7282.63738      -     -    0s\n",
            "     0     0 7283.97740    0    6          - 7283.97740      -     -    0s\n",
            "     0     0 7285.87733    0    -          - 7285.87733      -     -    0s\n",
            "     0     0 7287.47024    0    -          - 7287.47024      -     -    0s\n",
            "     0     0 7617.82768    0    -          - 7617.82768      -     -    0s\n",
            "*    0     0               0    7649.8377518 7649.83775  0.00%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Zero half: 4\n",
            "  Lazy constraints: 4\n",
            "\n",
            "Explored 1 nodes (104 simplex iterations) in 0.04 seconds (0.01 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 1: 7649.84 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 7.649837751833e+03, best bound 7.649837751833e+03, gap 0.0000%\n",
            "\n",
            "User-callback calls 121, time in user-callback 0.02 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 3 - Função Objetivo: R$ 7649.84\n",
            "\n",
            "Elapsed time: 0.048 seconds\n",
            "\n",
            "Rota: [0, 21, 30, 17, 2, 16, 20, 41, 6, 1, 29, 22, 19, 49, 28, 15, 45, 43, 33, 34, 35, 38, 39, 36, 37, 47, 23, 4, 14, 5, 3, 24, 11, 27, 26, 25, 46, 12, 13, 51, 10, 50, 32, 42, 9, 8, 7, 40, 18, 44, 31, 48]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance52']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance52'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances3.keys(),\n",
        "                     obj=pairwise_distances3,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances3.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "# Creating Restricted Space constraints\n",
        "# Define the maximum number of arcs allowed as a percentage of the total number of nodes\n",
        "max_arcs = 10 #/ 100 * (2 * num_nodes - 1)\n",
        "\n",
        "# Dictionary to store constraints\n",
        "constraints = {}\n",
        "\n",
        "# Access nodes specific to 'Instance29'\n",
        "nodes_instance29 = count_arcs_instances.get('Instance52', {})\n",
        "\n",
        "# Add constraints for 'Instance29'\n",
        "instance29_variables = [variables[i, j] for (i, j) in nodes_instance29.keys()]\n",
        "if len(instance29_variables) >= max_arcs:\n",
        "    constraint = m1.addConstr(gp.quicksum(instance29_variables) >= max_arcs, name=f'Constraint_Instance52')\n",
        "    constraints['Instance52'] = constraint\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Running optimization\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 3 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e7ac6861",
      "metadata": {
        "id": "e7ac6861"
      },
      "source": [
        "##### Instance IV"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e4df284e",
      "metadata": {
        "id": "e4df284e",
        "outputId": "64517f46-be76-4993-d98a-4dfd6674dbea"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 101 rows, 9900 columns and 10317 nonzeros\n",
            "Model fingerprint: 0x6814d07c\n",
            "Variable types: 0 continuous, 9900 integer (9900 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 2e+00]\n",
            "  Objective range  [1e+01, 4e+03]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 1e+01]\n",
            "Presolve time: 0.01s\n",
            "Presolved: 101 rows, 9900 columns, 10317 nonzeros\n",
            "Variable types: 0 continuous, 9900 integer (9900 binary)\n",
            "\n",
            "Root relaxation: objective 1.994385e+04, 147 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 19943.8534    0   18          - 19943.8534      -     -    0s\n",
            "     0     0 20015.7990    0    8          - 20015.7990      -     -    0s\n",
            "     0     0 20054.5297    0   17          - 20054.5297      -     -    0s\n",
            "     0     0 20076.4222    0   19          - 20076.4222      -     -    0s\n",
            "     0     0 20076.4222    0   19          - 20076.4222      -     -    0s\n",
            "     0     2 20076.4222    0   19          - 20076.4222      -     -    0s\n",
            "* 1975  1771              21    21892.251371 21842.4049  0.23%   5.0    1s\n",
            "* 2009  1644              20    21883.617856 21842.4049  0.19%   4.9    1s\n",
            "* 2011  1555              19    21864.367034 21861.2112  0.01%   4.9    1s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 14\n",
            "  MIR: 4\n",
            "  Flow cover: 4\n",
            "  Zero half: 8\n",
            "  Lazy constraints: 11\n",
            "\n",
            "Explored 2017 nodes (10063 simplex iterations) in 1.42 seconds (1.10 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 3: 21864.4 21883.6 21892.3 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 2.186436703366e+04, best bound 2.186436703366e+04, gap 0.0000%\n",
            "\n",
            "User-callback calls 4400, time in user-callback 0.36 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 4 - Função Objetivo: R$ 21864.37\n",
            "\n",
            "Elapsed time: 1.4236 seconds\n",
            "\n",
            "Rota: [0, 46, 92, 27, 66, 57, 60, 50, 86, 24, 80, 68, 63, 39, 53, 1, 43, 49, 72, 67, 84, 81, 94, 12, 75, 32, 36, 4, 51, 77, 95, 38, 29, 47, 99, 40, 70, 13, 2, 42, 45, 28, 33, 82, 54, 6, 8, 56, 19, 11, 26, 85, 34, 61, 59, 76, 22, 97, 90, 44, 31, 10, 14, 16, 58, 73, 20, 71, 9, 83, 35, 98, 37, 23, 17, 78, 52, 87, 93, 15, 21, 69, 65, 25, 64, 3, 96, 55, 79, 30, 88, 41, 7, 91, 74, 18, 89, 48, 5, 62]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance100']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance100'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances4.keys(),\n",
        "                     obj=pairwise_distances4,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances4.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "# Creating Restricted Space constraints\n",
        "# Define the maximum number of arcs allowed as a percentage of the total number of nodes\n",
        "max_arcs = 10 #/ 100 * (2 * num_nodes - 1)\n",
        "\n",
        "# Dictionary to store constraints\n",
        "constraints = {}\n",
        "\n",
        "# Access nodes specific to 'Instance29'\n",
        "nodes_instance29 = count_arcs_instances.get('Instance100', {})\n",
        "\n",
        "# Add constraints for 'Instance29'\n",
        "instance29_variables = [variables[i, j] for (i, j) in nodes_instance29.keys()]\n",
        "if len(instance29_variables) >= max_arcs:\n",
        "    constraint = m1.addConstr(gp.quicksum(instance29_variables) >= max_arcs, name=f'Constraint_Instance100')\n",
        "    constraints['Instance100'] = constraint\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Running optimization\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 4 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ebbd7e6f",
      "metadata": {
        "id": "ebbd7e6f"
      },
      "source": [
        "##### Instance V"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "69052123",
      "metadata": {
        "id": "69052123",
        "outputId": "51387653-f103-41ab-9062-e22dd8b20062"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter LazyConstraints to value 1\n",
            "Set parameter Heuristics to value 0\n",
            "Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (win64)\n",
            "\n",
            "CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]\n",
            "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
            "\n",
            "Optimize a model with 131 rows, 16770 columns and 17333 nonzeros\n",
            "Model fingerprint: 0x9de9984a\n",
            "Variable types: 0 continuous, 16770 integer (16770 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e+00, 2e+00]\n",
            "  Objective range  [7e+00, 9e+02]\n",
            "  Bounds range     [1e+00, 1e+00]\n",
            "  RHS range        [2e+00, 1e+01]\n",
            "Presolve time: 0.01s\n",
            "Presolved: 131 rows, 16770 columns, 17333 nonzeros\n",
            "Variable types: 0 continuous, 16770 integer (16770 binary)\n",
            "\n",
            "Root relaxation: objective 7.429119e+03, 225 iterations, 0.01 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 7429.11915    0   26          - 7429.11915      -     -    0s\n",
            "     0     0 7450.84097    0   24          - 7450.84097      -     -    0s\n",
            "     0     0 7459.32375    0   27          - 7459.32375      -     -    0s\n",
            "     0     0 7466.14773    0   19          - 7466.14773      -     -    0s\n",
            "     0     0 7466.14773    0   19          - 7466.14773      -     -    0s\n",
            "     0     2 7466.14773    0   19          - 7466.14773      -     -    0s\n",
            "* 1642  1478              12    7813.0549906 7813.05499  0.00%   6.1    1s\n",
            "\n",
            "Cutting planes:\n",
            "  Gomory: 5\n",
            "  Lift-and-project: 6\n",
            "  Flow cover: 1\n",
            "  Zero half: 9\n",
            "\n",
            "Explored 1642 nodes (10543 simplex iterations) in 1.86 seconds (1.48 work units)\n",
            "Thread count was 8 (of 8 available processors)\n",
            "\n",
            "Solution count 1: 7813.05 \n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 7.813054990600e+03, best bound 7.813054990600e+03, gap 0.0000%\n",
            "\n",
            "User-callback calls 3634, time in user-callback 0.39 sec\n",
            "________________________________________________________________\n",
            "\n",
            "Instância 5 - Função Objetivo: R$ 7813.05\n",
            "\n",
            "Elapsed time: 1.882 seconds\n",
            "\n",
            "Rota: [0, 15, 44, 4, 10, 75, 108, 60, 128, 123, 63, 68, 85, 87, 25, 6, 96, 69, 106, 126, 103, 42, 33, 16, 30, 26, 18, 14, 28, 23, 99, 115, 94, 78, 37, 11, 86, 80, 102, 76, 93, 88, 109, 97, 67, 62, 47, 24, 31, 35, 83, 112, 118, 110, 122, 100, 81, 56, 55, 64, 8, 51, 74, 73, 98, 72, 91, 105, 119, 52, 57, 48, 71, 90, 101, 5, 9, 13, 12, 66, 95, 121, 54, 59, 50, 41, 43, 92, 46, 39, 36, 21, 22, 32, 20, 17, 7, 2, 29, 58, 89, 124, 84, 65, 77, 120, 125, 107, 113, 82, 114, 27, 61, 104, 127, 111, 116, 129, 117, 79, 45, 19, 3, 34, 53, 1, 49, 70, 38, 40]\n",
            "________________________________________________________________\n"
          ]
        }
      ],
      "source": [
        "# Instance range\n",
        "inst1 = range(len(distance_matrices['Instance130']))\n",
        "# print(inst1)\n",
        "\n",
        "# Instance size\n",
        "num_nodes = distance_matrices['Instance130'].shape[0]\n",
        "node_ids1 = list(range(num_nodes))\n",
        "# print(nodesInst1)\n",
        "\n",
        "# Criando o modelo\n",
        "m1 = gp.Model()\n",
        "\n",
        "# Criando variávies: nesse caso todos o nós são adjacente entre si\n",
        "variables = m1.addVars(pairwise_distances5.keys(),\n",
        "                     obj=pairwise_distances5,\n",
        "                     vtype=gp.GRB.BINARY,\n",
        "                     name='x')\n",
        "\n",
        "# Direções simétricas: copiando objeto\n",
        "for u, v in pairwise_distances5.keys():\n",
        "    variables[v, u] = variables[u, v] # mesmo arco porém na direção oposta\n",
        "\n",
        "# Criando restrições\n",
        "# Em cada nó pode haver somente dois arcos selecionados - um nó de chegada e outro nó de saída\n",
        "cons = m1.addConstrs(variables.sum(node,'*') == 2 for node in inst1)\n",
        "\n",
        "# Creating Restricted Space constraints\n",
        "# Define the maximum number of arcs allowed as a percentage of the total number of nodes\n",
        "max_arcs = 10 #/ 100 * (2 * num_nodes - 1)\n",
        "\n",
        "# Dictionary to store constraints\n",
        "constraints = {}\n",
        "\n",
        "# Access nodes specific to 'Instance29'\n",
        "nodes_instance29 = count_arcs_instances.get('Instance130', {})\n",
        "\n",
        "# Add constraints for 'Instance29'\n",
        "instance29_variables = [variables[i, j] for (i, j) in nodes_instance29.keys()]\n",
        "if len(instance29_variables) >= max_arcs:\n",
        "    constraint = m1.addConstr(gp.quicksum(instance29_variables) >= max_arcs, name=f'Constraint_Instance130')\n",
        "    constraints['Instance130'] = constraint\n",
        "\n",
        "# Adicionando restrições de sub-rotas de Dantzig de forma iterativa\n",
        "# As restrições serão inseridas quando sub-rotas forem criadas\n",
        "def subtourelim(model, where):\n",
        "    # Importando biblioteca auxiliar\n",
        "    from itertools import combinations\n",
        "\n",
        "    if where == gp.GRB.Callback.MIPSOL:\n",
        "        # Criando uma lista de arcos selecionados\n",
        "        vals = model.cbGetSolution(model._vars)\n",
        "        selected = gp.tuplelist((i, j) for i, j in model._vars.keys()\n",
        "                             if vals[i, j] > 0.5)\n",
        "\n",
        "        # Encontrando a menor subrota criada dentro da atual lista de arcos selecionados\n",
        "        tour = subtour(selected)\n",
        "        if len(tour) < len(node_ids1):\n",
        "            # Adiciona restrição de eliminação de sub-rotas para o par de cidade na subrota\n",
        "            model.cbLazy(gp.quicksum(model._vars[i, j] for i, j in combinations(tour, 2))\n",
        "                         <= len(tour)-1)\n",
        "\n",
        "# Dada uma sublista de arcos encontre a menor subrota\n",
        "def subtour(edges):\n",
        "    unvisited = node_ids1[:]\n",
        "    cycle = node_ids1[:] # Proxy - será substituída\n",
        "    while unvisited:  # 'True' se a lista não estiver vazia\n",
        "        thiscycle = []\n",
        "        neighbors = unvisited\n",
        "        while neighbors:\n",
        "            current = neighbors[0]\n",
        "            thiscycle.append(current)\n",
        "            unvisited.remove(current)\n",
        "            neighbors = [j for i, j in edges.select(current, '*')\n",
        "                         if j in unvisited]\n",
        "        if len(thiscycle) <= len(cycle):\n",
        "            cycle = thiscycle # Nova menor subrota\n",
        "    return cycle\n",
        "\n",
        "# Resolvendo o modelo\n",
        "m1._vars = variables\n",
        "m1.Params.lazyConstraints = 1\n",
        "m1.setParam(\"Heuristics\", 0)\n",
        "\n",
        "#### Registrando tempo de início\n",
        "start_time = time.time()\n",
        "\n",
        "# Running optimization\n",
        "m1.optimize(subtourelim)\n",
        "\n",
        "#### Registrando tempo de término\n",
        "end_time = time.time()\n",
        "\n",
        "# Solução\n",
        "vals = m1.getAttr('x', variables)\n",
        "selected = gp.tuplelist((i, j) for i, j in vals.keys() if vals[i, j] > 0.5)\n",
        "\n",
        "# Rota a ser traçada\n",
        "rota1 = subtour(selected)\n",
        "\n",
        "#### Calculando tempo decorrido\n",
        "elapsed_time_DT1 = round(end_time - start_time, 4)\n",
        "\n",
        "print('________________________________________________________________')\n",
        "#### Retornando solução\n",
        "print('\\nInstância 5 - Função Objetivo: R$', round(m1.objVal, 2))\n",
        "\n",
        "#### Retornando tempo decorrido\n",
        "print(\"\\nElapsed time:\", elapsed_time_DT1, \"seconds\")\n",
        "print('\\nRota:', rota1)\n",
        "print('________________________________________________________________')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0201280e",
      "metadata": {
        "id": "0201280e"
      },
      "source": [
        "### End"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}