I am solving a p-median problem using genetic algorithm for facility location problem. However, after running the code I found that the fitness value come out is not stable. The code is as shown below:

# Set parameters
num_demand_points = 2130
num_stations = 46
p = 10
weights = [1 for _ in range(num_demand_points)]

# Helper function to repair an individual
def repair(individual):
    while sum(individual) < p:
        idx = random.choice([i for i, x in enumerate(individual) if x == 0])
        individual[idx] = 1
    while sum(individual) > p:
        idx = random.choice([i for i, x in enumerate(individual) if x == 1])
        individual[idx] = 0
    return individual

# Generate initial population with p stations
def generate_initial_population(size):
    population = set() 
    while len(population) < size:
        individual = [0] * num_stations
        indices = list(range(num_stations))
        random.shuffle(indices)
        selected_indices = indices[:p]
        for idx in selected_indices:
            individual[idx] = 1
        population.add(tuple(individual))
    return [list(ind) for ind in population]

# Fitness Function
def fitness(individual, distance_matrix):
    total_distance = 0
    
    for demand_point_index in range(num_demand_points):
        min_distance = float('inf')
        
        for station_index in range(num_stations):
            if individual[station_index] ==1:
                min_distance = min(min_distance, distance_matrix[demand_point_index][station_index])
        total_distance += weights[demand_point_index] * min_distance
    return total_distance

# Sort fitness
def sort_fitness(population, fitnesses):
    sorted_index = np.argsort(fitnesses)
    sorted_population = [population[i] for i in sorted_index]
    sorted_fitness = [fitnesses[i] for i in sorted_index]
    print (f"Fitness: {sorted_fitness[:5]}")
    return sorted_population, sorted_fitness

# Evaluation function
def evaluation(population):
    #Calculate fitness value for each chromosome
    pop_fitness = [fitness(individual, distance_matrix) for individual in population]
    
    #Sort population by fitness
    population_sorted, fitness_sorted = sort_fitness(population, pop_fitness)
    
    return population_sorted, fitness_sorted

# Roulette Wheel Selection
def roulette_wheel_select(population, fitnesses, selection_percentage):
    num_parents = int(len(population) * selection_percentage)
    num_parents = max(num_parents, 2)
    
    # Transform fitness values for minimization
    inverted_fitness = [1 / f for f in fitnesses]
    
    # Calculate total fitness
    total_fitness = sum(inverted_fitness)
    
    # Calculate selection probabilities
    selection_probs = [f / total_fitness for f in inverted_fitness]
    
    # Perform roulette wheel selection
    cumulative_probs = [sum(selection_probs[:i+1]) for i in range(len(selection_probs))]
   
    new_population = []
    for _ in range(num_parents):
        r = random.random()
        for i, individual in enumerate(population):
            if r <= cumulative_probs[i]:
                new_population.append(individual)
                break
                
    return new_population

# Three-parent crossover function
def three_parent_crossover(parent1, parent2, parent3):
    child = []
    for i in range(len(parent1)):
        child.append(parent1[i] if parent1[i] == parent2[i] else parent3[i])
    return repair(child)

# Binary Bit-Flipping Mutation function
def bit_flip_mutate(individual, mutation_rate=0.05):
    for i in range(len(individual)):
        if random.random() < mutation_rate:
            individual[i] = 1 - individual[i]
    return repair(individual)

# Check duplicate
def check_duplicate(list1, list2):
    for i in list2:
        if list1 == i:
            return False
    return True

def termination(no_improvement_count, generation):
    #Check if best fitness has not improved for certain number generations
    if no_improvement_count == max_no_improvement:
        return False
    if generation == max_iteration:
        return False
    return True

# Run Genetic ALgorithm
# Set Parameters
population_size = 240
max_iteration = 100
selection_percentage = 0.5
crossover_rate = 0.8
mutation_rate = 0.05
max_no_improvement = 50
elitism_count = 6
    
# Generate initial population
population = generate_initial_population(population_size)
    
# Evaluation
population_sort, fitness_sort = evaluation(population)
    
generation = 0
no_improvement_count = 0

while termination(no_improvement_count, generation):
    
    #Record best fitness before this generation
    best = fitness_sort[0]
        
    # Filter the best 50% chromosome from the population
    best_50_chromosome = population_sort[:int(population_size/2)]
    best_50_fitness = fitness_sort[:int(population_size/2)]
        
    # Create new child population
    children = []
    
    # Selection
    parent_pool = roulette_wheel_select(best_50_chromosome, best_50_fitness, selection_percentage)
        
    # Crossover
    for i in range (0, 4*len(parent_pool), 1):
        # Randomly select three distinct parents
        indices = random.sample(range(len(parent_pool)), 3)
        parent1 = parent_pool[indices[0]]
        parent2 = parent_pool[indices[1]]
        parent3 = parent_pool[indices[2]]
        
        if random.random() < crossover_rate:
            child = three_parent_crossover(parent1, parent2, parent3)
            
        else:
            child = parent1
            
        # Mutation        
        mutated_child = bit_flip_mutate(child, mutation_rate)
        
            
        # check duplicate
        if check_duplicate(mutated_child, best_50_chromosome) and check_duplicate(mutated_child, children) and len(children) < int(population_size/2):
            children.append(mutated_child)
                
    
    # Combine all population
    new_population = best_50_chromosome + children 
    
    # Elitism: carry over the top individuals
    new_population[:elitism_count] = population_sort[:elitism_count]
        
    # Evaluation
    population_sort, fitness_sort = evaluation(new_population)
    
    if len(new_population) > population_size:
        population_sort = population_sort[:population_size] 
        fitness_sort = fitness_sort[:population_size]    

    best_individual, best_fitness = population_sort[0], fitness_sort[0]
    print(f"Generation {generation} - Min: {best_fitness} - CS: {population_sort[:5]} - FS : {fitness_sort[:5]} - PL:{len(children)}")
        
    # Record improvement in best fitness
    new_best = fitness_sort[0]
    if best == new_best:
        no_improvement_count+=1
    else: 
        no_improvement_count = 0
    
    generation +=1
    
print(f"Best individual: {best_individual}, Best fitness: {best_fitness}")

What I am trying to get is the decreasing answer so that the convergence graph can look like this.
enter image description here

Besides, since I will have different mutation rate and selection operators for analysis, thus for now these are fix and should not be change. Thank you.