Since EI has an analytical formula, argmax(EI) should easily incorporate constraints. An example is weights, where sum(x's)=1, or points in a circle, where x^2+y^2<1. I've seen papers that use matlab's DIRECT algorithm.

As for python, scipy's COBYLA is a non-linearly constrained optimization algorithm. I'll try to see how it works when I have the time.