Obtaining wrong solutions according to the beta threshold in the search of hyper rectangular subspaces

Dear Gurobi team,

I'm solving a MIP for finding hyper rectangular subspaces with accuracy over the threshold (beta) and maximum number of points. Varying beta I obtain sometimes wrong solutions,  even if optimal solution was found: in these cases I have the notice "No other solutions better than .... ".

Here is the code for the model:, where c is the result of the classification task, v and p the final indicators of the instances in the range and s outside the range, u and l the final bounds.

Any advice is greatly appreciated.

Thanks in advance!

self.c = np.asarray(c)

self.threshold_constraint = self._mdl.addConstr(
    gp.quicksum([self.v[i] * self.c[i] for i in range(n_of_samples)])
    >= beta * self.v.sum(),
    name="threshold_constraint",
)

self.bigM_lower = self._mdl.addConstrs(
    (
        -M[k] * (1 - self.p[i, k]) + self.l[k] <= x[i][k] * self.p[i, k]
        for i in range(n_of_samples)
        for k in range(n_of_features)
    ),
    name="bigM_lower-",
)

self.bigM_upper = self._mdl.addConstrs(
    (
        M[k] * (1 - self.p[i, k]) + self.u[k] >= x[i][k] * self.p[i, k]
        for i in range(n_of_samples)
        for k in range(n_of_features)
    ),
    name="bigM_upper-",
)

self.eps_lower = self._mdl.addConstrs(
    (
        M[k] * (self.p[i, k] + self.sup[i, k])
        - range_model_accuracy.eps * (1 - self.p[i, k])
        + self.l[k]
        >= x[i][k]
        for i in range(n_of_samples)
        for k in range(n_of_features)
    ),
    name="eps_lower-",
)

self.eps_upper = self._mdl.addConstrs(
    (
        -M[k] * (self.p[i, k] + 1 - self.sup[i, k])
        + range_model_accuracy.eps * (1 - self.p[i, k])
        + self.u[k]
        <= x[i][k]
        for i in range(n_of_samples)
        for k in range(n_of_features)
    ),
    name="eps_upper-",
)

self.v_upper = self._mdl.addConstrs(
    (
        self.v[i] * n_of_features
        <= gp.quicksum(self.p[i, k] for k in range(n_of_features))
        for i in range(n_of_samples)
    ),
    name="v_upper-",
)

self.v_lower = self._mdl.addConstrs(
    (
        self.v[i]
        >= gp.quicksum(self.p[i, k] for k in range(n_of_features)) - n_of_features + 1
        for i in range(n_of_samples)
    ),
    name="v_lower-",
)

self.basic_ordering = self._mdl.addConstrs(
    (self.l[k] <= self.u[k] for k in range(n_of_features)),
    name="upper_bounds_on_upper_range-",
)

range_length = [
    abs(upper_original_value[i] - lower_original_value[i]) for i in range(n_of_features)
]
self.range_length = range_length
self.n_of_features = n_of_features
self._mdl.setObjective(
    1
    + self.v.sum()
    - gp.quicksum(
        [(self.u[i] - self.l[i]) / range_length[i] for i in range(n_of_features)]
    )
    / (1 + n_of_features),
    GRB.MAXIMIZE,
)