The micromechanically based network models are known to be superior over the purely phenomenological models in the analysis of unfilled rubber. However, various factors play a crucial role in the selection of the appropriate constitutive model for the analysis of technical rubbers. These are; (i) the number of available experiments under various loading conditions, (ii) maximum stretch level expected in the critical loading scenarios of the rubber component, and (iii) percentage of fillers and additives that might distort the mechanical response from that of an ideal rubber. Moreover, the number of material parameters to be identified should be as low as possible and the parameters should be physically interpretable. 12 years have passed since the celebrated work of Marckmann and Verron [Rubb. Chem. Tech. 79.5, 835-858, 2006] where they have compared 20 hyperelastic constitutive models for rubber. Since then, the rubber community witnessed emergence of many more hyperelastic models. Although the recent reviews have been very useful for the assessment of the strengths and weaknesses of various constitutive models, especially on the modeling of uncrosslinked ideal rubber, it is still a challenging endeavor for engineers to decide the best constitutive model for the specific rubber compound at hand. This paper presents a novel parameter identification toolbox based on various multi-objective optimization strategies for the selection of the best constitutive models from a given set of uniaxial tension, pure shear and equibiaxial tension experiments. The toolbox aims at providing an objective model selection procedure along with the material parameters for the rubber compound at hand. To this end, we utilize the multi-objective optimization using a successive genetic algorithm and gradient based search through the FMINCON utility in MATLAB. The novelty of our approach is (i) simultaneous fitting - the use of variable weight factors for uniaxial, equibiaxial, and pure shear data which accommodates for the (in)capability of the specific constitutive model to the specific deformation modes (ii) the sorting of the models based on an objective normalized quality of fit metric and (iii) simultaneous fittings are conducted for the models that are inherently incapable of catching large deformations, starting from zero stretch, and with 5% stretch increments until the model breaks where its quality of fit value goes above a predetermined set value. This allows identification of their validity region, i.e., the region they can work within an acceptable quality of fit margin. Accordingly, over 40 hyperelastic models are sorted with respect to the experimental data of Treloar and Kawabata.