To evaluate the grand canonical partition function, as is done for no interaction, you need the microstates $E_i$. In the noninteracting case the energy of each microstate is simply the sum of the single particle energies, i.e. you could write $E_i=\sum_j \epsilon_j n_{i_j}$, with $n_{i_j}$ the occupation of the single particle state $\epsilon_j$ in the microstate $i$.
This is decomposition of many-body eigenstates into products of single particle eigenstates is no longer accurate when the interaction is nonzero. There is no way around finding all the exact eigenenergies $E_i$ of the many-body problem. Needless to say, this is extremely hard.
Perturbation theory is therefore often used to approximate
the energy of each microstate
the spectrum itself, i.e. the types of possible excitations
The procedure is: compute a many-body ground state, do linear response on it, aka Bogolyubov, evaluate the partition function and hope for the best.
Caveat: There are two main problems with this approach:
related to 1: the higher you go in energy, the worse your linear response approximation for the $E_i$ becomes.
related to 2: from exactly solvable models it is known that Bogolyubov excitations constitute only a fraction of the full excitation spectrum, see e.g. Exact Analysis of an Interacting Bose Gas. II. The Excitation Spectrum. So you do not evaluate (not even approximately) the full partition function when you rely on Bogolyubov theory. See also the plots of the excitation spectrum on Scholarpedia.
Nevertheless, if you are still interested, this question
Partition function of weakly interacting Bose gas discusses a Bogolyubov treatment.