Penalty functions not working as expected in PyGAD

[FoxtrotF87]

I'm trying to use PyGAD to solve a minimization problem which is part of a bigger algorithm I'm trying to implement, but I have no idea on the appropiate way to include the constraints. I've tried to add if structures to the fitness function to add a penalty value to the objective, but it isn't always working.

What I'm trying to solve is a binary decision problem, which is meant to place a rectangle in a set of points that define a grid inside a bigger rectangle in which some rectangles might have been already placed. The current function just takes into account the distance to the already placed rectangles, but I'll probably add more parameters in the form of rewards and penalties. The constraints to avoid exceeding the site boundaries, the overlap between rectangles and placing the rectangle avoiding or inside some specific zones, make the problem quite complex, so thats why I'm trying to use a discrete set of possible locations and the GA to solve it.

I think I can rework some constraints, also have to add others. I tried using 1/Distance and -Distance, due to the fact that PyGAD always tries to maximize and I want to minimize said Distance. My current fitness function is the following:

def fitness_func(solution, solution_idx):
    Distance = 0
    xi = sum(LookUpList[i][0]*solution[i] for i in range(len(LookUpList)))
    yi = sum(LookUpList[i][1]*solution[i] for i in range(len(LookUpList)))
    alphai = sum(LookUpList[i][2]*solution[i] for i in range(len(LookUpList)))
    LXi = Stages[Current]['LX']
    LYi = Stages[Current]['LY']
    VXi = (LXi/2*(alphai-1)**2,LYi/2*alphai)
    VYi = (LYi/2*(alphai-1)**2,LXi/2*alphai)
    Penalty = np.inf
    
    # Only placed in one spot
    if sum(solution) > 1:
        Distance += Penalty
    
    # Site boundary constraint
    if xi+VXi[0]+VXi[1] >= Site[0] or xi-VXi[0]-VXi[1] <= 0 or yi+VYi[0]+VYi[1] >= Site[1] or yi-VYi[0]-VYi[1] <= 0:
        Distance += 1000*max(Site)
    
    # Avoid overlap between facilities
    for p in Previous:
        xp = Stages[p]['X']
        yp = Stages[p]['Y']
        alphap = Stages[p]['Alpha']
        LXp = Stages[p]['LX']
        LYp = Stages[p]['LY']
        VXp = (LXp/2*(alphap-1)**2,LYp/2*alphap)
        VYp = (LYp/2*(alphap-1)**2,LXp/2*alphap)
        if xi+VXi[0]+VXi[1] <= xp+VXp[0]+VXp[1] and xi-VXi[0]-VXi[1] >= xp-VXp[0]-VXp[1] and yi+VYi[0]+VYi[1] <= xp+VYp[0]+VYp[1] and yi-VYi[0]-VYi[1] >= xp-VYp[0]-VYp[1]:
            Distance += Penalty
    
    # Zones where a certain facility can't be placed
    for e in ExclusionZones:
        if Current == ExclusionZones[e]['Facility']:
            xp = ExclusionZones[e]['X']
            yp = ExclusionZones[e]['Y']
            LXp = ExclusionZones[e]['LX']
            LYp = ExclusionZones[e]['LY']
            if xi+VXi[0]+VXi[1] <= xp+LXp/2 and xi-VXi[0]-VXi[1] >= xp-LXp/2 and yi+VYi[0]+VYi[1] <= xp+LXp/2 and yi-VYi[0]-VYi[1] >= xp-LXp/2:
                Distance += Penalty
    
    for p in Previous:
        Distance += abs(xi-Stages[p]['X']) + abs(yi-Stages[p]['Y'])
    fitness = - Distance
    return fitness

The configuration of the GA and it's execution are done as follows:

num_generations = 150
num_parents_mating = 2

sol_per_pop = 10
num_genes = 2*(Grid[0]-1)*(Grid[1]-1)
gene_type = int

init_range_low = random_mutation_min_val = 0
init_range_high = random_mutation_max_val = 2

parent_selection_type = 'sss'
keep_parents = 1

crossover_type = 'single_point'

mutation_type = 'random'
mutation_by_replacement = True
mutation_percent_genes = 10

save_solutions = False

ga_instance = pygad.GA(num_generations = num_generations,
                       num_parents_mating = num_parents_mating,
                       fitness_func = fitness_func,
                       sol_per_pop = sol_per_pop,
                       num_genes = num_genes,
                       gene_type = gene_type,
                       init_range_low = init_range_low,
                       init_range_high = init_range_high,
                       parent_selection_type = parent_selection_type,
                       keep_parents = keep_parents,
                       crossover_type = crossover_type,
                       mutation_type = mutation_type,
                       mutation_by_replacement = mutation_by_replacement,
                       mutation_percent_genes = mutation_percent_genes,
                       save_solutions = save_solutions
                       )

ga_instance.run()

I have a function that evaluates the dictionary called Stages, which stores all the data relevant to the algorithm, that gives me the final cost value. It matches the one I get after running the PyGAD instance, but when plotting the solution with another function (I don't think is relevant, just a matplotlib figure with shapes drawn) I can see the solution isn´t always in the feasible space. I can understand some overlap due to the grid being finite so placing the new facility in one spot would be the best solution if this spot was placed a little bit lower, upper or to the side. However, this could be adjusted if the penalty function took into account how much it overlaps, so it doesn't bother me that much.

What I dont understand is why the cost function just gives me the distance, not including the penalty value that should be added as it's definitely violating the constraint stated in the condition. Should I find another way of stating constraint violation? Also, I'm starting to find that some runs don't give a solution at all.

[ahmedfgad]

So, if you have Distance=0.5 for example and one of the conditions hold, then the value returned by the fitness function is still 0.5 not Inf? Correct me if I am wrong please.