Constrained optimization in python

Question

In Python, using SciPy, I need to find the maximum of a function Q*((1+y)*2-3*Q-0.1)-y**2 Given the restriction (1-z)*(Q*((1+y)*2-3*Q-0.1))-y**2=0. For z I would like to input some values to find the arguments that maximize the function given that value of z.

I have tried many ways to use SciPy optimization functions but I cannot figure out how to go about doing this. I did succeed in doing it using WolframAlpha, but this does not provide me with an answer to questions that follow up on this one.

Attempt:

from scipy.optimize import minimize

def equilibrium(z):
    #Objective function
    min_prof = lambda(Q,y): -1*(Q*((1+y)*2-3*Q-0.1)-y**2)

    #initial guess
    x0 = (0.6,0.9)

    #Restriction function
    cons = ({'type': 'eq', 'fun': lambda (Q,y): (1-z)*(Q*((1+y)*2-3*Q-0.1))-y**2})

    #y between 0 and 1, Q between 0 and 4
    bnds = ((0,4),(0,1))

    res = minimize(min_prof,x0, method='SLSQP', bounds=bnds ,constraints=cons)

    return res.x

from numpy import arange

range_z = arange(0,1,0.001)

print equilibrium(range_z)

Error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-527013574373> in <module>()
     21 range_z = arange(0,1,0.01)
     22 
---> 23 print equilibrium(range_z)

<ipython-input-20-527013574373> in equilibrium(z)
     14     bnds = ((0,4),(0,1))
     15 
---> 16     res = minimize(min_prof,x0, method='SLSQP', bounds=bnds ,constraints=cons)
     17 
     18     return res.x

/Users/Joost/anaconda/lib/python2.7/site-packages/scipy/optimize/_minimize.pyc in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    456     elif meth == 'slsqp':
    457         return _minimize_slsqp(fun, x0, args, jac, bounds,
--> 458                                constraints, callback=callback, **options)
    459     elif meth == 'dogleg':
    460         return _minimize_dogleg(fun, x0, args, jac, hess,

/Users/Joost/anaconda/lib/python2.7/site-packages/scipy/optimize/slsqp.pyc in _minimize_slsqp(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, **unknown_options)
    324             + 2*meq + n1 + ((n+1)*n)//2 + 2*m + 3*n + 3*n1 + 1
    325     len_jw = mineq
--> 326     w = zeros(len_w)
    327     jw = zeros(len_jw)
    328 

ValueError: negative dimensions are not allowed

The above link provides an example of performing constrained optimization. If you attempt to solve your problem and run into an issue, please share your code attempt, and the specific error you run into and we can help from there. — Cory Kramer, Jun 12 '17 at 13:06
@CoryKramer I rephrased my problem, I get the solution now for having a particular value for `z` (which was 0.3), now I would like to study how these optima react to the value of `z`. — Joost Bouten, Jun 12 '17 at 13:53

score 0 · Answer 1 · answered Jun 14 '17 at 15:20

You need to evaluate your function for one z at a time. A minimal modification to make work your code is as follows:

print [equilibrium(z) for z in z_range]

In your current code the function describing the constraints returns a vector instead of a scalar, which leads to the error message.

Instead of optimizing numerically you might note that your problem is solvable analytically:

a = 0.1
Q = (6-3*a+3**.5 *(4-4*a+a**2-4*z+4*a*z-a**2 *z)**.5)/(6*(2+z))
y = Q*(1-z)+(Q*(-1+z)*(-2+a+Q*(2+z)))**.5

You can test this and convince yourself that it gives the same result (up-to numerical precision) then the numerical optimization. (I've tested for z=0.745 - you need to check the 2nd derivative to select the correct maximum generally. But this is doable.)

Constrained optimization in python

1 Answers1