Scipy.optimize.minimize function to determine multiple variables

Question

I would like to find out optimal parameters like bo,ho,t1 and t2 such that the error function (f) is minimum (i.e., zero or close to zero). Full code is shown below,

import numpy
import scipy.optimize as optimize

#OPTIMISATION USING SCIPY
def Obj_func(x):
    bo,ho,t1,t2=x
    f=-321226.4817 + (10400000*bo*ho**3 - 10400000*(bo - 2*t2)*(ho - 2*t1)**3)/(125 + (10400000*bo*ho**3 - 10400000*(bo - 2*t2)*(ho - 2*t1)**3)/(1563920*t1*(bo - 2*t2)))
    return f

initial_guess=[2,2,0.2,0.2]
cons = ({'type': 'eq', 'fun': lambda bo,ho,t1,t2: (bo*ho - (bo - 2*t2)*(ho - 2*t1)-7.55)})
bnds = ((0, None), (0, None),(0, None), (0, None)) #all four variables are positive and greater than zero
#Always t1 and t2 should always be lesser than bo and ho
#res=optimize.minimize(Obj_func, method='SLSQP',initial_guess, bounds=bnds,constraints=cons)
res=optimize.minimize(Obj_func,initial_guess, bounds=bnds,constraints=cons)
print ("Result",res)

I am facing some issue here, Q1. lambda function created for cons is giving me an error which states`"TypeError: () missing 3 required positional arguments: 'ho', 't1', and 't2'". Why this is happening? even though i followed the syntax correctly, some thing i am missing. Q2. I need to minimize the objective function i.e., here its minimizing the error function (f) to zero or close to zero. I have mentioned optimize.minimize() without mentioning the type of optimization method. Is it ok? Will it give me the correct answer i am looking for?

Answer 1

Q1 (syntax): Your constraint must take a single vector argument:

import numpy
import scipy.optimize as optimize

#OPTIMISATION USING SCIPY
def Obj_func(x):
    bo,ho,t1,t2=x
    f=-321226.4817 + (10400000*bo*ho**3 - 10400000*(bo - 2*t2)*(ho - 2*t1)**3)/(125 + (10400000*bo*ho**3 - 10400000*(bo - 2*t2)*(ho - 2*t1)**3)/(1563920*t1*(bo - 2*t2)))
    return f

initial_guess=[2,2,0.2,0.2]
cons = ({'type': 'eq', 'fun': lambda bhtt: (bhtt[0]*bhtt[1] - (bhtt[0] - 2*bhtt[3])*(bhtt[1] - 2*bhtt[2])-7.55)})
bnds = ((0, None), (0, None),(0, None), (0, None)) #all four variables are positive and greater than zero
#Always t1 and t2 should always be lesser than bo and ho
#res=optimize.minimize(Obj_func, method='SLSQP',initial_guess, bounds=bnds,constraints=cons)
res=optimize.minimize(Obj_func,initial_guess, bounds=bnds,constraints=cons)
print ("Result",res)

Output:

Result      fun: -15467.04696553346
     jac: array([  165915.125     ,   147480.57421875,  1243506.8828125 ,
        -130354.05859375])
 message: 'Positive directional derivative for linesearch'
    nfev: 6
     nit: 5
    njev: 1
  status: 8
 success: False
       x: array([ 2. ,  2. ,  0.2,  0.2])

Q2 (maths): The solver doesn't seem to like the problem as it is posed. I managed to get it spit out a solution by eliminating the equality constraint and one variable. The default solver still didn't work but COBYLA did:

import numpy
import scipy.optimize as optimize
from operator import itemgetter

#OPTIMISATION USING SCIPY
def Obj_func_orig(x):
    bo,ho,t1,t2=x
    f=-321226.4817 + (10400000*bo*ho**3 - 10400000*(bo - 2*t2)*(ho - 2*t1)**3)/(125 + (10400000*bo*ho**3 - 10400000*(bo - 2*t2)*(ho - 2*t1)**3)/(1563920*t1*(bo - 2*t2)))
    return f

def Obj_func(x):
    bh, h, t = x
    btht = bh - 7.55
    D = bh*h*h - btht*(h-t)*(h-t)
    f = -321226.4817 + D / (125/10400000 + D*(h-t)/(781960*t*btht))
    return f

def cons(x):
    bh, h, t = x
    btht = bh - 7.55
    return (btht * h - bh * (h-t)) * (t-h)

def ge(i):
    return {'type': 'ineq', 'fun': itemgetter(i)}

initial_guess=[1,1,1]
initial_guess=[64,8,0.5]
cons = ({'type': 'ineq', 'fun': cons}, {'type': 'ineq', 'fun': Obj_func}, ge(0), ge(1), ge(2))
res=optimize.minimize(lambda x: Obj_func(x)**2,initial_guess,constraints=cons,method='COBYLA')
print ("Result",res)
bh, h, t = res.x
bo, ho, t1, t2 = bh / h, h, t/2, (bh / h - (bh - 7.55) / (h-t))/2
print('bo, ho, t1, t2, f', bo, ho, t1, t2, Obj_func_orig([bo,ho,t1,t2]))

Answer with OP's initial guess:

# Result      fun: 285138.38958324661
#    maxcv: 0.0
#  message: 'Optimization terminated successfully.'
#     nfev: 47
#   status: 1
#  success: True
#        x: array([  6.40732657e+01,   8.29209901e+00,   6.03312595e-02])
# bo, ho, t1, t2, f 7.72702612177 8.29209901408 0.0301656297535 0.430273240999 533.983510591

Please be advised that these are clearly local minima that depend on the initial guess. For example if I use 1,1,1 for the transformed variables:

# Result      fun: 5829.8452437661899
#    maxcv: 0.0
#  message: 'Optimization terminated successfully.'
#     nfev: 48
#   status: 1
#  success: True
#        x: array([ 4.39911932,  0.93395108,  0.89739524])
# bo, ho, t1, t2, f 4.71022456896 0.933951079923 0.44869761948 45.4519279253 76.3534232616