scipy.optimize minimize inconsistent results

I am getting some very weird results when running the minimize function from scipy optimize.

Here is the code

from scipy.optimize import minimize

def objective(x):
    return - (0.05 * x[0] ** 0.64 + 0.4 * x[1] ** 0.36)

def constraint(x):
    return x[0] + x[1] - 5000

cons = [{'type':'eq', 'fun': constraint}]

when running

minimize(objective, [2500.0, 2500.0], method='SLSQP',  constraints=cons)

i get allocation 2500 for each element of x. with fun: -14.164036415985395

With a quick check, this allocation [3800, 1200] gives -14.9

It is highly sensitive to the initial conditions also .

Any thoughts as to what i am doing wrong

the two functions plotted

UPDATE It actually returns the initial conditions.

If i try this though

def objective(x):
    return - (x[0] ** 0.64 + x[1] ** 0.36)

def constraint(x):
    return x[0] + x[1] - 5000.0

cons = [{'type':'eq', 'fun': constraint}]

minimize(objective, [2500, 2500], method='SLSQP', constraints=cons)

returned results seem to be just fine (i have changed the objective function)

This has a few upvotes, so I thought it deserves an answer.

Any thoughts as to what i am doing wrong

Nothing, really. This is just an issue with problem scaling and the default tolerances. An important default termination tolerance of method='slsqp' is 1e-6 (_minimize_slsqp). By setting that to 1e-7:

minimize(objective, [2500.0, 2500.0], method='SLSQP', constraints=cons, tol=1e-7)

You get the solution you're expecting:

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: -14.913839470893286
       x: [ 3.901e+03  1.099e+03]
     nit: 15
     jac: [-1.630e-03 -1.630e-03]
    nfev: 45
    njev: 15

Notice that at termination, x is on the order of 1e3 and jac is on the order of 1e-3. While the details are more complicated (it probably also has something to do with relationship between the scales of the constraint Jacobian and the gradient), I don't think the 6 orders of magnitude scale difference is a coincidence. If we rescale your problem by making a variable substitution x = x * scale in the objective and constraints:

import numpy as np
from scipy.optimize import minimize

scale = 1.5  # it doesn't take much

def objective(x):
    x = x*scale
    return - (0.05 * x[0] ** 0.64 + 0.4 * x[1] ** 0.36)

def constraint(x):
    x = x*scale
    return x[0] + x[1] - 5000

cons = [{'type':'eq', 'fun': constraint}]

x0 = np.asarray([2500.0, 2500.0])
minimize(objective, x0/scale, method='SLSQP', constraints=cons)

Then we don't change the initial values of the objective or constraint, yet the solver continues on and eventually converges to the expected answer, even without changing the termination tolerance.

You'll find that if you change scale to 10, you bring the final x and jac two orders of magnitude closer together, and you can increase tol to 1e-5 and still converge to the desired solution, but with tol=1e-4 it doesn't. More generally, if you set tol = 1e-7 * scale**2, it converges to the desired solution, but with tol = 1e-6 * scale**2, it doesn't (for scale within a few orders of magnitude of 1).

So it looks like your problem in the original "units" is just poorly scaled in the sense that the gradient is very small with respect to the decision variables (or again, the constraint Jacobian might matter, too). Bringing these scales a little closer together gives the solver the nudge it needs to find a better solution.