What are the constraints on the two electron integral matrix?

In semi-empirical methods, the two electron integral matrix is simplified and the entries are replaced by actual values instead of integrals to evaluate. However, while doing that, the variational principle does not seem to directly hold; in other words, you can put values in the two electron integral matrix that make the energy go down infinitely. To make the variational principle hold, what are the contraints that the constants you put in the two electron integral matrix obey?

The basic premise of semi-empirical methods is to replace the computationally expensive integrals with values fit to experimental data or, in some cases, to high level ab initio calculations of the relevant quantity [1]. You are completely correct to point out that making this substitution will break the theoretical validity of the variational principle.

The short answer to your main question then is that the constraints are given by the experimental values to which one fits. That is to say, if you fit the values of your integrals to experimental data and all of a sudden you get out energies which are unbounded, discontinuous, or have some other undesirable and un-physical feature, then the problem is either in the fitting or in the corollary between the experiment and integral you are fitting. So, as in many physical contexts, the answer is circular. Sure, we could have any value for the integral, but the answers which experiments give live in a very small region of all possible answers to the question.

Approximations:

First of all, there are tons of semi-empirical methods which differ slightly in how they do the fitting, what approximations are made, etc. One important thing to note is that no semi-empirical method is really claiming to get the physics of these interactions correct for the right reason. They all rely on cancellation of errors, but this is built into the method on purpose, so this cancellation tends to be rather robust. One can achieve this because the behavior of the integrals arising in electronic structure is quite systematic, so if you are careful to know what you neglected in one part, you can try to neglect something of the opposite sign in another part. Reading ref. [1] is a good place to start to see how this is done in practice.

DFTB:

I should also point out that density function tight binding (DFTB) is often described as a semi-empirical method, but DFTB does obey a variational principle[2, 3]. That is, one can expand the total energy from regular density functional theory with respect to variations of the charge density (the central object in DFT). To zeroth order, one gets a non-iterative theory which depends only on properties of the individual atoms in the molecules and then corrects for interatomic repulsions. This is trivially variational.

If one includes terms up to second-order in the density variations, then you get self-consistent-charge DFTB (SCC-DFTB). This still defines an approximate Hamiltonian, so one can variationally solve for the energy. So, you might ask why this is considered semi-empirical. I'll just quote from ref. [2] as they give a nice summary.

Comparing SCC-DFTB with MNDO-type methods, the underlying philosophy and the derivation of the applied approximations seem quite different at first sight. If one focuses on the actual working equations and their implementation, there are, however, many similarities: All these methods are valenceelectron SCF-MO treatments with a minimal valence AO basis set, only one-center and two-center terms are included, the twocenter two-electron integrals are represented by damped Coulomb interactions with the correct limit at large and small distances, and there are repulsive atom-pair terms that correct for deficiencies in the formalism. The methods are semiempirical in the sense that they employ empirical parameters which are adjusted to reproduce selected reference data, i.e., DFT energies and potential curves in the case of SCC-DFTB, and experimental heats of formation and other experimental data in the case of the MNDO-type methods.

Hopefully I've addressed your main concerns in the original question and left you with some resources so you better approach further questions.

References:

[1]: Thiel, W. (2014). Semiempirical quantum–chemical methods. Wiley Interdisciplinary Reviews: Computational Molecular Science, 4(2), 145-157.

[2]: Otte, N., Scholten, M., & Thiel, W. (2007). Looking at self-consistent-charge density functional tight binding from a semiempirical perspective. The Journal of Physical Chemistry A, 111(26), 5751-5755.

[3]: Seifert, G. (2007). Tight-binding density functional theory: an approximate Kohn− Sham DFT scheme. The Journal of Physical Chemistry A, 111(26), 5609-5613.