You can't ensure positivity and sum consistency of hierarchical forecasts if you use
What I personally found useful was to set up the summation matrix (see the original publication by Hyndman et al., 2011) and then solve the relevant (weighted) least squares problem with additional nonnegativity constraints. This will give you nonnegative and sum consistent forecasts. I have found repeatedly that this approach still results in better forecasts on all levels in the hierarchy.
This approach also allows including equality constraints, or constraints that are more general than just "$\geq 0$". I have had applications where some forecasts needed to be larger than a certain number (because of existing orders), or needed to be constrained to be equal to a given value, which you can model using two inequality constraints.
One possible tool that solves (weighted) least squares with linear constraints is the
pcls() function in the
mgcv package for R. (Note that this is slightly more general than your use case:
mgcv::pcls() allows for linear constraints, but you and the use cases I outline in the previous paragraph only need box constraints.) However, this is of course also not optimized to leverage the specific structure of potential forecast hierarchy matrices, so your performance may be significantly worse than if you use