## 1 Introduction
[PERSON] and [PERSON] [29] introduced new effective values into the yield condition related to the structure and interactions of dislocations as a major plasticity mechanism. [PERSON] [30] expressed the egg-shaped yield function by decomposing deviatoric stress into its collinear and orthogonal parts with respect to the back stress of kinematic hardening. [PERSON] et al. [31] developed a non-associated elastoplastic anisotropic hardening model, which is coupled with a damage model based on [PERSON]'s approach [30]. [PERSON] et al. [32] suggested a model to represent the Bauschinger effect with the kinematic hardening component, while other strain-path change effects could be expressed through the distortion of the yield surface. Besides the kinematic hardening or combined kinematic-distortional hardening, several anisotropic hardening models are solely based on the distortional hardening approach. Some of the distortional hardening models were developed to express the yield surface evolution by using the anisotropic coefficients as a function of the plastic work or equivalent plastic strain [33; 34; 35; 36]. However, these models do not take into account the loading path change effect. [PERSON] and one the of co-authors of the present study [37; 38; 39] proposed a series of anisotropic hardening models without yield surface translation, which they named the homogeneous yield function based anisotropic hardening (HAH) models. The main concept of the original HAH model consists of using the distortion of the yield surface along a designated loading path, which is called a microstructure deviator, and the distortional amounts are controlled via adequately defined plastic state variables. Later, the model was extended for reproducing latent hardening, work hardening stagnation [38], as well as cross-hardening and softening under more general loading conditions [39]. This improved model was named the extended HAH model (or e-HAH, in this study). The performance of the distortional hardening models was validated through many applications, including the springback in U-draw bending and industrial S-rail forming [40; 41; 42; 43; 44; 45; 46]. The model is also applicable to subjects for which non-linear strain path effects are important, such as the fracture behaviors of metal after pre-deformation [47; 48; 49]. As the constitutive models of plasticity have advanced, especially the hardening laws, their implementation into FE simulations have become increasingly challenging. This is mainly attributed to the advanced constitutive models having more state variables, which are non-linearly cross-related. These challenges have led to more robust numerical formulations and implementations of the constitutive models, which eventually determine the accuracy and robustness of FE simulations. Numerous numerical algorithms were proposed to take into account the stress integration or stress update, using the elastic-plasticity constitutive models. The details of the general theoretical studies on these stress-integration algorithms for nonlinear plasticity are well documented in a book by [PERSON] [50]. The basic principle of the stress update algorithm for the classical rate-independent elastoplastic model consists of locating the stress state on the yield surface described in the six-dimensional stress space, which is consistent with the material hardening. The hardening of the material is often represented by a uniaxial stress-plastic-strain curve as a reference stress state. Most of the stress update algorithms were developed in the elastic-predictor and plastic-corrector schemes. For example, the closest point projection method (CPPM) [51; 52; 53] and the cutting-plane method (CPM) [54] are employed as popular stress integration algorithms in FE models. The CPPM is often based on the Euler backward method (EBM), and thus it is an implicit method requiring the first and second derivatives of the yield surface to satisfy both the consistency and flow rule (or normality rule).