Is it possible to find the ground state of generalized Ising models? Is there a general solver (or a theoretical algorithm) for obtaining the ground state configuration of the extended Ising model, which involves an arbitrary lattice, arbitrary coordination number (i.e. $n$-body interactions for arbitrary $n$), arbitrary Hamiltonian, and arbitrary spin of each particle (-1,0,1 and so forth)? Except exact enumeration, what could we do better?
I know such a requirement is too stringent. But is there some research or established knowledge that aims to accomplish this?
 A: Yes there is. Such Ising models are called 'Potts models' (due to the arbitrary number of possible spins). This can be mapped exactly to a graph or network, for which partitioning algorithm or 'community detection' algorithms are useful. For appropriate choices of objective function, the optimal partition will correspond to the ground state of the Potts model. See, for example,

Global disorder transition in the community structure of large-q Potts systems. Peter Ronhovde, Dandan Hu and Zohar Nussinov. EPL 99 (2012) 38006. arXiv:1204.3649.

A: As in your question the stress was on the word general, I have some bad news:
an efficient "general solver (or a theoretical algorithm) for (...) extended Ising models, which involves an arbitrary lattice" does not exists. 
Of course, one can invent algorithms that, in principle, could find the ground state. The most trivial would be checking the energy of all configurations (I guess this is what you referred to by exact enumeration). However, the time used by this algorithm would scale exponentially with the number of sites. One may ask - as you did - whether one can do something substantially better, e.g., finding an algorithm where the "used time" scales polynomially with the number of sites. Unfortunately, there is no such algorithm (if we assume that NP$\ne$P). Already for the ordinary Ising Hamiltonian
$$H= \sum_{ab} J_{ab}S_aS_b $$
with $J_{ab}\in\{+1,-1,0\}$ and with connectivity on an $L\times L\times 2$ cubic lattice, it was shown that finding the ground state is an NP hard problem. The proof of this can be found here: 

F. Barahona, On the computational complexity of Ising spin glass models, J. Phys. A: Math. Gen., 15  3241 (1982). 

Of course, if you don't consider the general case, but restrict your attention to a set of "easier" lattices or graphs (e.g. to planar graphs), then there could be a polynomial time algorithms (depending on the structure of the specific restricted cases).
