Why a metaheuristic? I think those methods are used way too often. Do you have a good motivation go only consider metaheuristics?

There are several things to think about before selecting an algorithm for such a problem:

1. How fast can you calculate one (1) objective value? I suppose you use a simulation tool for that purpose, so it then depends on how accurate and fast that simulation is.

2. So: in your problem, how accurate is that objective value output? That is, can we say that the simulation provides an accurate representation of the ACTUAL objective value?

If your simulation is fairly accurate then you may utilize a radial basis function (RBF) representation of the (unknown) function, and optimize that explicit function through some good global optimization solver. (These days they are fairly good.) If your "real" function is well-behaved (that is, Lipschitz continuous and lower bounded) then you have a fair chance of being able to solve the problem eventually. If it is not, then I fear you have entered the ugly world of the so-called "No Free Lunch Theorem" - which basically says that you no clue at all of what is going on, and there is no guarantee that any method will give you anything of interest. If the latter is the case, I for one would not even bother, but if this is your case and you NEED to solve it, then, and only then, would I throw every meta-heuristic at it and simply type out the best that you get, without ANY claims that you have an optimal solution of course.

The world is not fair you know - some of the most interesting problems are also very, very hard ... :-)

For my point of view, for black-box function optimization the very good method is Cross-Entropy (see, e.g., Article: "The Cross-Entropy Method for Continuous Multi-Extremal Optimization",  https://www.researchgate.net/profile/Sergey_Porotsky/publications )

During tuning performing (control parameter selection) you can change both sample size, elite part size, smmothing parameters, etc. and to investigate speed of convergence and accuracy. Certainly, it is (i.e. tuning) rather art, not scince, but it is really worked!

I agree with Michael Patriksson that if your black-box function is not well-behaved then it is very difficult to get clues about its topology, structure or general features. Trying a lot of heuristics could be the best way to go. But perhaps you could benefit from applying (meta)-heuristics in organized manner. For instance:

1- You can start by determining whether your function is multi-modal. Use some simple heuristics which are known for being good on uni-modal problems (e.g. Hill Climbing, Nelder-Mead, Evolutionary Strategies). Restart them with different initial solutions, by measuring whether they converge to different/distant local optima, how much they improve the initial solution and how fast they converge you may get some idea on whether your problem is uni-modal, multi-modal or highly multi-modal.

2- Use algorithms which exploit separability (e.g. Coevolutionary Computation) and compare their results to other heuristics. This should allow you to get some clues of whether your problem is separable, partially separable or non-separable at all.

3- Use algorithms which detect correlations between variables (eg. Estimation Distribution Algorithms or CMA-ES). You may use different algorithms with the ability to detect different kinds of correlations: linear, non-linear, non-monotonic, etc.

4- Use algorithms which are good on highly multi-modal problems with good global structure (e.g Minimum Population Search) to have some clue about the multi-modality and conditioning of your problem.

5-… and so on.

From a practical point of view black-box optimization is most interesting, as very many practical(!) systems can only be described procedurally as very complex response functions. A bit in contrast to Chin-Chang Yang's comment (which is of course true from a purely systematic perspective) most practical systems adhere to comparably smooth general functional behaviour (according to Einstein: "Nature doesn't jump"), so in the vast majority of cases I do not expect sudden isolated and sharp optimality peaks. I rather see the main problem of black-box systems in a non-mathematical definition of partial boundaries and internal mode switches to different kinds of responses that lead to non-differentiable jumps to different target value domains that cannot really be anticipated. A practical example: trail breakouts of different cars on different street surfaces in curves, depending on velocity and acceleration. Beyond general functional friction considerations (resulting in rather smooth behaviours in the differing friction domains) only the black-box response in a given situation (tire, surface, wetness, car agility ...) will be able to describe the situation. So the "true" optima (e.g. maximum curve velocity) will quite likely to lie on the edges of restriction regimes that cannot be seriously anticipated without numerous "function calls" to the black box.

The classical engineering approach for characterization of 'black boxes' is impulse response resp. step response.

This primarily applies to "boxes" with 1 input and 1 output, but boxes with more than 1 input can be characterized by scanning the responses of 1 input after the other.

(Boxes with more than 1 output can be considered as more than 1 box.)

Overall analysis and optimization will become anything between 'complicated' and impossible if your black box contains storage elements as they will interact 'unpredictably'  (in terms of the 'black' box) with the input signals.

Thus black box optimization must be restricted to relatively simple 'boxes'. Otherwise you will have to 'open the box'. At least partially...

Memetic algorithms and iterated local search have some mechanisms you described.  If we refer to a mathematical function, we would use a systematic method to find the characteristics and behaviours.  

Thinking aloud, a penalising the solutions for a diverging from an optimum could be one way to tackle the problem. Perhaps the fitness function can help in this aspect.

In response to U. Dreher's comment: I think the initial asker had a "computer generated" response function in mind. You are right w/r/to storaging or delayed response effects. I blanked this out in my answer and assumed implicitly a kind of equilibrium conditions (in the sense of: neither storage nor hysteresis effects).

An approach to storaging effects etc. may be one that I took in my "optimization" of heat transfer understanding in structured opaque wall elements (see my featured publications, available only in German, though): Set up a (kind of, potentially simplified) physical simulation model of the process in question, containing/modelling respective effects, and optimize the time resolved theoretical response function w/r/to its deviation from an observation time series by adapting the modelling parameters that influence the storaging effects as well.

This is not the basic "black box" approach, but rather a "simulated black box by educated guesses" optimization, though.  :-)

Using input output data only........

The only way to analyze such functions is firstly characterized them. This is with a data divern approach to extract an alternative representation. With this respresentation, then, it is posible to analize some ofe its properties.

All the answer are great