We develop large-scale numerical optimization methods for optimal control of steady incompressible Navier-Stokes flows. The control is affected by the suction or injection of fluid on portions of the boundary, and the objective function represents the rate at which energy is dissipated in the fluid. We use reduced Hessian sequential quadratic programming methods that avoid converging the flow equations at each iteration. Both quasi-Newton and Newton variants are developed and compared to the approach of eliminating the flow equations and variables, which is effectively the generalized reduced gradient method. Optimal control problems are solved for two-dimensional flow around a cylinder and three-dimensional flow around a sphere.