python - 如何绘制多元积分函数的 3D 图，并找到其全局最小值

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

如何在 3D 图中绘制此函数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 ])- 我怀疑这是真正的答案，也几乎与提供的初始猜测相同。我在最小化方面做错了什么还是有其他方法可以做到这一点？所以我想绘制成本函数并自己寻找。

如果我能有一个互动的情节来玩，那就太好了

标签： pythonnumpymatplotlibscipyplotly

Answer 1

我的回答只与绘图有关，但最后我会评论极小极大的问题。

对于您需要的 3D 曲面图，恕我直言，矫枉过正，我将向您展示如何使用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中有两条局部最小值（下线）和局部最大值（上线）的亚水平线。