Important mathematical models in science and technology are based on first-order symmetric hyperbolic systems of differential equations whose solutions must satisfy certain constraints. When the models are restricted to bounded domains, the problem of well-posed, constraint-preserving boundary conditions arises naturally. However, for numerical solutions, finding such boundary conditions may represent just a step in the right direction. Including the constraints as dynamical variables of a larger, unconstrained system associated to the original one could provide better numerical results, as the constraints are kept under control during evolution. One of the main goals of this work is to investigate this idea in the case of constrained constant-coefficient first-order symmetric hyperbolic systems of differential equations subject to maximal nonnegative boundary conditions.