如何用Scipy.minimize结合ε数值求解函数上确界

我需要数值求解如下函数的上确界：

其中p(x)是满足p(0+)>0的单调递减概率函数，p(x)的示例曲线为类似0.5 - 0.5*erf((10/√2)*(α/σ)*log₁₀(x/200))的形式（x范围0到600）。

我尝试了以下方法，但均存在问题：

• 使用scipy.optimize.minimize求解最大值：ε取值难以适配所有p(x)实例，ε过小会导致minimize返回无穷大。
• 自行编写循环脚本：部分场景表现较好，但当σ/α趋近0时效果不佳。
• 通过二阶导数找极值点：结果不符合需求。

请问我的思路存在什么问题？有没有其他可行方法？

附相关代码 #

基于scipy.minimize的实现 #

import scipy as sp
import math
import numpy as np

def density_upper_scaled(alpha, sigma):
    epsilon = 0.1
    constraint = {'type': 'ineq', 'fun': constraint_function, 'args': (alpha, sigma)}
    guess = np.array([1])
    result = sp.optimize.minimize(density_upper_function, guess, bounds=[(epsilon, None)], constraints=constraint)
    max_value = -result.fun  

    return (1.437 * 200**2 / max_value)
#, {'type': 'ineq', 'fun': lambda x: epsilon * x**2}

def density_upper_function(r):
    return - r**2 * 0.1


def constraint_function(r, alpha, sigma):
    return probability_function(r, alpha, sigma) - 0.1


def probability_function(distance, alpha, sigma):  
    ro0 = 200    # theoretical communication distance for sigma = 0
    if sigma != 0:
        p = 0.5 - 0.5 * math.erf((10 / math.sqrt(2)) * (alpha / sigma) * math.log10(distance / ro0))
    else:
        p = int(distance <= ro0)
    return p

自行编写的循环脚本 #

condition = probability_function(x, alpha, sigma)
prior = -1
epsilon = 0.2
while condition >= epsilon and condition * x ** 2 >= epsilon:
    prior = condition
    condition = round(probability_function(x, alpha, sigma),5)

    x += 1
print (x, condition)
print(x,condition, 200**2 * 1.437 / (condition * x **2) )

问题分析与解决方案 #

现有思路的核心问题 #

1. scipy.minimize实现的致命错误

   • 目标函数density_upper_function写死了系数0.1，完全偏离了最大化x²p(x)的核心需求，应该直接返回-x²*probability_function(x, alpha, sigma)，让minimize求解等价的最小值问题。
   • 硬编码ε=0.1，没有根据p(x)的特性动态调整：当σ/α趋近0时，p(x)趋近于阶跃函数（x≤200时p=1，x>200时p=0），固定ε会导致优化器找不到有效可行域，直接返回无穷大。

2. 循环脚本的局限性

   • 固定步长为1+截断概率值的逻辑，会在σ/α趋近0时错过x=200附近的极值点；终止条件未直接针对最大化x²p(x)的目标，容易提前退出循环。

3. 二阶导数方法的失效场景

   • p(x)包含erf和对数函数，二阶导数解析形式复杂；当σ/α趋近0时，p(x)变为不可导的阶跃函数，二阶导数方法完全失效。

可行解决方案 #

1. 修正scipy.minimize实现

• 修正目标函数：直接映射到最大化x²p(x)的等价问题

  def density_upper_function(x, alpha, sigma):
      return - (x**2 * probability_function(x, alpha, sigma))
  

• 动态调整ε与约束：根据σ/α的比例设置极小值ε，避免硬编码；σ=0时直接返回x=200处的结果（阶跃函数的极值点）。
• 优化初始猜测：将初始值设为200（ro0），让优化器更快收敛到极值区域。

修正后的示例代码：

import scipy as sp
import math
import numpy as np

def probability_function(distance, alpha, sigma):  
    ro0 = 200
    if sigma != 0:
        p = 0.5 - 0.5 * math.erf((10 / math.sqrt(2)) * (alpha / sigma) * math.log10(distance / ro0))
    else:
        p = 1.0 if distance <= ro0 else 0.0
    return p

def density_upper_function(x, alpha, sigma):
    return - (x**2 * probability_function(x, alpha, sigma))

def density_upper_scaled(alpha, sigma):
    if sigma == 0:
        return 1.437  # 直接返回x=200处的计算结果
    
    ratio = sigma / alpha
    epsilon = 1e-8 if ratio > 1e-3 else 1e-4  # 动态调整ε
    
    def constraint(x):
        return probability_function(x, alpha, sigma) - epsilon
    
    bounds = [(1e-5, None)]
    guess = np.array([200.0])
    result = sp.optimize.minimize(density_upper_function, guess, args=(alpha, sigma), 
                                 bounds=bounds, constraints={'type': 'ineq', 'fun': constraint},
                                 method='L-BFGS-B')
    
    if result.success:
        max_value = -result.fun
        return (1.437 * 200**2) / max_value
    else:
        # 优化失败时 fallback 到x=200处的结果
        return (1.437 * 200**2) / (200**2 * probability_function(200, alpha, sigma))

2. 黄金分割法（单峰函数专用一维搜索）

由于f(x)=x²p(x)是单峰函数（先增后减），黄金分割法无需计算导数，稳定性更强，适合非光滑或导数复杂的场景：

def golden_section_search(alpha, sigma, tol=1e-5):
    ro0 = 200
    
    # 自动确定搜索右端点
    def find_right_bound():
        x = ro0
        while probability_function(x, alpha, sigma) > 1e-8:
            x *= 1.5
        return x
    
    a = 1e-5
    b = find_right_bound()
    gr = (math.sqrt(5) - 1) / 2  # 黄金分割比例
    
    c = b - gr * (b - a)
    d = a + gr * (b - a)
    fc = c**2 * probability_function(c, alpha, sigma)
    fd = d**2 * probability_function(d, alpha, sigma)
    
    while abs(c - d) > tol:
        if fc > fd:
            b = d
            d = c
            fd = fc
            c = b - gr * (b - a)
            fc = c**2 * probability_function(c, alpha, sigma)
        else:
            a = c
            c = d
            fc = fd
            d = a + gr * (b - a)
            fd = d**2 * probability_function(d, alpha, sigma)
    
    max_x = (a + b) / 2
    max_value = max_x**2 * probability_function(max_x, alpha, sigma)
    return (1.437 * 200**2) / max_value

这种方法无需设置ε，在σ/α趋近0时会自动收敛到x=200附近，稳定性优于梯度优化器。

内容的提问来源于stack exchange，提问作者Sinured