如何绘制一个多变量积分函数的三维图形，并找到它的全局最小值？

Question

我有一个成本函数 f(r, Q), 在下面的代码中得到。成本函数 f(r, Q) 是两个变量的函数 r和Q. 我想在下面给出的范围内，绘制所有r和Q值的成本函数值，并找到全局最小值，即 f(r, Q).

r和Q的范围分别为 。

0 < r < 5000 5000 < Q < 15000。

该图应在r、Q和f(r,Q)轴上。

成本函数的代码。

from numpy import sqrt, pi, exp
from scipy import optimize
from scipy.integrate import quad
import numpy as np

mean, std = 295, 250
l = 7
m = 30
p = 15
w = 7
K = 100
c = 5
h = 0.001  # per unit per  day


# defining Cumulative distribution function
def cdf(x):
    cdf_eqn = lambda t: (1 / (std * sqrt(2 * pi))) * exp(-(((t - mean) ** 2) / (2 * std ** 2)))
    cdf = quad(cdf_eqn, -np.inf, x)[0]
    return cdf


# defining Probability density function
def pdf(x):
    return (1 / (std * sqrt(2 * pi))) * exp(-(((x - mean) ** 2) / (2 * std ** 2)))


# 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

如何绘制这个函数 f(r,Q) 在三维图中，还可以得到 全局最小值 和价值的 r和Q 在那个特殊的点上。

此外，我已经尝试使用 scipy.optimize.minimize 来最小化成本函数f(r, Q)，但我面临的问题是，它输出的结果--几乎与优化.最小化参数中给出的初始猜测相同。下面是最小化函数的代码。

initial_guess = [2500., 10000.]
result = optimize.minimize(f, initial_guess, bounds=[(1, 5000), (5000, 15000)], tol=1e-3)
print(result)

输出：

      fun: 2712.7698818644253
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([-0.01195986, -0.01273293])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 6
      nit: 1
   status: 0
  success: True
        x: array([ 2500.01209628, 10000.0127784 ])

输出 x: array([ 2500.01209628, 10000.0127784 ]) - 我怀疑这是不是真正的答案，而且也和最初提供的猜测差不多。我是不是在最小化的过程中做错了什么，还是有其他方法？所以我想绘制成本函数，自己看看。

如果能有一个交互式的图，让我玩玩就好了。

Answer 1

我的回答只涉及到绘图，但最后我会对最小值的问题进行评论。

对于你所需要的三维曲面图，我认为是矫枉过正的，我将向你展示，而不是使用 contourf 和 contour 以便对你的函数有一个很好的了解。

首先，代码--关键点。

• 你的代码不能在向量环境下执行 所以我写了一个显式循环来计算值。
• 由于Matplotib的设计，矩阵数据的x轴是关联在列上的，这一点必须要考虑到。
• 的结果 countour 和 contourf 必须保存，因为标签和颜色条分别需要它们。
• 没有标签或图例，因为我不知道你在做什么。

也就是说，这里是代码

import matplotlib.pyplot as plt
import numpy as np
from numpy import sqrt, pi, exp
from scipy.integrate import quad

mean, std = 295, 250
l, m, p = 7, 30, 15
w, K, c = 7, 100, 5
h = 0.001  # per unit per  day

# defining Cumulative distribution function
def cdf(x):
    cdf_eqn = lambda t: (1 / (std * sqrt(2 * pi))) * exp(-(((t - mean) ** 2) / (2 * std ** 2)))
    cdf = quad(cdf_eqn, -np.inf, x)[0]
    return cdf
# defining Probability density function
def pdf(x):
    return (1 / (std * sqrt(2 * pi))) * exp(-(((x - mean) ** 2) / (2 * std ** 2)))
# 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
# pulling it all together
def f(r, Q):
    TC = G(r, Q)/CL(r, Q) + I(r, Q)
    return TC

nr, nQ = 6, 11
r = np.linspace(0, 5000, nr)
Q = np.linspace(5000, 15000, nQ)

z = np.zeros((nr, nQ)) # r ←→ y, Q ←→ x

for i, ir in enumerate(r):
    for j, jQ in enumerate(Q):
        z[i, j] = f(ir, jQ)
    print('%2d: '%i, ','.join('%8.3f'%v for v in z[i]))

fig, ax = plt.subplots()
cf = plt.contourf(Q, r, z)
cc = plt.contour( Q, r, z, colors='k')
plt.clabel(cc)
plt.colorbar(cf, orientation='horizontal')
ax.set_aspect(1)
plt.show()

而这里的执行结果

$ python cost.py
 0:  4093.654,3661.777,3363.220,3120.073,2939.119,2794.255,2675.692,2576.880,2493.283,2426.111,2359.601
 1:  4072.865,3621.468,3315.193,3068.710,2887.306,2743.229,2626.065,2528.934,2447.123,2381.802,2316.991
 2:  4073.852,3622.443,3316.163,3069.679,2888.275,2744.198,2627.035,2529.905,2448.095,2382.775,2317.965
 3:  4015.328,3514.874,3191.722,2939.397,2758.876,2618.292,2505.746,2413.632,2336.870,2276.570,2216.304
 4:  3881.198,3290.628,2947.273,2694.213,2522.845,2394.095,2293.867,2213.651,2148.026,2098.173,2047.140
 5:  3616.675,2919.726,2581.890,2352.015,2208.814,2106.289,2029.319,1969.438,1921.555,1887.398,1849.850
$ 

我可以补充一下，全局最小值和全局最大值在角上，而在近似区域内有两条局部最小值（下线）和局部最大值（上线）的次水平线。r ≈1000和 r ≈ 2000.