In the standard formulation of the population balance equation that consists of breakage terms, significant loss of mass is observed for the dispersed phase. This mass loss is actually caused by the inexact conservation property reflected by many breakage kernels; hence, incorrect physical interpretations of the model simulations may be drawn. In this work, a constrained method is developed enforcing mass conservation. This numerical property is accomplished by adding an extra restriction to the original population balance equation in terms of the dispersed phase continuity equation through the Lagrange multipliers strategy. The discretized system resulting from applying the method to a two-phase population balance equation problem is symmetric and pseudopositive definite. Numerical experiments are carried out simulating the motion of a two-phase mixture passing through a 2D domain. The results obtained by the modified least-squares spectral element method show that the mass is conserved everywhere in the domain with high accuracy.