Wind Turbine Blade Design Optimization

Wind Turbine Blade Design Optimization Anonymous MIT Students Massachusetts Institute of Technology, Cambridge, MA, 02139, USA We develop a methodology for analyzing wind turbine blade geometries and pitch con­ trol schemes over a range of incoming wind speeds. We use this model for an orthogonal array design of experiments, a gradient-based sequential quadratic programming optimiza­ tion, and a multi-objective genetic algorithm to maximize the expected power output while minimizing the blade volume and structural stress violations. Design of experiments gen­ erates good results with little expense. Sequential quadratic programming with Hessian re-scaling and multiple starting points generates good design vectors with a large compu­ tational expense, and heuristic algorithms such as the genetic algorithm are not suited to this problem. The best design point achieves between 60 and 70% of the Betz limit for efficiency for a large range of incoming wind speeds. Nomenclature R Qmax t k T F C β x J(x, param) PE P Q Vblades ξ(x, param) r σmax σallowable cweibull kweibull ω Ft /ΔR,Ft /ΔR v0 H Blade radius, m Maximum generator torque, N-m Blade shell thickness, m Cut-off velocity in standard deviations above the mean wind speed Twist distribution vector, radians Foil shape parameter distribution, non-dimensional Chord length distribution, m Pitch control curve, radians Design vector Objective function Expected power output, W Turbine power at a particular wind speed, W Turbine torque at a particular wind speed, N-m Blade material volume, m3 Penalty function Maximum blade stress, MPa Maximum allowable blade stress, MPa Scale parameter of the Weibull distribution of incoming wind speeds Shape parameter of the Weibull distribution of incoming wind speeds Angular velocity, rad/s Tangential and axial aerodynamic load distributions on the blade, N/m Incoming wind speed ahead of the blade, m/s Hessian matrix, matrix of second derivatives I. Motivation As renewable energies become a growing part of the energy portfolio, focus is being put on increased performance and efficiency of proven sources such as horizontal axis wind turbines. Modern commercial power wind turbines are predominantly horizontal-axis, three-bladed behemoths, with a fixed blade design ∗ Graduate Student, Computation for Design and Optimization 1 of 11 American Institute of Aeronautics and Astronautics that is adapted to varying wind conditions by a blade pitch control mechanism. They are complex systems whose design requires the integration of many engineering disciplines including aerodynamics, structures, controls, and electrical engineering. Previous wind turbine design optimization techniques have focused on specific regions of interest, including optimal control,1 optimal blade shape, and site-specific performance increases.2 Our goal was to combine all of these with a simplified analysis model to make computation tractable, but to design a blade from scratch with optimal performance over a site-specific wind profile probability distribution. II. A. Problem Formulation Goal The goal of our project was to design an electricity-generating wind turbine in a manner that maximizes power output over all expected wind conditions while minimizing construction costs. The predominant expense of wind turbine construction is the manufacture of the turbine blades.3 As such, we focus on minimizing the amount of material required to make the blade, while maximizing power output. B. System Boundary The considered disciplines include aerodynamics, structures, and control. We consider a range of incoming velocities that represent the possible operating conditions of the turbine, calculating the expected power output and extreme structural load over this range. We model the generator as a resistive torque and a constraint on the maximum torque allowed, and design a control system to keep the blades at the optimal pitch angle for a given incoming wind velocity. We do not consider tower design, nor do we consider nacelle shape. However, we assume there is no wind within 20% of the blade radius, roughly the expected nacelle size. To further limit our design space, we assumed a three-bladed design with aluminum blades. We choose three blades because an initial design of experiments (DOE) test showed that three blades outperformed both four and five-bladed designs in almost all cases. Although modern blades are often being constructed with composites, our knowledge of structural behavior in composites is limited, and we therefore elected to assume aluminum blades. Finally, we assume that our wind turbine has a control system that allows it to feather the blades in a manner that the structural stresses are alleviated when the wind speed is too high. C. Design Vector In order to perform numerical optimization, the design must be broken down into a list of decision variables and parameters that completely define the design. In order to reduce the order of our model, we used a minimal number of design variables to define the blade shape but then mapped the design vector into a higher dimensional discretization for analysis by using Piecewise Cubic Hermite Interpolating Polynomials (PCHIP). The bounds were chosen to mirror the physical constraints of the wind turbine design. For example, the blade radius is limited to 16.15 m, because this is the standard flatbed truck length of 52 feet, and each blade would need to be transported to the site. Also, the sum of the bounds for the twist angle and pitch control angles is equal to a range from zero to ninety degrees. This ensures a realistic range of feasible geometric blade angles. The design vectors and bounds are shown in Table 1. In addition, a set of parameters that include physical constants such as air density, or fixed design parameters such as the number of turbine blades, are defined in Table 2. For the twist distribution, the twist at the hub is assumed to be zero, so only two decision variables for the value of twist at the midpoint and tip of the blade specify the entire distribution. Similarly, the foil shape parameter specifies a linear combination of two airfoil shapes and characteristics. The foil at the root is the S814 foil, and foil at the tip is the S813 foil, as specified by.4 The value of the foil shape is considered to be zero at the root and one at the tip, and therefore, only a single decision variable is required to specify the parameter at the midpoint. Both the chord and control curve distributions are specified by three points. Also note that angles are specified in radians in our codes, but are displayed in degrees in figures for convenience. D. Objective The performance of the turbine can be measured as the cost of producing a unit amount of power, or inversely, the power output obtained by a given initial cost. We chose to measure the performance of the 2 of 11 American Institute of Aeronautics and Astronautics Table 1: Decision Variables Decision variable Symbol Dimension Blade radius Maximum generator torque Blade shell thickness Maximum design wind speed Twist distribution Foil shape distribution Chord length distribution Pitch control curve R Qmax t k T F C β 1 1 1 1 2 1 3 3 Lower bnd. Upper bnd. Units 5.00 1000 0.004 0 0 0 0.1 0 16.15 20000 0.020 4 π/4 1 10 π/4 meters Newton-meters meters std. dev. above mean radians non-dimensional meters radians Table 2: Selected Design Parameters Parameter Symbol Value Units Blade material yield strength Weibull distribution scale parameter Weibull distribution shape parameter Points in blade discretization Number of wind speeds in Weibull integration σY cweibull kweibull n nW 20 5 2.2 11 7 MPa ­ ­ ­ ­ turbine as the expected power output over a set of incoming wind velocities given by a Weibull distribution with site-specific shape coefficients divided by the volume of material required to construct the turbine blades Objective = J(x, param) = Expected Power Output PE (x, param) = , Blade Material Volume Vblades (x, param) (1) where x is the turbine design vector as defined in Table 1, param is a set of fixed parameters as defined in Table 2, PE is the expected power in Watts, and Vblades is the combined material volume of all blades in m3 . Together, x and param completely define a wind turbine design and all the constant parameters required to derive the expected power output, blade volume, and blade structural stresses. E. Constraints In order to make our code as general as possible, and to allow for use with heuristic algorithms, we decided to model our structural stress constraints in terms of a penalty function as opposed to including them explicitly. While logarithmic barrier functions are extremely attractive in terms of convergence properties, they are typically more useful for convex feasible sets, and it is important to change the barrier coefficient at the correct rate so that the solution converges quickly and never becomes infeasible.5 Since we have little intuition or knowledge to predict when and how the structural stresses will exceed their limits and the feasible space is likely non-convex, we elected instead to employ a square-term penalty function of the form �� � r �2 � ξ(x, param) = ρpenalty (σmax − σallowable )+ , (2) where ξ(x, param) is the penalty function of a particular design, (arg)+ means max(arg, 0), ρpenalty is the r is the maximum stress in MPa at a particular penalty coefficient applied to each blade stress violation, σmax blade section r over all operating wind speeds, and σallowable is a fraction of the yield stress, we chose 70%. Because the stress violations are all in the same units, we applied one common penalty coefficient to all blade sections. Note that when σmax < σallowable for all blade radii, ξ(x, param) = 0. For aluminum with σY = 20 MPa, σallowable = 70% × 20 MPa = 14 MPa. 3 of 11 American Institute of Aeronautics and Astronautics III. A. Methodology N2 Diagram As a concise way of summarizing our analysis code, its modules, the interdependencies between modules, and the inclusion of optimization, we present an N2 diagram of our analysis code in Figure 1. Here Q is the turbine torque, P is the power at a particular point, Ft /ΔR and Fa /ΔR are the tangential and axial load distributions on the blade, ω is the angular velocity, and v0 is the incoming wind velocity ahead of the turbine disk. An inter-module optimization is performed by the Control module, which chooses the optimum ω for a given point by calling the Aero module successively and performing a line search. The global optimization is performed by calling the Wrapper module, and the Aero module performs its own root-finding technique to approximate the wind velocity at the turbine disk. Inputs x, param x, param Wrapper Max(σmax(v0)) cweibull, kweibull Expected Power x, param x, param x, param Output PE, Vblades, Max(σmax(v0)) v0 ω, v0 ω Control ω, v0 Q, P, Ft/ΔR, Fa/ΔR Q, P, Ft/ΔR, Fa/ΔR Aero σmax(v0) ω, Ft/ΔR, Fa/ΔR Structure Figure 1: N2 Diagram. B. Wrapper Module The Wrapper module takes as inputs the design vector and parameters and outputs the expected power, blade material volume, and maximum structural stress, thereby serving as the wrapper for the entire function evaluation. It transforms the decision variables representing distributions into PCHIP splines, which it then evaluates at a given number of points, in order to discretize the blade for analysis. The Wrapper module calls the Expected Power module. C. Expected Power Module As is documented in the literature, wind profiles are commonly approximated by Weibull distributions, with data for specific sites available in6 and.2 The Expected Power module discretizes the control curve spline, β, evaluates the performance of the turbine at each incoming wind speed by calling the Control module, and performs a Simpsons Rule integration over the Weibull distribution on the turbine power output in order to derive the expected value. It also calls the Structures module to find the extreme stress value in each blade section over all incoming wind speeds. D. Control Module The Control module optimizes power output for a set of blade geometry decision variables and a given control pitch angle by choosing the optimal rotational speed without violating the maximum allowable torque, Qmax, constraint on the generator. The algorithm assumes a concave power curve with respect to rotational speed, ω, and finds the point of maximum power by MATLAB's fminbnd routine. If that operating point exceeds Qmax, the rotational speed is incrementally increased until the operating torque matches Qmax by employing MATLAB's fzero routine. E. Aero Module The Aero module calculates the turbine torque and power for a given rotational speed, control pitch angle, and blade geometry details by iteratively solving for the velocity behind the blade based on the momentum 4 of 11 American Institute of Aeronautics and Astronautics taken out of the blade at the turbine disk. Forces at the turbine disk are calculated by integrating lift and drag, calculated by lookup table from the airfoil property curves, over the blade radius. The iteration is enhanced by a relaxation factor, which is iteratively lowered to counteract instability. The airfoil lift and drag coefficients are taken from a series of wind tunnel tests described in7 and8 and extrapolated to post-stall regions by assuming flat plate properties as described in.9 F. Structure Module The Structure module calculates the structural stresses based on the aerodynamic loading and rotational speed of the blade. The airfoil section is replaced with an equivalent rectangle with the same area and chord length as the airfoil, and the entire blade is treated as a cantilevered beam. The aerodynamic loads acting at the quarter chord point at each station are transferred to the center of the rectangular shell with appropriate moments. The bending moments and the centripetal force due to the rotation of the beam cause normal stresses, which are combined with shear stresses to calculate the maximum stress in the blade section. The shell thickness at each station is the decision variable, t. A spar box is also added to the blade geometry to stiffen the blade. The breadth of the spar box is quarter of the chord at that section and the height is the inner height of the rectangular shell. The box thickness is optimized over the range of 5 to 20% of the breadth of the box in order to minimize cross-sectional area without exceeding the maximum allowable stress. If the thickness is tight at 20% at any section, the maximum allowable stress is violated at that section. G. Validation Unfortunately, we did not have a real wind turbine design to compare our codes with. However, we checked the Structure module against analytical results for a rectangular cantilevered beam, and found it to output identical values. We attempted to use higher fidelity models such as Qprop and VABS to check our codes against, but were unable to make working input files in time. IV. Optimization Our formulation has both a nonlinear, non-convex objective function and a nonlinear, non-convex feasible region. Function evaluations are non-trivial and can range from several seconds to ten minutes, depending on the design vector. However, we employ only continuous decision variables, and believe both our objective and constraints are smooth functions. Because of the multiple local optima and expensive function evaluations, our problem lends itself to design space exploration and the use of multi-start gradient methods. While heuristics offer algorithms that can work well to find global optima, they are typically very expensive to run, and in practice we found that we never had enough computational resources to fully converge with heuristic methods including genetic algorithms and simulated annealing. A. Design Space Exploration We performed a design space exploration by generating a space-filling Design of Experiments (DOE) using the orthogonal array method as described in.10 We chose factors and levels as presented in Table 3 and Figure 2. The choice of factors and levels was derived from extensive point-testing of our analysis code and common-sense. For example, the levels for Factor 1, corresponding to blade radius, were chosen to fill only the upper portion of the allowable range for the blade radius as a larger radius increases the total amount of wind energy entering the turbine disk. For other factors that represent distributions, we chose levels to either test different distribution shapes, or relative magnitudes. At least one level in each factor was chosen to correspond with optimal decision variables from a gradient-based optimization we performed early in our project. 1. DOE Results Surprisingly, no DOE test point exceeds the allowable stress. The main effects of each factor are shown in Table 4. To avoid using the penalty method on the constraint, we calculated a main effect for both the maximum stress and the expected power over blade volume, which can be combined to estimate the objective function and maximum stress of any combination of factors. 5 of 11 American Institute of Aeronautics and Astronautics Table 3: DOE Factors and Levels 1 2 3 4 5 6 7 8 Variable Levels R Qmax t T F C β k 3 3 3 3 3 3 3 2 Twist Distribution Factor Levels 20 10 0.4 0.6 0.8 Non−dimensional blade radius, r/R Level 3 16.15 12000 0.004 (0.6, 0.7) 0.25 (2, 1.4, 0.4) (0.7, 0.5, 0.3) 3 14.15 16000 0.01 (0.1, 0.7) 0.5 (0.5, 0.35, 1) (0.4, 0.4, 0.4) 4 12.15 20000 0.02 (0.35, 0.7) 0.75 (1, 0.7, 0.2) (0.8, 0.7, 0.05) - 1 Chord Length Distribution Factor Levels 1 0.8 0.4 1 0.5 0.2 0 0.2 0.4 0.6 0.8 Non−dimensional blade radius, r/R 1 0 0.2 Magnitude Magnitude Magnitude Shape Shape Magnitude Shape Magnitude 50 Level 1 Level 2 Level 3 1.5 0.6 Focus Pitch Control Curve Factor Levels 2 Level 1 Level 2 Level 3 Chord Length, m Foil shape parameter, non−dimensional Twist, degrees Level 1 Level 2 Level 3 30 0 0.2 Level 2 Foil Shape Distribution Factor Levels 50 40 Level 1 Global Blade Pitch, deg Factor 0.4 0.6 0.8 Non−dimensional blade radius, r/R Level 1 Level 2 Level 3 40 30 20 10 1 0 0 5 10 15 Incoming Velocity, m/s (scale varies with k) Figure 2: Graphical representation of factor levels for distribution design variables. Table 4: DOE Factor Main Effects r σmax PE /Vblades Factor 1 2 3 4 5 6 7 8 Level 1 Level 2 Level 3 Level 1 Level 2 Level 3 -1286 -114 3374 4148 -105 -670 251 -250 1438 865 -783 -2919 8 1862 -2319 250 -152 -751 -2591 -1229 96 -1192 2069 - 0.112 -0.962 1.895 0.955 0.322 -2.410 0.087 -0.535 0.983 0.685 -0.778 -0.183 0.297 3.592 -0.855 0.535 -1.095 0.277 -1.117 -0.772 -0.618 -1.182 0.768 - 6 of 11 American Institute of Aeronautics and Astronautics 20 2. Main Effect Extrapolation To find the best combination of main effects, without considering two-factor interaction, we chose the factor levels that corresponded to the highest main effects for the objective function. This combination yielded the levels (2,2,1,1,3,2,3,2), with an expected objective value of PE /Vblades = 20686 W/m3 and a maximum stress r of σmax = 12.22 MPa. Trying this point in our analysis code yielded an objective function value roughly twice as large, PE /Vblades = 44527 W/m3 , and a slightly higher maximum stress, which coincided with the r = 14.00 MPa. maximum allowable σmax B. 1. Gradient-based Algorithms Sequential Quadratic Programming (SQP) For a smooth objective function like ours, gradient-based algorithms offer attractive convergence to local optima. SQP is based on solving a series of sub-problems designed to minimize a quadratic model of the objective subject to linearization of the constraints.11 It mimics Newtons method for constrained problems, and for unconstrained problems, it reduces to Newtons method. At each major iteration an approximation of the Hessian of the Lagrangian function is done using quasi-Newton updates.12 We implemented SQP with MATLAB®'s fmincon function. 2. Scaling SQP convergence is highly sensitive to the condition number of the Hessian. The condition number is defined as the ratio of the highest Eigenvalue to the lowest Eigenvalue of a matrix. A high condition number squashes the level sets of the objective function and can lead to cross-stitching behavior that can stall gradient-based algorithms. In order to alleviate this problem, we halt the optimization process when the function value does not change between iterations within some tolerance, and scale the design variables. To scale the design variables, we estimate the diagonal elements of the Hessian matrix by the standard central difference finite difference scheme Hii = ∂2J J(x1 , . . . , xi + Δx, . . . , xn ) − 2J(x1 , . . . , xi , . . . , xn ) + J(x1 , . . . , xi − Δx, . . . , xn ) ≈ 2 ∂xi Δx2 (3) 1/2 and subsequently apply a scaling factor of Hii to each design variable before passing it to the SQP algorithm. The SQP algorithm is then restarted working with the scaled variables. This process is repeated until the condition number of the Hessian at the converged solution is lower than 10. 3. Multi-start Because our objective function and feasible regions are not convex, we have no guarantee of finding the global optimum with gradient-based algorithms, which only find local optima. Therefore, we sampled the design space by starting the SQP algorithm with multiple start vectors. For simplicity, we used the space-filling orthogonal array points generated for the DOE as starting points for SQP. 4. Convergence Figure 3 shows the improvement in objective function values, using the penalty-method objective, over successive iterations for the SQP algorithm using several starting points. This plot is intended to show the difference in not only convergence rates, but also in solution values for different starting points. In practice, we tried more starting points than could be easily displayed. V. A. Optimality and Post-Optimality Analysis Best Design Vector The best single-objective design vector, found by SQP, is x = (14.13, 20000, 0.004, 4, π/4, π/4, 0.25, 1.91, 0.79, 0.24, π/4, 0.75, 0.36162). The blade shape and performance of this design point is graphically depicted in Figures 4, 5, and 6. 7 of 11 American Institute of Aeronautics and Astronautics 0 −10 Objective, PE/Vblades + ξ −20 −30 −40 −50 −60 −70 −80 0 5 10 15 20 25 Iteration 30 35 40 45 Figure 3: Convergence for selected SQP runs. B. 1. Post-Optimality Analysis Tight Constraints In our best design point, Qmax, T , k, and βcutin were at their upper bounds, and the blade shell thickness, t, was at its lower bound. These tight box constraints make sense because they allow for the biggest generator to be fitted, with the most aggressive pitch angles at slow speeds and close to the hub, and the thinnest blade shell thickness, which minimizes volume. 2. Sensitivity Analysis Table 5 shows the sensitivity of all three components of the objective function with respect to the design variables at the best design, normalized with respect to the best design variable values. As expected, the blade volume is independent of the twist angle and pitch control curve. Also, the expected power depends mainly on the radius, twist angles, and pitch control curve, which heavily influence the aerodynamic forces. Finally, the stresses are mainly dependent on the twist angles, which heavily influence torsional strength. The values in Table 5 could be used to estimate the change in the objective function for points close to the best design point. Table 5: Normalized Jacobian at optimal design vector. J1 = PE J2 = Vblades J3 = σmax J1 = PE J2 = Vblades J3 = σmax 3. ∂J ∂R ∂J ∂Qmax ∂J ∂t ∂J ∂k ∂J ∂Tmid ∂J ∂Ttip ∂J ∂Fmid 1.510 1.0000 0.0823 0.041 0 -0.0161 0 0.8557 0.0434 0.367 0 0.1015 0.858 0 0.1661 0.2966 0 0.1856 -0.007 0.0024 -0.0005 ∂J ∂Croot ∂J ∂Cmid ∂J ∂Ctip ∂J ∂βcutin ∂J ∂βmean -0.122 0.5382 0.0260 0.368 0.5638 0.0024 0.150 0.0486 -0.0010 0.067 0 -0.0706 1.787 0 0.0791 ∂J ∂βcutof f 0.045 0 0.0115 Multi-objective Optimization In practice, it may be more valuable to identify trade-offs between multiple objective functions. In our case, we decided to find the design trade-offs between the expected power and the cost, where the cost was simply a scalar multiple of the blade volume. In this case, we had little luck with gradient-based methods such 8 of 11 American Institute of Aeronautics and Astronautics Foil Parameter Distribution 0 = Root Foil (S814) 1 = Tip Foil (S813) Blade Twist Distribution Foil Parameter 1 40 20 0.4 0.6 0.8 Non−dimensional Blade Radius (x) Chord Length Distribution Chord Length (m) 2 1.5 1 0.5 0 0.2 0.4 0.6 0.8 Non−dimensional Blade Radius (x) 1 Spar Box Thickness (boxt) Structural Spar Thickness wrt. Spar Box Width Distribution Max Box Thickness = 0.062 0.2 0.15 0.1 0.05 0 0.2 0.4 0.6 0.8 Non−dimensional Blade Radius (x) 1 0.5 0 0.2 1 Global Pitch Angle (beta) (deg) 0 0.2 0.4 0.6 0.8 Non−dimensional Blade Radius (x) 1 Global Pitch Angle Control Curve 60 40 20 0 2 4 1.6 6 8 10 12 Incoming Velocity (v0) (m/s) 14 Maximum Stress Distribution Max Stress = 1.40E+007 7 Maximum Stress (N/m2) Blade Twist (deg) 60 x 10 1.4 1.2 1 0.8 0.2 0.4 0.6 0.8 Non−dimensional Blade Radius (x) 1 Figure 4: Best design vector decision variables and outputs. Power (W) Expected Power = 19403 5 Power (W) 10 Turbine Power Output Energy in Air Betz Limit 5 0 2 3 4 5 6 7 8 9 Incoming Wind Velocity (v0) (m/s) 10 11 12 10 11 12 10 11 12 10 11 12 Efficiency (Percentage of Betz Limit Achieved) Blade Radius = 14.13 m, Qmax = 20000 N−m 80 60 40 20 0 Turbine Output Torque (N−m) Turbine Angular Velocity (omega) (rad/s) Efficiency as a Percentage of Betz Limit (%) x 10 2 3 4 5 4 2 6 7 8 9 Incoming Wind Velocity (v0) (m/s) Torque (N−m) x 10 1.5 1 0.5 0 2 3 4 5 6 7 8 9 Incoming Wind Velocity (v0) (m/s) Turbine Angular Velocity 5 4 3 2 1 2 3 4 5 6 7 8 9 Incoming Wind Velocity (v0) (m/s) Figure 5: Best design vector performance diagram. 9 of 11 American Institute of Aeronautics and Astronautics 12 12 10 10 Blade Radius (m) Blade Radius (m) Blade Shape Along Blade Radius Blade Twist Along Blade Radius 8 6 4 8 6 4 2 2 0 0 −0.47751.4325 Chord Length (m) 0 20 40 60 Twist Angle (deg) Figure 6: Best design vector blade sections and twist curve. as weighted-sum and normal boundary intersection methods, and turned to the Multi-Objective Genetic Algorithm (MOGA) heuristic method, as implemented in MATLAB's Optimization Toolbox. The resulting Pareto front is shown in Figure 7, with the best point from SQP shown as a circle on the plot. 800 Cost, $ = Vblades × 1000$/m3 700 600 500 400 300 200 100 0 0 0.5 1 Expected Power 1.5 2 4 x 10 Figure 7: Pareto front for cost and expected power output. In Figure 7, the utopia point is located in the bottom right hand corner, at the maximum power output and minimum cost. Clearly, the SQP result is much better than the MOGA Pareto front. Unfortunately, we were unable to solve our implementation issues with the weighted-sum SQP approach, or to run the MOGA long enough to find the true Pareto front during the term. If the shape of the Pareto front happens to hold when expanded to the actual front, the marginal cost to increase power appears to be linear at low power levels, with a steep increase in cost at higher power levels. VI. Conclusions and Future Work Our work shows promising results from the DOE and SQP runs with a final design that achieves over 60% of the Betz limit for power output for most of the incoming wind distribution. In the future, we would like to extend the analysis to include more decision variables in the input distributions, extend the code to use higher fidelity analysis codes such as VABS for the structures and Qprop for the aerodynamics, and thoroughly validate our codes. Finally, we would like to use a more powerful optimization package, such as DAKOTA to perform the global optimization. 10 of 11 American Institute of Aeronautics and Astronautics References 1 Casas, V.D., F.L. Pena, and R.J.Duro, Automatic Design and Optimization of Wind Turbine Blades, Integrated Group for Engineering Research, University of Corunna (Spain). 2 Kenway, G. and J.R.R.A. Martins, Aerostructural Shape Optimization of Wind Turbine Blades Considering Site-Specific Winds, 12th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Victoria, British Columbia, Canada, 2008. 3 Wang, Q., Personal Interviews, MIT, Spring, 2010. 4 http://wind.nrel.gov/airfoils/Shapes/S814_Shape.html. 5 Freund, R., Interior-Point Theory for Convex Optimization, MIT, Cambridge, 2010. 6 Seguro, J.V. and T.W. Lambert, Modern Estimation of the Parameters of the Weibull Distribution for Wind Energy Analysis, Journal of Wind Engineering and Industrial Aerodynamics, Elsevier, 2000. 7 Ramsay, R.R. and G.M. Gregorek, Effects of Grit Roughness and Pitch Oscillations on the S813 Airfoil, NREL, 1996. 8 Janiszewska, J.M.,R.R. Ramsay, M.J. Hoffman, and G.M. Gregorek, Effects of Grit Roughness and Pitch Oscillations on the S814 Airfoil, NREL, 1996. 9 Tangler, J. and J.D. Kocurek, Wind Turbine Post-Stall Airfoil Performance Characteristics for Blade-Element Momen­ tum Methods, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2005. 10 Kuhfeld, W.F., Orthogonal Arrays, S.A.S., http://support.sas.com/techsup/technote/ts723.html. 11 http://en.wikipedia.org/wiki/Sequential_quadratic_programming. 12 Sequential Quadratic Programming, MATLAB Optimization Toolbox Documentation. 11 of 11 American Institute of Aeronautics and Astronautics