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 an interaction-matrix-based form for latent hardening. The algorithmic formulation and numerical implementation treats the displacement and glide-system slips as the primary field variables. Example 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 compared and correlated features of the macroscopic stress-strain response such as stress relaxation.