Minimizing energy consumption and reducing pollutant emissions during the carbon anode baking process are critically important for the aluminum industry. The present study investigates the effects of oxidizer inlet temperature, inlet oxygen concentration, equivalence ratio, refractory wall thermal conductivity, and refractory wall emissivity on the baking process using unsteady Reynolds-averaged Navier–Stokes (URANS)-based simulations in conjunction with the presumed probability density function method. Numerical results are combined with a response surface methodology (RSM) to optimize the anode baking process. The advantage of the coupled method is that it can adequately provide information on interactions of different input parameters. It is remarked that the significance level of the studied parameters varies drastically for different outputs. It is noted that diluting inlet oxygen concentration (from 23% in atmospheric air to 15%) at an elevated oxidizer temperature leads to enhanced furnace fuel efficiency, more uniform temperature distribution, and lower pollutant emissions. A linear model is detected to be adequate for response surface modeling of the anode baking furnace NOx formation. On the other hand, furnace soot formation is modeled with a higher-order model due to the quadratic behavior of the response.