We address bilevel programming problems when the derivatives of both the upper- and the lower-level objective functions are unavailable. The core algorithms used for both levels are trust-region interpolation-based methods, using minimum Frobenius norm quadratic models when the number of points is smaller than the number of basis components. We take advantage of the problem structure to derive conditions (related to the global convergence theory of the underlying trust-region methods, as far as possible) under which the lower level can be solved inexactly and sample points can be reused for model building. In addition, we indicate numerically how effective these expedients can be. A number of other issues are also discussed, from the extension to linearly constrained problems to the use of surrogate models for the lower-level response. One important application of our work appears in the robust optimization of simulation-based functions, which may arise due to implementation variables or uncertain parameters. The robust counterpart of an optimization problem without derivatives falls into the category of the bilevel problems under consideration here. We provide numerical illustrations of the application of our algorithmic framework to such robust optimization examples. © 2012 Copyright Taylor and Francis Group, LLC.