In this paper the quasi-static initial and boundary value problem for an elasto-plastic mixed hardening material is reformulated within the constitutive framework of small strains. The plastic factor plays the basic role in describing the rate independent evolution equations for the plastic strain and hardening variables. The plastic factor is equivalently represented as the solution of an appropriate local inequality involving the yield function. The main idea was to introduce the variational inequality at any time t to be solved for the velocity field and the complementary plastic factor. There is the plastic factor in a strain-driven process. The solution procedure proposed here to solve the initial and boundary value problem is based on the solutions of the variational inequality at time t, coupled with an update algorithm in order to evaluate the current state of the material for an incremental deformation process. This time the return mapping algorithm is avoided as the values of the plastic factor and the velocity are known at time t. As we developed a procedure to simultaneously solve the equilibrium equation coupled with the rate-independent evolution equations, no necessity to compute the algorithmic elasto-plastic tangent moduli occurs. The numerical simulations are done for the mixed hardening elasto-plastic model involving Armstrong–Frederick kinematic hardening. To validate the proposed numerical algorithms, we compare the solutions based on the variational inequality and those based on return mapping algorithm, computed for the same Prager kinematic hardening law.