Adding criteria to objective function does not change the solution

Question

I have a MINLP problem that could be described with the following:

• a decision matrix NxM with binary variables
• a weight matrix T sized NxM
• the main criteria is to maximize sum of selected weights (the ones with 1 corresponding place in decision matrix)
• subject to hard constraint: each row in decision matrix must have exactly S entries equal to 1

The above is trivial to solve; however, I want to introduce a penalty that depends on each pair of ones in each column that "overlap". There is an overlap_matrix[N][N] that defines penalty values for each pair and the total penalty value (multiplied by penalty weight) reduces the total objective value. Since I cannot query pyomo variables, I multiply the penalty value with the model decision variable for all pairs. I have defined a model in pyomo and tried solving with different solvers (couenne, ipopt, bonmin). However, no matter what penalty factor I choose, the resulting decision matrix is always the same (as in the case with no penalty). The objective value gets lower when increasing the penalty factor, but the solution does not change. (I have verified that a better solution does exist with a metaheuristic.)

def get_penalty(m):
    total = 0
    for day in range(M):
        for i in range(N):
            penalty = 0
            for j in range(i + 1, N):
                penalty += overlap_matrix[i, j] * m.x[i, day]
            total += penalty
    return total

def get_reward(m):
    total = 0
    for day in range(M):
        reward = 0
        for i in range(N):
            reward += T[i, day] * m.x[i, day]
            for j in range(i + 1, N):
                reward -= overlap_matrix[i, j] * m.x[j, day]
        total += reward
    return total

def objective(m):
    # main objective
    total = sum(T[i, j] * m.x[i, j] for i in range(N) for j in range(M))

    # secondary: reward or penalty
    penalty = reward = 0
    penalty = get_penalty(m) * penalty_w
    reward = get_reward(m) * reward_w

    obj = total
    obj += (-penalty + reward) 
    return obj


model = pyo.ConcreteModel()
model.x = pyo.Var(range(N), range(M), within=pyo.Binary)
model.obj = pyo.Objective(rule = objective, sense = pyo.maximize)
# Constraint
def limit_constraint(model, i):
    return sum(model.x[i, j] for j in range(M)) == S
model.limit = pyo.Constraint(range(N), rule=limit_constraint)
solver = pyo.SolverFactory('couenne')
result = solver.solve(model, tee = True)

For a workround, I also tried to use a reward (instead of penalty), which is increased for a single selected variable, but decreased for each overlapping pair (in the code above). However, the same happens: the solution does not budge from the "optimal" one with no penalty or reward, although the final objective value is clearly affected. Am I doing something wrong?...