I think one way you can do this, although probably at an increased computational cost relative to the energy-penalty strategy you mentioned (which is what I would try first) is to do the following:
- define a Hamiltonian MPO (assuming this is a 1D problem) in the usual Hilbert space which includes all configurations including ones which violate your constraints
- define "isometry" tensors which map three neighboring site indices into one "mega site" index, sending to zero any states which violate your constraints
- since these isometries will overlap, try to resolve them into a MERA-like tensor network by introducing unitary two-site maps (gates) that allow you to then use two-site isometries on the indices resulting after the unitary maps
- finally, apply these MERA-like tensors to your MPO and use an algorithm like DMRG to find the ground state in your modified MPO
This is just an idea off the top of my head which I haven't worked out in any detail, so I'm not sure if step #3 is possible in the way I'm guessing it might be. It's just an idea to get started with. Hope it helps though-
Another off-the-top-of-my-head idea is to work out the most general form of MPS that obeys your constraints #1 and #2. It might have a block-sparse structure, similar to quantum number conserving MPS. Then you could formulate an MPS with such a constrained structure, write down the energy expectation value expression in terms of the MPS blocks (naming them like "A", "B", "C") and then use gradient descent or some other optimization technique to minimize the energy. Once you found the optimal blocks, you could "load" them back into your MPS to compute other properties numerically.