Strongly correlated materials are now under intense development, and natural orbital functional (NOF) methods seem to be able to capture the physics of these systems. We present a benchmark based on the Hubbard model for a class of commonly used NOF approximations (also known as reduced density matrix functional approximations). Our findings highlight the importance of imposing ensemble N-representability conditions in order to obtain consistent results in systems with either weak or strong electronic correlation, such as the Hubbard system with a varying two-particle interaction parameter. Based on the accuracy of the results obtained using PNOF7, which retrieves a large amount of the total strong nondynamic correlation, the Hubbard model points out that N-representability gives solid foundations for NOF development.