如何使用python在短时间内获得具有两个变量积分表达式的成本函数的最小值？
How to get the minimum value of a cost function, having two variable integration expression, in short time using python?
我想找到成本函数 T 的最小值。成本函数 T 在两个变量（Q 和 r）中有一个表达式。我还需要找到成本函数 T 达到全局最小值的 Q 和 r 的值。 （如果有多个全局最小值 - 那么全部）Q 和 r 的界限是：0 < Q < 15000 ； 0 < r < 5000
这是等式
我正在使用 Sympy library to generate the equations. and using the minimize function of scipy.optimize.minimize 来查找最小值。这些函数的代码是：
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)
但是报错如下。
---------------------------------------------------------------------------
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.
另外，用其他方法，需要很长时间运行，30多分钟左右，最后会报错。 如何在很短的时间内找到全局最小值以及 Q 和 r 的值。最好1-5分钟左右。
Posting on behalf of my Friend
未来的注意事项：在您的函数 f 中，如果您将 r 和 Q 设置为某些内容，它不会更改您之后使用的 SymPy 表达式，因为它们之前已经为符号变量定义了。
您的工作似乎大量使用数字，事实上，由于您的答案不需要符号，您最好进行非符号积分。 SymPy 是纯粹的 Python ，它可能会很慢，特别是对于集成来说，而 SciPy 被设计得很快。这就是为什么我将所有内容都转换为 SciPy 的原因：
编辑：我知道我对 r=0 收敛的第一个答案是可疑的。在@VishalAnand 对从 -inf 开始的 cdf 积分进行校正后，我再次尝试 运行 该程序。 T 的单次迭代花费了约 15 秒，但没有找到解决方案；可能是由于现在存在一个非常复杂的表面。
cdf 也产生了错误的值；例如，quad(pdf, -np.inf, 50000)[0] 产生了一个非常接近 0 的数字，而它应该接近 1。这破坏了最小化，所以我尝试了类似 quad(pdf, -1000000, 50000)[0] 的东西，结果产生了与 sympy.N(sympy.erf((x-mean)/(sqrt(2)*std)))/2 + 1/2 结果证明计算速度更快。
问题是 SciPy 最小化函数无法收敛，而是产生了 ABNORMAL_TERMINATION_IN_LNSRCH。所以我给了它一个特定的使用方法：Nelder-Mead。这收敛了。但最终值非常令人担忧，因为它们在 inf 和 -1.793193606659277e+19 之间跳跃。 Python 不知道溢出错误（至少据我所知）所以我能想到的唯一可能的解释是函数 C 有一个根导致 T 在 r 和 Q 的某些值处有一个渐近线。
这远远超出了我的范围，所以我将在这里留下我的最新成果：
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)
在大约 15 秒内产生以下输出：
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])
希望更有资格的人能解释一下。对我造成的任何不便或错误结论深表歉意。
我想找到成本函数 T 的最小值。成本函数 T 在两个变量（Q 和 r）中有一个表达式。我还需要找到成本函数 T 达到全局最小值的 Q 和 r 的值。 （如果有多个全局最小值 - 那么全部）Q 和 r 的界限是：0 < Q < 15000 ； 0 < r < 5000 这是等式
我正在使用 Sympy library to generate the equations. and using the minimize function of scipy.optimize.minimize 来查找最小值。这些函数的代码是：
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)
但是报错如下。
---------------------------------------------------------------------------
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.
另外，用其他方法，需要很长时间运行，30多分钟左右，最后会报错。 如何在很短的时间内找到全局最小值以及 Q 和 r 的值。最好1-5分钟左右。
Posting on behalf of my Friend
未来的注意事项：在您的函数 f 中，如果您将 r 和 Q 设置为某些内容，它不会更改您之后使用的 SymPy 表达式，因为它们之前已经为符号变量定义了。
您的工作似乎大量使用数字，事实上，由于您的答案不需要符号，您最好进行非符号积分。 SymPy 是纯粹的 Python ，它可能会很慢，特别是对于集成来说，而 SciPy 被设计得很快。这就是为什么我将所有内容都转换为 SciPy 的原因：
编辑：我知道我对 r=0 收敛的第一个答案是可疑的。在@VishalAnand 对从 -inf 开始的 cdf 积分进行校正后，我再次尝试 运行 该程序。 T 的单次迭代花费了约 15 秒，但没有找到解决方案；可能是由于现在存在一个非常复杂的表面。
cdf 也产生了错误的值；例如，quad(pdf, -np.inf, 50000)[0] 产生了一个非常接近 0 的数字，而它应该接近 1。这破坏了最小化，所以我尝试了类似 quad(pdf, -1000000, 50000)[0] 的东西，结果产生了与 sympy.N(sympy.erf((x-mean)/(sqrt(2)*std)))/2 + 1/2 结果证明计算速度更快。
问题是 SciPy 最小化函数无法收敛，而是产生了 ABNORMAL_TERMINATION_IN_LNSRCH。所以我给了它一个特定的使用方法：Nelder-Mead。这收敛了。但最终值非常令人担忧，因为它们在 inf 和 -1.793193606659277e+19 之间跳跃。 Python 不知道溢出错误（至少据我所知）所以我能想到的唯一可能的解释是函数 C 有一个根导致 T 在 r 和 Q 的某些值处有一个渐近线。
这远远超出了我的范围，所以我将在这里留下我的最新成果：
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)
在大约 15 秒内产生以下输出：
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])
希望更有资格的人能解释一下。对我造成的任何不便或错误结论深表歉意。