After a bunch of trial and error, i've found a good solution (at least for my 50-node solution. I'll list all the changes I did:

1. I've gone back and simplified my evaluation function, now it takes the value of an edge, without penalization, as I found out my penalization was causing too much of an effect on the result.
2. I've added local search (memetization) on the algorithm used.
3. I've decreased the HallOfFame size to 1.
4. I added stagnation detection. If the algorithm doesn't improve after 25 generations, I increase the mutation rate. After 5 increases, I regenerate the population using the best individual and reset the mutation rate to the original. After 3 cycles of this, the genetic algorithm halts.

After this improvements, I finally got the following graph:

Looking at the evolution graph, you can clearly see the population resets that lead to a better path.

Any further tips will be greatly appreciated! :)

Final code:

# Previous code
    def eaSimpleWithElitism(self,
                            population,
                            toolbox,
                            cxpb,
                            mutpb,
                            ngen,
                            stats=None,
                            halloffame=None,
                            verbose=__debug__):
        """This algorithm is similar to DEAP eaSimple() algorithm, with the modification that
        halloffame is used to implement an elitism mechanism. The individuals contained in the
        halloffame are directly injected into the next generation and are not subject to the
        genetic operators of selection, crossover and mutation.
        """
        logbook = tools.Logbook()
        logbook.header = ['gen', 'nevals'] + (stats.fields if stats else [])

        # Evaluate the individuals with an invalid fitness
        invalid_ind = [ind for ind in population if not ind.fitness.valid]
        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit

        if halloffame is None:
            raise ValueError("halloffame parameter must not be empty!")

        halloffame.update(population)
        hof_size = len(halloffame.items) if halloffame.items else 0

        record = stats.compile(population) if stats else {}
        logbook.record(gen=0, nevals=len(invalid_ind), **record)
        if verbose:
            print(logbook.stream)

        best_val = float('inf')
        gens_stagnated = 0
        mut_exploder = 1
        cicles = 0

        # Begin the generational process
        for gen in range(1, ngen + 1):

            # Select the next generation individuals
            offspring = toolbox.select(population, len(population) - hof_size)

            # Vary the pool of individuals
            offspring = algorithms.varAnd(offspring, toolbox, cxpb, mutpb)

            # Evaluate the individuals with an invalid fitness
            invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
            fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
            for ind, fit in zip(invalid_ind, fitnesses):
                ind.fitness.values = fit

            elite = halloffame.items
            for i, e in enumerate(elite):
                ie = self.local_search_2opt(e)
                e[:] = ie[:]
                e.fitness.values = self.evaluate_tsp(e)

            # add the best back to population:
            offspring.extend(elite)

            # Update the hall of fame with the generated individuals
            halloffame.update(offspring)

            # Replace the current population by the offspring
            population[:] = offspring

            # Append the current generation statistics to the logbook
            record = stats.compile(population) if stats else {}
            logbook.record(gen=gen, nevals=len(invalid_ind), **record)
            if verbose:
                print(logbook.stream)

            val = halloffame[0].fitness.values[0]
            if val < best_val:
                best_val = val
                gens_stagnated = 0
            else:
                gens_stagnated += 1
            if gens_stagnated >= 25:
                print("Stagnated")
                if mut_exploder < 5:
                    toolbox.register("mutate", 
                                     tools.mutShuffleIndexes, 
                                     indpb=1/(self.graph.nodes - mut_exploder))
                    mut_exploder += 1
                else:
                    print("Reseting...")
                    for i, ind in enumerate(population):
                        population[i] = halloffame.items[0]
                    mut_exploder = 1
                    toolbox.register("mutate", 
                                     tools.mutShuffleIndexes, 
                                     indpb=1/(self.graph.nodes))
                    cicles += 1
                gens_stagnated = 0
            
            if cicles >= 3: break

        return population, logbook

    def run_ga_tsp(self,
                   ngen: int = 3000,
                   cxpb: float = 0.7,
                   mutpb: float = 0.2,
                   pop_size: int = 1000,
                   dir: str | None = None,
                   idx: int = 0,
                   vrb: bool = True) -> tuple[list[int], float]:
        """Runs the Genetic Algorithm for the Traveling Salesman Problem.
        
        This function calls the wrapper functions that define the creator, 
        toolbox and the attributes for the Genetic Algorithm designed to solve
        the Traveling Salesman Problem. It then runs the Genetic Algorithm and 
        returns the best path found and its total value, while also calling the
        wrapper function to plot the results.

        Args:
            ngen (optional): The number of generations. Defaults to 100.
            cxpb (optional): The mating probability. Defaults to 0.9.
            mutpb (optional): The mutation probability. Defaults to 0.1.
            pop_size (optional): The size of the population. Defaults to 200.
            dir (optional): The directory where the plots should be saved. 
                Defaults to None, in which case the plot(s) won't be saved.
            idx (optional): The index for the plot to save. Defaults to 0.
            vrb: (optional): Run the algorithm in verbose or non-verbose mode.
                Defaults to True.

        Returns:
            A tuple containing the best path found and its total value.
        """
        random.seed(169)
        if not self.graph.distances: self.graph.set_distance_matrix()
        
        creator = self._define_creator()
        toolbox = self._define_toolbox()
        population, stats, hof, = self._define_ga(toolbox, pop_size)

        population, logbook = self.eaSimpleWithElitism(population,
                                                       toolbox,
                                                       cxpb=cxpb,
                                                       mutpb=mutpb,
                                                       ngen=ngen,
                                                       stats=stats,
                                                       halloffame=hof,
                                                       verbose=vrb)

        best = [i for i in hof.items[0]]
        best += [best[0]]
        total_value = self.evaluate_tsp(best)[0]

        if vrb:
            print("-- Best Ever Individual = ", best)
            print("-- Best Ever Fitness = ", total_value)

        if dir:
            self._plot_ga_results(best, logbook, dir, idx)
        else:
            self._plot_ga_results(best, logbook).show()

        return best, total_value

    def _define_toolbox(self) -> base.Toolbox:
        """Defines a deap toolbox for the genetic algorithms.
        
        The ``deap.base.createor`` module is part of the DEAP framework. It's 
        used as a container for functions, and enables the creation of new
        operators by customizing existing ones. This function extracts the
        ``toolbox`` instantiation from the ``run_ga_tsp`` function so the code
        is easier to read and follow. 
        
        In the ``toolbox`` object is where the functions used by the genetic
        algorithm are defined, such as the evaluation, selection, crossover
        and mutation functions.

        Returns:
            The toolbox defined for the genetic algorithm.
        """
        toolbox = base.Toolbox()
        toolbox.register("random_order", 
                         random.sample, 
                         range(self.graph.nodes), 
                         self.graph.nodes)
        toolbox.register("individual_creator", tools.initIterate,
                         creator.Individual, toolbox.random_order)
        toolbox.register("population_creator", tools.initRepeat, list,
                         toolbox.individual_creator)

        toolbox.register("evaluate", self.evaluate_tsp)
        toolbox.register("select", tools.selTournament, tournsize=2)

        toolbox.register("mate", tools.cxOrdered)
        toolbox.register("mutate",
                         tools.mutShuffleIndexes,
                         indpb=1.0/self.graph.nodes)
        toolbox.register("clone", self._clone)

        return toolbox
# Rest of code