While I do not know how you have modelled the problem you can probably add a small penalty term of the form

\mu 1/x

(with \mu being a very small positive scalar) 

to your objective. There are also other penalty functions with a similar behavior, such as

- \mu log x

In both cases, the penalty term tends to + infinity if you approach the boundary (x=0), and it is positive for feasible solutions. Also, both functions are strictly convex in x, so they have nice theoretical properties. In standard penalty algorithms you would run your favorite algorithm for a fixed positive value of \mu, then rerun the algorithm from the solution in the previous iteration while using a value of \mu which is, let's say, divided by 10 compared to the previous iteration. Experiments may be needed to find the best value to start with AND how fast you should decrease that penalty value. 

"Another option I have been considering is, in each iteration of the optimization, to set all negative variables to zero."

----- I believe many iterative(optimization) reconstruction methods simply enforce non-negativity this way, but these methods are not Lagrangian methods. 

The attached paper(ADMM method) discusses the inclusion of nonnegativity constraint. I am not very familiar with the Lagrangian methods, so I cannot tell for now if the discussion in this paper is helpful for you.

Article A Splitting-Based Iterative Algorithm for Accelerated Statis...

You can try transformation of the variables: (x_i)^2 instead of x_i (which must satisfy x_i >= 0). No additional variables. The loss of convexity (in case of its existing) is not crucial for numerical algorithms (in case of smooth objective function)