I have an MILP problem with an objective of the type $\max z(x,g(x))= X - \varepsilon\cdot g(X)$, where $\varepsilon$ is a penalty term determined through an iterative process. What I want is to find for a given $z(x,g(x))$ the optimal solution (out of all the optimal solutions) that maximizes $X$ (without interfering with $\varepsilon$). All the above to be implemented in GUROBI using python.
I'm not a Gurobi user, but what you want to look at is their "solution pool" feature. If I understand the question, you want to set the solution pool to find as many solutions as possible within a tolerance of zero (or perhaps something close to zero) of the optimal value, subject to time limits etc. This is using $z(x,g(x))$ as the objective function. Once Gurobi is done doing its thing, retrieve the solutions from the pool and retain only the solution(s) with maximal $x$ value.
Assuming your iterative process happened prior to formulating the problem, simply maximizing it with Gurobi (assuming that Gurobi can solve it to global optimality) will deterministically provide the $X$ point that provides the best value for your objective, since $\epsilon$ is not a variable in your problem.
This seems a bit too straightforward though, so maybe I'm missing some nuance in your question that doesn't quite come through.