In this study, a two-dimensional rate-dependent gradient crystal plasticity model for non-convex energetic hardening is formulated and applied to the simulation of inelastic microstructure formation. In particular, non-convex hardening is modeled via a Landau-Devonshire potential for self-hardening and two interaction-matrix-based forms for latent hardening. The algorithmic formulation and the numerical implementation treats the displacement and the glide-system slips as the primary field variables. The numerical simulations are carried out for the case of tensile loading with periodic displacement and slip boundary conditions. The results for the formation of inelastic microstructures and their evolution under mechanical loading are illustrated together with the macroscopic stress-strain responses.