Suitable conditions for thermal extraction of semi-volatile organic compounds have largely been arrived at by univariate optimization or based on the recommendations provided by the manufacturers of the extraction equipment. Herein, we demonstrated the multivariate optimization of Tenax TA-thermal extraction for determining organophosphate esters in the gas phase fraction of air samples. Screening and refining experiments were performed using the eighth fraction factorial and Box-Behnken designs, respectively, and satisfactory models were obtained. Subsequently, the process was optimized by Derringer's desirability function and the global desirability was 0.7299. Following optimization, the analytes were desorbed at 290 °C for 10 minutes at a helium flow of 95 mL min-1, with the transfer line set at 290 °C. The analytes were then cryofocused at 20 °C and then cryodesorbed into the chromatographic column at 295 °C for 6 minutes. Method validation exhibited high linearity coefficients (>0.99), good precision (CV < 14%) and low detection limits (0.1-0.5 ng m-3). The method was tested by pumping 0.024 m3 of real indoor environment air through Tenax TA sorbent tubes. Furthermore, with multivariate optimization, analysis time and other resources were significantly reduced, and information about experimental factor interaction effects was investigated, as compared to the univariate optimization and other traditional methods.