A data-driven modeling method with feature selection capability is proposed for the combustion process of a station boiler under multi-working conditions to derive a nonlinear optimization model for the boiler combustion efficiency under various working conditions. In this approach, the principal component analysis method is employed to reconstruct new variables as the input of the predictive model, reduce the over-fitting of data and improve modeling accuracy. Then, a k-nearest neighbors algorithm is used to classify the samples to distinguish the data by the different operating conditions. Based on the classified data, a least square support vector machine optimized by the differential evolution algorithm is established. Based on the boiler key parameter model, the proposed model attempts to maximize the combustion efficiency under the boiler load constraints, the nitrogen oxide (NOx) emissions constraints and the boundary constraints. The experimental results based on the actual production data, as well as the comparative analysis demonstrate: (1) The predictive model can accurately predict the boiler key parameters and meet the demands of boiler combustion process control and optimization; (2) The model predictive control algorithm can effectively control the boiler combustion efficiency, the average errors of simulation are less than 5%. The proposed model predictive control method can improve the quality of production, reduce energy consumption, and lay the foundation for enterprises to achieve high efficiency and low emission.