Constrained problems are frequently encountered in classical and quantum optimization. Particle conservation, in particular, is commonly imposed when studying energy spectra of chemical and solid state systems. Though particle number-constraining techniques have been developed for fermionic (e.g. molecular electronic structure) Hamiltonians, analogous techniques are lacking for non-binary and non-fermionic problems, as in the case of bosonic systems or classical optimization problems over integer variables. Here we introduce the binary encoded multilevel particles circuit ansatz (BEMPA) – an ansatz which preserves particle count by construction – for use in quantum variational algorithms. The key insight is to build the circuit blocks by carefully positioning a set of symmetry-preserving 2- and 3-qubit gates. We numerically analyze the problem of finding the ground state eigenvalues – via the Variational Quantum Eigensolver (VQE) algorithm – of the Bose-Hubbard Hamiltonian. For a range of model parameters spanning from Mott insulator to superfluid phase, we demonstrate that our proposed circuit ansatz finds the ground state eigenvalues within drastically shorter runtimes compared to penalty-based strategies methods. Finally, we analyze the potential resource benefits of changing the qubit encoding at the end of the optimization routine. Our results attest to the efficacy of BEMPA for simulating bosonic problems for which particle number is preserved.