Using a penalty function to embed a non-negativity constraint can be done. You might want to take a look at e.g. Kim/Kollmann/Kim (2010) Solving the incomplete market model with aggregate uncertainty using a perturbation method) or dynare.org/wp-repo/dynarewp027.pdf