Nonetheless, I think working with inequality constraints in models solved by perturbation methods is at least no more wrong that working with inequality constraints in log-linearized models, which is something everybody does (see the zero lower bound literature).
In both cases you're basically assuming that the polynomial approximation remains valid up to the point at which the constraint binds.
The last appendix to my paper
http://users.ox.ac.uk/~ball1377/tholden_jmp.pdf contains a particularly convenient way of generating IRFs for linear models with inequality constraints using Dynare. You basically treat the constraint binding as a news shock which arrives at the same time as the impulse, which just leaves a quadratic optimization problem to solve in matlab. If you want to simulate models with inequality constraints though, or you want higher order approximations, you won't be able to rely on Dynare to do most of the work for you though. (See the zero lower bound literature instead.)