The dynamics of quantum systems can be approximated by the time propagation of Gaussian wave packets. Applying a time dependent variational principle, the time evolution of the parameters of the coupled Gaussian wave packets can be calculated from a set of ordinary differential equations. Unfortunately, the set of equations is ill behaved in most practical applications, depending on the number of propagated Gaussian wave packets, and methods for regularization are needed. We present a general method for regularization based on applying adequate nonholonomic inequality constraints to the evolution of the parameters, keeping the equations of motion well behaved. The power of the method is demonstrated for a nonintegrable system with two degrees of freedom.