Fictitious domain methods enable to solve physical problems on unfitted grids, thereby avoiding time-consuming and error prone meshing phases. However, an accurate integration of the weak formulation is still mandatory, leading to the need for efficient quadrature strategies in the elements that are partially located in the physical domain. Various methodologies have been proposed to this end. Among them, the design of element-specific moment-fitting quadrature rules seems promising. However, some of the resulting weights are usually negative and the points located out of the domain which may lead to a lack quadrature stability. In this contribution, we aim to construct quadrature rules whose weights are all positive and points all located in the physical domain. Such rules can be constructed based on a non-negative least square resolution of the moment-fitting equations. The resulting quadrature formulas are compared to empirical quadrature rules and different variants of classical moment-fitting quadrature rules in both 1D and 2D. Benchmarks show that non-negative moment-fitting quadrature rules are robust and efficient although their setup cost can be significantly higher than classical moment-fitting. Finally, their application to linear elastic and small-strain elasto-plastic problems highlight their robustness for engineering applications.