This paper presents an aerodynamic shape optimization methodology for rotor blades of horizontal axis wind turbines. Genetic Algorithm and Blade Element Momentum (BEM) Theory are implemented for maximization of the power production at a given wind speed, rotor speed and rotor diameter. The potential flow solver with a boundary layer model, XFOIL, provides sectional aerodynamic loads. Optimization variables are selected as the sectional chord length, the sectional twist and the blade profiles at root, mid and tip regions of the blade. The blade sections may be defined by the NACA four digit airfoil series or by arbitrary airfoil profiles defined by a Bezier curve. Firstly, validation studies are performed for the method developed. Then, several optimization studies are performed on the Risoe wind turbine. Finally, a design optimization for a 1 MW wind turbine is performed.