The description of the complex interplay between deformation-induced twinning and dislocation slip, typical for metals showing an hcp structure such as magnesium, is of utmost importance for understanding their deformation behavior. In the present paper, an incremental energy principle is presented for that purpose. Within this principle, dislocation slip is modeled by crystal plasticity theory, while the phase decomposition associated with twinning is considered by a mixture theory. This mixture theory naturally avoids the decomposition of the twinning process into so-called pseudo- dislocations followed by a reorientation of the total crystal. By way of contrast, the proposed model captures the transformation of the crystal lattice due to twinning in a continuous fashion by simultaneously taking dislocation slip within both, possibly co-existent, phases into account. The shear strain induced by twinning as well as the deformation history are consistently included within the twinned domain by an enhanced multiplicative decomposition of the deformation gradient. Kinematic compatibility between the different phases is enforced by a Hadamard-type compatibility condition, while compatibility with respect to the boundary conditions requires the introduction of a boundary layer. The evolution of all state variables such as the twinning volume and the plastic strains associated with dislocation slip follow jointly and conveniently from minimizing the stress power of the total crystal. This canonical variational principle is closely related to the postulate of maximum dissipation and guarantees thermodynamical consistency of the resulting model. Particularly, the second law of thermodynamics is fulfilled. In contrast to previous models suitable for the analysis of the deformation systems in magnesium, the Helmholtz energy of the twinning interfaces and that of the aforementioned boundary layer are considered. Analogously, the energy due to twinning nucleation and that related to twinning growth are accounted for by suitable dissipation functionals. By doing so, the number of twinning laminates becomes an additional unknown within the minimization principle and thus, the thickness of the lamellas can be computed. Interestingly, by interpreting this thickness as the mean free path of dislocations, a size effect of Hall–Petch-type can naturally be included within the novel model. The predictive capabilities of the resulting approach are finally demonstrated by analyzing the channel die test. For that purpose, a certain rank-two laminate structure is considered. However, it bears emphasis that the proposed framework is very general and consequently, it can also be applied to other materials.