This study utilizes the genetic algorithm (GA) and Levenberg–Marquardt (L–M) algorithm to optimize the parameter acquisition process for two commonly used viscoelastic models: 2S2P1D and Havriliak–Negami (H–N). The effects of the various combinations of the optimization algorithms on the accuracy of the parameter acquisition in these two constitutive equations are investigated. Furthermore, the applicability of the GA among different viscoelastic constitutive models is analyzed and summarized. The results indicate that the GA can ensure a correlation coefficient of 0.99 between the fitting result and the experimental data of the 2S2P1D model parameters, and it is further proved that the fitting accuracy can be achieved through the secondary optimization via the L–M algorithm. Since the H–N model involves fractional power functions, high-precision fitting by directly fitting the parameters to experimental data is challenging. This study proposes an improved semi-analytical method that first fits the Cole–Cole curve of the H–N model, followed by optimizing the parameters of the H–N model using the GA. The correlation coefficient of the fitting result can be improved to over 0.98. This study also reveals a close relationship between the optimization of the H–N model and the discreteness and overlap of experimental data, which may be attributed to the inclusion of fractional power functions in the H–N model.