The theory and methodology of a three-level orthogonal array design as a chemometric method for the optimization of analytical procedures were developed. In the theoretical section, firstly, the matrix of a three-level orthogonal array design is described and orthogonality is proved by a quadratic regression model. Next, the assignment of experiments in a three-level orthogonal array design and the use of the triangular table associated with the corresponding orthogonal array matrix are illustrated, followed by the data analysis strategy, in which significance of the different factor effects is quantitatively evaluated by the analysis of variance (ANOVA) technique and the percentage contribution method. Then, a quadratic regression equation representing the response surface is established to estimate each factor that has a significant influence. Finally, on the basis of the quadratic regression equation established, the derivative algorithm is used to find the optimum value for each variable considered. In the application section, microwave dissolution for the determination of selenium in biological samples by hydride generation atomic absorption spectrometry is employed, as a practical example, to demonstrate the application of the proposed three-level orthogonal array design in analytical chemistry.