We develop a model for wind farm power optimization while considering the wake interaction among wind turbines. The proposed model is a Navier-Stokes equations-constrained optimization model with the objective of maximizing the total power where the operating points of the turbines are the decision variables, and the three-dimensional vorticity-velocity Navier-Stokes equations of wind speed are among the constraints. Moreover, we develop an efficient numerical algorithm to solve the optimization model. This algorithm is based on the pattern search method, the actuator line method and a time-stepping scheme which is used to solve the vorticity-velocity Navier-Stokes equations. In the numerical experiments, we first compute the power generation of a commercial wind turbine called WindSpot for different wind speed. It is shown that the computed power is in a good agreement with the measurements. Then, in the case of two turbines, we find that by optimizing the turbines’ operation while considering the wake effect, we can gain an additional 8.11% in the total power when the incoming wind speed on the boundary is 10 m/s.