Question: I am having problems minimizing a simple if slightly idiosyncratic function. I have scipy.optimize.minimize but I can't get consistent results. Here is the full code: from math import log, exp, sqrt from bisect import bisect_left from scipy.optimize import minimize from scipy.optimize import Bounds import numpy as np def new_inflection(x0, x1): return log((exp(x0)+exp(x1) + sqrt(exp(2*x0)+6*exp(x0+x1)+exp(2*x1)))/2) def make_pairs(points): new_points = [] for i in range(len(points)): for j in range(i+1, len(points)): new_point = new_inflection(points[i], points[j]) new_points.append(new_point) return new_points def find_closest_number(numbers, query): index = bisect_left(numbers, query) if index == 0: return numbers[0] if index == len(numbers): return numbers[-1] before = numbers[index - 1] after = numbers[index] if after - query < query - before: return after else: return before def max_distance(target_points): pair_points = make_pairs(target_points) target_points = sorted(target_points) dists = [] return max(abs(point - find_closest_number(target_points, point)) for point in pair_points) num_points = 20 points = np.random.rand(num_points)*10 print("Starting score:", max_distance(points)) bounds = Bounds([0]*num_points, [num_points] * num_points) res = minimize(max_distance, points, bounds = bounds, options={'maxiter': 100}, method="SLSQP") print([round(x,2) for x in res.x]) print(res)  Every time I run it I get quite different results. This is despite the output saying Optimization terminated successfully. An example output:  message: Optimization terminated successfully success: True status: 0 fun: 0.4277378933292031 x: [ 5.710e+00 1.963e+00 ... 1.479e+00 6.775e+00] nit: 15 jac: [ 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00] nfev: 364 njev: 15  Sometimes I get a result as low as 0.40 and other times as high as 0.51. Is there any way to optimize this function properly in Python?

Answer: The problem here is that you are exploring immense non convex search space by trying to optimize 20 variables at the same time, so usual optimization methods such as gradient descent related etc. will potentially be trapped in a local minima. In your case, as you are starting from a random coordinates each time, you end in a different local minima each time. If the problem has no analytical solution, there is no optimizer (at least that I know) that can solve this kind of problem and guarantee you that you are in the global minimum, but that doesn't mean that we can get a very good approximations. Approach One Try different types of algorithms more suited for non convex high dimension space optimization, within scipy you have a global optimizer function: from scipy import optimize optimize.shgo(eggholder, bounds)  In my case that was very slow, but maybe it can help you. Update: The global optimizer basinhopping seems to give good results faster: from scipy.optimize import basinhopping res = basinhopping(max_distance, points, minimizer_kwargs={"method": "SLSQP", "bounds": bounds}, niter=100) print([round(x,2) for x in res.x]) print(res)  With this optimizer you can also reach a fitness of 3.9 consistently. Approach Two I would give a try to genetic algorithms using pygad, here is my trial and it reaches 3.9 fitness and consistently bellow 4.1 (though I change the sign for the optimizer). import pygad def fitness_func(solution, solution_idx): return -max_distance(solution) fitness_function = fitness_func num_generations = 500 num_parents_mating = 5 sol_per_pop = 100 num_genes = len(points) init_range_low = 0 init_range_high = 20 parent_selection_type = "sss" keep_parents = 5 crossover_type = "single_point" mutation_type = "random" mutation_percent_genes = 10 ga_instance = pygad.GA(num_generations=num_generations, num_parents_mating=num_parents_mating, fitness_func=fitness_function, sol_per_pop=sol_per_pop, num_genes=num_genes, gene_space = {'low': 0, 'high': 20}, init_range_low=init_range_low, init_range_high=init_range_high, parent_selection_type=parent_selection_type, keep_parents=keep_parents, crossover_type=crossover_type, mutation_type=mutation_type, mutation_percent_genes=mutation_percent_genes) ga_instance.run()  To show the solution: solution, solution_fitness, solution_idx = ga_instance.best_solution() print("Parameters of the best solution : {solution}".format(solution=solution)) print("Fitness value of the best solution = {solution_fitness}".format(solution_fitness=solution_fitness)) # plots the optimization process ga_instance.plot_fitness(title="PyGAD with Adaptive Mutation")  Approach Three Probably the worst idea but also maybe is everything you need, Just iterate many times the algorithm you already have, with different initial conditions saving the best solution. Update As the OP commented, somewhat similar to this is what the algorithm basinhopping is doing (with a bit more complexity) but it archives very good results, see approach one. Conclusion There is very little you can do for having consistent results in the optimization algorithm in a problem like this one (if not fixing the seed, which I strongly discourage), but at least you can choose an algorithm that is more suited for the task maximizing the search space, so you get closer and closer to the global minimum. On the other hand you should notice that there may be multiple/infinite solutions with the same fitness, but that can look very different, specially if there are hidden symmetries in the problem, this problem for instance seems invariant under permutations, so it would make sense to sort the solution list. Also in this case little changes in one element of the 20 numbers doesn't necessarily change the fitness.