This work addresses many problems associated with designing aerodynamic shapes using computational fluid dynamics (CFD) codes. The investigation focuses in the transonic flow regime where shock waves may have an adverse effect on the convergence of the optimization process. In particular, the interaction of the flow discontinuity with the discrete representation of the design problem may cause the objective function to be non-smooth. Methods for robust optimization of the nonsmooth functions are presented. The dissertation is divided into two parts, The first part investigates a simple model problem involving quasi-one-dimensional flow in a duct. The flow field computation is simple and contains many of the elements present in more complicated fluid flow problems. The optimization involves finding the cross-sectional area distribution of a duct that produces velocities which closely match a targeted velocity distribution containing a shock wave. The objective function which quantifies the difference between the targeted and calculated velocity distribution becomes nonsmooth due to the presence of the shock in the discretized field. Two techniques for derivative-based optimization are offered to resolve the difficulties associated with the non-smoothness of the objective function. The first technique, shock-fitting, involves careful integration of the objective function through the shock wave. The second technique, coordinate straining with shock penalty, uses a coordinate transformation to align the calculated shock with the target and then adds a penalty proportional to the square of the distance between shocks. These techniques are evaluated and tested using several methods to compute the derivatives, including finite-differences, direct and adjoint methods. The above two techniques rely on accurate estimations of the shock position, which may not be available for the general case. In the second part of the dissertation, we present an optimization method to solve the difficult model design problem requiring no information about the shock. The optimization begins with the construction of a response surface that smoothly approximates the objective function. Here the response surface is a least squares polynomial fit to carefully selected design points. By minimizing the response surface we can obtain a first guess for a reasonable design. Optimization may continue in one of two ways. In the first method, we probe a small region of the design space around the minimum and perform another response surface minimization. In the second method we switch to a derivative-based method assuming that in the small region around the minimum the function is smooth. In addition to the one-dimensional duct problem, two other shape design problems involving two-dimensional flow are solved to demonstrate the efficacy and robustness of the response surface method. One involves the inverse design of a bump in a transonic channel flow. The other involves the design of an airfoil for transonic flight. / Ph. D.