Abstract: The continuous-time rational approximation (CRA) function is an important numerical model to describe the dynamic impedance function of foundations in the time domain for the time history analysis of soil-structure interaction systems. The stability and accuracy of the identified CRA function determine the analysis stability and accuracy of the soil-structure interaction system in the time domain. So far, the reported identification methods for the CRA function cannot guarantee accuracy, stability and calculation efficiency simultaneously. In view of the linear system control theory, the CRA function can be decomposed as the combination of a series of first- and second-order subsystems. Furthermore, the stable boundaries of these identified parameters are built theoretically according to the stability principles of these subsystems. Based on this concept, a new parameter identification method was proposed by using the genetic algorithm and sequential quadratic programming for the CRA function. The case studies for different frequency responses show that the proposed method provides the same high accuracy with existing methods when the function is relatively simple. Moreover, the proposed method ensures the stability and also saves 35% the time cost of existing methods. For complex functions, the proposed method has more than 95% accuracy for higher order rational functions with 25% the time cost of the existing methods. The simulation results of these cases prove that the proposed method in this work not only ensures the stability of CRA function, but also enhances the identification accuracy and efficiency. Therefore, the proposed method extends the applicability of the CRA function for the time history analysis of soil-structure interaction systems.