A computational fluid dynamics (CFD)-based mathematical model is applied to study emission formation in an industrial-scale bubbling fluidized bed boiler burning a fuel mixture consisting mainly of biomass. A modification of the standard k-epsilon model is used to model turbulence, global mechanisms are used to model chemical reactions, the eddy-dissipation combustion model and the eddy-dissipation concept are used to model turbulence-chemistry interaction, and the finite-volume method (discrete ordinates) together with the weighted-sum-of-grey-gases model are used to model radiative heat transfer. A simplified approach is used to take into account the bubbling bed in the overall CFD model. Experimental data available on the temperature and the concentrations of NO and NH3 are presented to validate the predictions to a certain extent. It is emphasized that the CFD-based study of a real-world industrial boiler contains significant uncertainties and that the results must be evaluated with great care. The error caused by the computational grid may easily dominate the errors caused by simplified models used in industrial-scale boiler modeling. It is shown that over two million computational cells may be required for the grid-independent solution for only a single jet. On the other hand, local grid refinement is shown to be a good option for improving the predictions without an excessive increase in the number of computational cells. The predictions are shown to be strongly dependent on the chosen ammonia chemistry mechanism. Mechanisms suitable for modeling ammonia injection, i.e. the selective noncatalytic reduction (SNCR) process, are identified, and it is shown that the correct emission formation trends can be predicted as a function of the SNCR process load. Due to the different conditions prevailing in the lower and upper part of the boiler, a combination of more than one mechanism may be a good approach in modeling the whole boiler. After the CFD-based model has been validated to a certain extent, it is used for optimization. There are nine design variables (nine distinct NH3 injections in the SNCR process) and two conflicting objective functions (minimize NO and NH3 emission in flue gas). The multiobjective optimization problem is solved using the reference-point method involving an achievement scalarizing function. An interactive approach is used to incorporate the preferences of decision makers and to generate Pareto optimal solutions. Two inherently different optimization algorithms, viz. a genetic algorithm and Powell's conjugate-direction method, are applied in the solution of the resulting optimization problem. The constraints are incorporated using a penalty function. With respect to the current operating point, approximately 12% decrease in the NO emission is obtained, while maintaining the NH3 emission at an acceptable level. It is shown that the optimization connected with CFD is a promising and scientific approach for combustion optimization. The strengths and weaknesses of the proposed approach and of the interactive multiobjective optimization method used are discussed from the point of view of the complex real-world optimization problem. Also, the performance of the two optimization algorithms is evaluated. Some techniques related to the reduction of overall computational cost are applied in the present study, and some other possible techniques are reviewed.