## 1 Introduction  
The use of numerical codes that account for plasticity is essential when designing geotechnical structures, as they are functional decision-making tools. However, the implementation of a good plasticity model requires the development of powerful algorithms to resolve the numerical difficulties arising from complex plasticity laws such as those used for soil modeling. For example, [PERSON] et al. [1] present several algorithms for solving plasticity problems (that are incorporated into the HYPLAS calculation code written in FORTRAN 77). The classical Newton-Raphson scheme is very often used to solve the nonlinearities of the plasticity models since it converges at the second order, but under the condition that the initial solution is close to the true solution. Its implementation requires a rewriting of the load function as a function of the plastic deformation increment. Furthermore, it should be ensured that during the iterative procedure, this plastic increment is not less than zero, which could happen for certain plasticity models. However, the Newton-Raphson scheme can be improved by using the Homotopy-Newton-CPPM or CG-Newton-CPPM algorithms proposed by [PERSON] et al. [2]. [PERSON] et al. [3] proposed a line-search algorithm with the Newton-Raphson method making it possible to force the convergence when the first guessed solution is far from the sought numerical solution. Similarly, [PERSON] and [PERSON] [4] described the theoretical foundations of several formulations and algorithms for plasticity based on the Newton-Raphson procedure. In particular, they show that the construction of the elastoplastic tangent matrix, which is closely related to the rate of convergence in the Finite Element Method (FEM) and guarantees the consistency of the algorithmic integration process [5], is a crucial step for the robustness of these algorithms, hence the difficulty of implementing certain constitutive laws (such as the Hardening Soil Model law [6]). However, for the perfectly plastic behavior laws, particularly [PERSON]'s law, some authors have been able to establish this matrix analytically, as with [PERSON] et al. [1], who determined this matrix in the principal basis before returning it to the global basis (Cartesian) using spectral projection. Indeed, the arising functions and their derivatives expressed in the global basis are complicated, and so the eigenvectors must therefore be determined numerically [7]. This difficulty can be avoided by keeping the stresses in the main base when the plasticity or failure criterion has a very simple shape compared to its formulation in the general stress space [8]. Another approach was proposed in [9], where the authors preferred to express the yield function in a non-invariant form. This enables the derivatives and the plastic matrix to be expressed in a fairly simple algebraic form. In general, great care should be taken to avoid errors in deriving the consistent elastoplastic tangent analytically, as the calculations can become very complicated. Recently, another approach for solving local constitutive equations using conic optimization is discussed [10]. This method avoids blind guessing in the search for the elastic stress tensor (trial stresses) for the constitutive scheme. Indeed, the solution of the problem can be searched as the optimal pair (isotropic and linear behavior and linear hardening) for the convex optimization problem. Moreover, this procedure also allows to lighten the construction of the tangent operator. This operator can also be decomposed into elastoplastic and elastic tangent operators [11]. Commercial finite element software systems are available to solve geotechnical problems, including PLAXIS, ZSOIL and FLAC. However, coupling them with calibration or uncertainty analysis algorithms can be difficult due to the impossibility of accessing source codes.