This paper addresses the problem of finding a stationary point of a nonlinear dynamical system whose state variables are under inequality constraints. Systems of this type often arise from the discretization of PDEs that model physical phenomena (e.g., fluid dynamics) in which the state variables are under realizability constraints (e.g., positive pressure and density). We start from the popular pseudo-transient continuation method and augment it with nonlinear inequality constraints. The constraint handling technique does not help in situations where no steady-state solution exists, for example, because of an under-resolved discretization of PDEs. However, an often overlooked situation is one in which the steady-state solution exists but cannot be reached by the solver, which typically fails because of the violation of constraints, that is, a non-physical state error during state iterations. This is the shortcoming that we address by incorporating physical realizability constraints into the solution path from the initial condition to steady state. Although we focus on the DG method applied to fluid dynamics, our technique relies only on implicit time marching and hence can be extended to other spatial discretizations and other physics problems. We analyze the sensitivity of the method to a range of input parameters and present results for compressible turbulent flows that show that the constrained method is significantly more robust than a standard unconstrained method while on par in terms of cost. Copyright (c) 2015 John Wiley & Sons, Ltd.