A mixed-integer linear programming (MILP) formulation is presented for parameter estimation of the Potts model. Two algorithms are developed; the first method estimates the parameters such that the set of ground states replicate the user-prescribed data set; the second method allows the user to prescribe the ground states multiplicity. In both instances, the optimization process ensures that the bandgap is maximized. Consequently, the model parameter efficiently describes the user data for a broad range of temperatures. This is useful in the development of energy-based graph models to be simulated on Quantum annealing hardware where the exact simulation temperature is unknown. Computationally, the memory requirement in this method grows exponentially with the graph size. Therefore, this method can only be practically applied to small graphs. Such applications include learning of small generative classifiers and spin-lattice model with energy described by Ising hamiltonian. Learning large data sets poses no extra cost to this method; however, applications involving the learning of high dimensional data are out of scope.