How to get the minimum value of a cost function, having two variable integration expression, in short time using python?
I want to find the minimum value of the cost function T. The cost function T has an expression in two variables (Q and r). I also need to find values of Q and r at which the cost function T reaches the global minimum. (if multiple global minimum values - then all) The bounds of Q and r are : 0 < Q < 15000 ; 0 < r < 5000 Here are the equations
I am using Sympy library to generate the equations. and using the minimize function of scipy.optimize.minimize to find the minimum value. The code for the functions are:
from sympy import *
from scipy.optimize import root_scalar
mean, std = 291, 253
l = 7 #
m = 30
#Q = mean*(lead_time + shelf_life)
p = 5
w = 2
K = 100
c = 5
h = 0.001 #per unit per day
x = symbols("x")
t = symbols("t")
r = symbols("r")
Q = symbols("Q")
#defining Cumulative distribution function
def cdf():
cdf_eqn = (1/(std*sqrt(2*pi)))*exp(-(((t-mean)**2)/(2*std**2)))
cdf = Integral(cdf_eqn, (t,-oo,x)).doit()
return cdf
#defining Probability density function
def pdf():
pdf = (1/(std*sqrt(2*pi)))*exp(-((( (x - mean)**2)/(2*std**2)))).doit()
return pdf
pdf = pdf()
cdf = cdf()
#getting the equation in place
G = K + c*Q + w*(Integral(cdf , (x, 0, Q)) + Integral(cdf.subs(x, (r + Q - x))*cdf , (x, 0, r)))\
+ p*(mean*l - r + Integral(cdf , (x, 0, r)))
CL = (Q - r + mean*l - Integral(cdf , (x, 0, Q)) - Integral(cdf.subs(x, (r + Q - x))*cdf , (x, 0, r)) + Integral(cdf , (x, 0, r)))/mean
I = h*(Q + r - mean*l - Integral(cdf , (x, 0, Q)) - Integral(cdf.subs(x, (r + Q - x))*cdf , (x, 0, r)) + Integral(cdf , (x, 0, r)))/2
#TC.free_symbols
#optimising the cost function
from scipy import optimize
def f(params):
r, Q = params
TC = G/CL + I
return TC
initial_guess = [2500., 10000.]
result = optimize.minimize(f, initial_guess, tol=1e-6 )
if result.success:
fitted_params = result.x
print(fitted_params)
else:
raise ValueError(result.message)
But it throws an error as below.
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
/usr/local/lib/python3.6/dist-packages/scipy/optimize/optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
699 try:
--> 700 df = df.item()
701 except (ValueError, AttributeError):
AttributeError: 'Zero' object has no attribute 'item'
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
5 frames
<ipython-input-6-e9bb4190fef5> in <module>()
39 return TC
40 initial_guess = [2500., 10000.]
---> 41 result = optimize.minimize(f, initial_guess, tol=1e-6 )
42 if result.success:
43 fitted_params = result.x
/usr/local/lib/python3.6/dist-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
602 return _minimize_cg(fun, x0, args, jac, callback, **options)
603 elif meth == 'bfgs':
--> 604 return _minimize_bfgs(fun, x0, args, jac, callback, **options)
605 elif meth == 'newton-cg':
606 return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
/usr/local/lib/python3.6/dist-packages/scipy/optimize/optimize.py in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
1007 else:
1008 grad_calls, myfprime = wrap_function(fprime, args)
-> 1009 gfk = myfprime(x0)
1010 k = 0
1011 N = len(x0)
/usr/local/lib/python3.6/dist-packages/scipy/optimize/optimize.py in function_wrapper(*wrapper_args)
325 def function_wrapper(*wrapper_args):
326 ncalls[0] += 1
--> 327 return function(*(wrapper_args + args))
328
329 return ncalls, function_wrapper
/usr/local/lib/python3.6/dist-packages/scipy/optimize/optimize.py in approx_fprime(xk, f, epsilon, *args)
763
764 """
--> 765 return _approx_fprime_helper(xk, f, epsilon, args=args)
766
767
/usr/local/lib/python3.6/dist-packages/scipy/optimize/optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
700 df = df.item()
701 except (ValueError, AttributeError):
--> 702 raise ValueError("The user-provided "
703 "objective function must "
704 "return a scalar value.")
ValueError: The user-provided objective function must return a scalar value.
Additionally, with other methods, it takes a long time to run, more than 30 minutes or so and ends up throwing an error. How can I find the global minima and also the values of Q and r in a very short time. Preferably 1-5 minutes or so.
Posting on behalf of my Friend
Just a note for the future: in your function f, if you set r and Q to something, it does not change the SymPy expressions that you use afterwards since they were already previously defined for symbolic variables.
Your work seems heavily numerical and in fact, since your answers don't need symbols, you're probably better doing non-symbolic integration. SymPy is pure Python which can be slow especially for integration while SciPy is designed to be fast. That's why I converted everything to SciPy things:
Edit: I knew that my first answer with the r=0 convergence was fishy. After @VishalAnand's correction for the integral of the cdf starting from -inf, I tried running the program again. It took ~15 seconds for a single iteration of T and it was not finding a solution; probably due to a now very complex surface that now exists.
The cdf was also producing wrong values; for example, quad(pdf, -np.inf, 50000)[0] produced a number very close to 0 when instead it should have been close to 1. This ruined the minimization and so I tried something like quad(pdf, -1000000, 50000)[0] which turned out to produce similar values to sympy.N(sympy.erf((x-mean)/(sqrt(2)*std)))/2 + 1/2 which turned out to be faster to compute.
The problem is that the SciPy minimize function could not converge but rather produced ABNORMAL_TERMINATION_IN_LNSRCH. So I gave it a specific method to use: the Nelder-Mead. This converged. But the final values were very concerning since they were jumping between inf and -1.793193606659277e+19. Python is not known for overflow errors (at least to my knowledge) so the only possible explanation I can think of is that the function C has a root that causes T to have an asymptote at some values of r and Q.
This is far beyond me so I'll just leave my updated efforts here:
from numpy import sqrt, pi, exp, inf
from sympy import erf, N
from scipy import optimize
from scipy.integrate import quad
mean, std = 291, 253
l = 7
m = 30
# Q = mean*(lead_time + shelf_life)
p = 5
w = 2
K = 100
c = 5
h = 0.001 # per unit per day
# defining Probability density function
def pdf(x):
return (1 / (std * sqrt(2 * pi))) * exp(-(((x - mean) ** 2) / (2 * std ** 2)))
# defining Cumulative distribution function
def cdf(x):
# cdf1 = quad(pdf, -1000000, x)[0] # slow
# cdf2 = quad(pdf, -inf, x)[0] # slow and produces wrong values at hugh positive x
cdf3 = N(erf((x-mean)/(sqrt(2)*std)))/2 + 1/2
return cdf3
# getting the equation in place
def G(r, Q):
return K + c * Q \
+ w * (quad(cdf, 0, Q)[0] + quad(lambda x: cdf(r + Q - x) * cdf(x), 0, r)[0]) \
+ p * (mean * l - r + quad(cdf, 0, r)[0])
def CL(r, Q):
return (Q - r + mean * l - quad(cdf, 0, Q)[0]
- quad(lambda x: cdf(r + Q - x) * cdf(x), 0, r)[0]
+ quad(cdf, 0, r)[0]) / mean
def I(r, Q):
return h * (Q + r - mean * l - quad(cdf, 0, Q)[0]
- quad(lambda x: cdf(r + Q - x) * cdf(x), 0, r)[0]
+ quad(cdf, 0, r)[0]) / 2
def f(params):
r, Q = params
TC = G(r, Q)/CL(r, Q) + I(r, Q)
return TC
initial_guess = [2343.70601496, 239.89137499]
result = optimize.minimize(f, initial_guess, bounds=[(0, 5000), (0, 15000)], method="Nelder-Mead")
# result = f(initial_guess) # single check
print(result)
Resulting in the following output in ~15 seconds:
final_simplex: (array([[2343.70594323, 257.01581672],
[2343.70594323, 257.01581672],
[2343.70594323, 257.01581672]]), array([-1.79319361e+19, -1.79319361e+19, -1.79319361e+19]))
fun: -1.793193606659277e+19
message: 'Optimization terminated successfully.'
nfev: 360
nit: 155
status: 0
success: True
x: array([2343.70594323, 257.01581672])
Hopefully someone more qualified can explain this. Apologies for any inconvenience or false conclusions from myself.
- matplotlib 3D scatter plot alpha varies when viewing different angles
- How to write very long string that conforms with PEP8 and prevent E501
- Getting Home Directory with pathlib
- how to avoid bot detection on websites using selenium python
- Python mock to create a fake object return a dictionary when any of its attributes are used
- Polars vs. Pandas: size and speed difference
- How to mock.patch a class imported in another module
- Python - error cannot determine truth value of Relational (Newton-Raphson)
- ProcessPoolExecutor logging fails to log inside function on Windows but not on Unix / Mac
- SQLAlchemy ORM Insert or Update when importing from JSON
- django managers vs proxy models
- Pytroch clamp for complex values
- For every identifier select only rows with largest order column
- truth value for Expr is ambiguous in with_columns ternary expansion on dates
- Remove equal characters from two python strings
- Python pyad module can't set UPN
- Macro VS Micro VS Weighted VS Samples F1 Score
- Printing a Tree data structure in Python
- How to fix/reset decreasing timestamps while preserving gaps in time-series data for CNN training?
- Test that module is NOT imported
- Pyserial module isn't installed on PATH
- Print a multiplication table in Python
- Python: ModuleNotFoundError: No module named 'xyz'
- Receiving Import Error: No Module named ***, but has __init__.py
- PyQt5 QProgressBar border radius issue
- URL-encoding and -decoding a string in Python
- Fastest way to find the smallest possible sum of the absolute differences of pairs within a single array?
- Flask: Update Code Reference for: current_app._get_current_object()
- Export Charts from Excel as images using Python
- Align yaxis label spanning two axes with yaxis labels of one axes in subplots