  We develop a procedure which systematically generates all conserved operators
in the disordered models of interacting fermions. Among these operators, we
identify and count the independent and local integrals of motion (LIOM) which
represent the hallmark of the many-body localization (MBL). The method is
tested first on the prototype disordered chain of interacting spinless
fermions. As expected for full MBL, we find for large enough disorder
$N_M=2^M-1$ independent and quasi-local LIOM with support on $M$ consecutive
sites. On the other hand, the study of the disordered Hubbard chain reveals
that $3^M-1< N_M < 4^M/2$ which is less than required for full MBL but much
more than in the case of spinless fermions.
