The thermodynamic modeling of various thermoelectric devices does not consider the influence of nonlinear Thomson effect. It leads to incomplete solutions of the heat transfer equations and serious analytical errors. On the other hand, appropriate balance among various performance parameters of thermoelectric devices is also required to improve its operating characteristics. In this context, the thermodynamic modeling on the basis of first/second laws for multielement single- and 2-stage (series/parallel) thermoelectric heat pumps, including the influence of Thomson effect in combination with Fourier heat conduction and Joule effects, is done. The optimization of the heat pumps has been carried out to obtain the optimal values of 4 input parameters by using the second version of nondominated sorting genetic algorithm in matrix laboratory. The optimal values from Pareto frontier of dual/triple objectives are obtained through 3 decision makings viz. Shannon's entropy, Fuzzy Bellman-Zadeh, and TOPSIS. It is observed that triple-objective optimization gives much lower difference between ideal and obtained solution, termed as deviation index, as compared with the single/dual one. Additionally, sensitivity analysis has been carried out to study the influence of Thomson effect on heating capacity and coefficient of performance of the pumps. To validate the evolutionary algorithm, the optimal values are compared with analytical ones from previous literature.