Elastoplastic models with multiple plastic surfaces is an alternative to represent the nonlinear behavior of materials such as concrete, rocks and soils. The nonlinear response of these materials depends highly on the stress state. However, these models require the definition of many parameters which do not always have physical meaning. In addition, the implementation of elastoplastic models represented by multiple plastic surfaces brings additional complexities to the analysis. The use of this type of model requires a robust numerical integration scheme of the elastoplastic evolution equations. This work presents two contributions. The first contribution is a robust return mapping algorithm for a multisurface plasticity model in general stress space known as closest point projection algorithm. The return mapping algorithm is based on a numerical method for unconstrained optimization. In this scenario, it is adopted the Newton-Raphson method with line search. A consistent tangent modulus for multisurface plasticity is also proposed. The second contribution is a methodology for parameter calibration. This methodology is formulated as an optimization problem, with the solution obtained through a genetic algorithm. A parametric study is developed in order to better undestand specific parameters of the algorithm, solving global optimization problems. Robustness and effectiveness of the proposed algorithm are evaluated through numerical examples applied to a constitutive model used for modelling concrete, rocks and soils: Cap Model. Applications available in the literature are analysed. Lastly, the parameters of this model are calibrated using experimental data avaliable in the literature. Thus, this work aims at improving the feasibility of the use of complex elastoplastic models in engineering problems.