Suppose I have a set $S=\{1,2,\dots,500\}$ and some function $f(\sigma)$ from the permutations $\operatorname{Perm}(S) \rightarrow \mathbb{R}$ to be minimized.
The function is complex (simulation based), black box, non-smooth, and since dealing with a discrete space there is no gradient defined. However there is reason to believe that there are strong relationships between the positions of some elements and the value of $f(\sigma)$ for example when $17 \in S$ is "earlier" in the permutation then this contributes to minimizing $f$ and when $305\in S$ is "later" in the permutation then this also helps minimize $f$. There are other relationships such as if $s \in S$ is earlier than $s' \in S$ in the permutation then this is also "good".
In other words we are seeking the optimal ranking of the elements where the objective function is really defined on the ranks, not the permutation indices.
A first attempt I made at optimisation is to define entity weights $\theta\in[0,1]^{500}$ such that the ranks within $\theta$ define the permutation/ranks of the entities, inspired by the "random keys" approach in genetic algorithms. Then I fed $\min f(\operatorname{ranks}(\theta))$ into Nelder-Mead function to minimize and also compared to uniform random sampling within the parameter space. Nelder-Mead typically does not out-perform random sampling (given same limit on function evaluations) and often gets stuck in a poor local minimum. Neither approach leverages the ranked structure in the parameter space which has a huge degree of symmetry (e.g. any perturbation of $\theta$ which does not alter the ranking will not change the objective value).
What is a more natural approach to optimising the rankings of a discrete, finite set of entities with respect to some black box objective function?
P.S. I am looking for practical and effective methods to apply - global optimality/analytic solutions are not necessary. A reasonable function call limit in this context would be 10000-50000.