SHGO优化返回None问题：蒙特卡洛收益变异性影响及解决方案咨询

我正在最小化目标函数：

-100 * fill_rate(floors) + l1 * ||epsilon|| + l2 * ||gamma||

约束条件包括：

1. 收益约束：

0.95 * budget ≤ revenue(floors) ≤ budget
（收益函数通过蒙特卡洛采样计算，存在变异性）

2. 参与率约束：

0 ≤ participation rate(floor) + gamma(floor) ≤ 1

同时floors、epsilon、gamma均为3维数组（epsilon和gamma为松弛变量），且各自存在取值边界。已知至少存在一个最优解，但使用SHGO优化时输出为None，尝试过sobol、halton、simplicial三种采样方法均无效，而COBYLA和BOBYQA优化器可正常运行。请问是否是收益评估的变异性导致SHGO失效？恳请提供让SHGO正常工作的建议。

import time
import numpy as np
from scipy.stats import norm, bernoulli
from shgo import shgo

"""
Supply = BQ1 U BQ2 U BQ3
"""
# means of the bqs
mean1 = 5
mean2 = 6
mean3 = 3
# variance of bqs
std1 = 5
std2 = 4
std3 = 1
# probabilities that bqs wont participate in the auction due to some reasons
err1 = 0.6
err2 = 0.7
err3 = 0.2
# current floors of each bqs
# do the auction and collect data to get ao, fr, cpm


def get_bq_samples(pr, loc, scale, size):
    """
    get real samples of bqs when they participate in the auction and collect the data. The real distributions of bqs are
    unknown to us. Once we have collected data, we will use that to evaluate our objective function
    """
    samples1 = norm.rvs(loc=loc, scale=scale, size=size)
    samples2 = bernoulli.rvs(p=pr, size=size)
    return np.multiply(samples1, samples2)


def simulate_auction_collect_bids(bq_budgets, floors, num_samples):
    """
    Simulate the first price auction for 3 advertisers(bqs) and collect the bids data
    """
    X1, X2, X3 = [], [], []
    R1, R2, R3 = 0, 0, 0
    bids_above_floor1 = 0
    bids_above_floor2 = 0
    bids_above_floor3 = 0
    samples1 = get_bq_samples(err1, mean1, std1, num_samples)
    samples2 = get_bq_samples(err2, mean2, std2, num_samples)
    samples3 = get_bq_samples(err3, mean3, std3, num_samples)

    bq1_budget_exhausted = False
    bq2_budget_exhausted = False
    bq3_budget_exhausted = False
    served_ao = 0
    supply_rev = 0
    for b1, b2, b3 in zip(samples1, samples2, samples3):
        if bq1_budget_exhausted and bq2_budget_exhausted and bq3_budget_exhausted:
            print("Budget exhausted but inventory is there")
            X1.append(0)
            X2.append(0)
            X3.append(0)
            continue
        ab1 = b1
        ab2 = b2
        ab3 = b3
        bids_above_floor1 += 1
        bids_above_floor2 += 1
        bids_above_floor3 += 1
        if not bq1_budget_exhausted:
            X1.append(b1)
        else:
            X1.append(0)
        if not bq2_budget_exhausted:
            X2.append(b2)
        else:
            X2.append(0)
        if not bq3_budget_exhausted:
            X3.append(b3)
        else:
            X3.append(0)
        if (b1 < floors[0]) or (bq1_budget_exhausted == True):
            ab1 = 0
            bids_above_floor1 -= 1
        if (b2 < floors[1]) or (bq2_budget_exhausted == True):
            ab2 = 0
            bids_above_floor2 -= 1
        if (b3 < floors[2]) or (bq3_budget_exhausted == True):
            ab3 = 0
            bids_above_floor3 -= 1
        actual_bids = np.array([ab1, ab2, ab3])
        Y = np.amax(actual_bids)
        Z = np.argmax(actual_bids)
        if Y > 0:
            served_ao += 1
            supply_rev += Y
            if Z == 0:
                R1 += Y
            elif Z == 1:
                R2 += Y
            elif Z == 2:
                R3 += Y
            if R1 >= bq_budgets[0]:
                bq1_budget_exhausted = True
            if R2 >= bq_budgets[1]:
                bq2_budget_exhausted = True
            if R3 >= bq_budgets[2]:
                bq3_budget_exhausted = True

    supply_fr = served_ao / num_samples
    bq1_fr = bids_above_floor1 / num_samples
    bq2_fr = bids_above_floor2 / num_samples
    bq3_fr = bids_above_floor3 / num_samples
    return supply_rev, supply_fr, [R1, R2, R3], [bq1_fr, bq2_fr, bq3_fr], np.array([X1, X2, X3])


def get_participation_rate(max_bids, floors):
    """
    Given bids data, get the participation rate for all bqs
    """
    N = len(max_bids[0])
    pr = []
    for bq_bids, f in zip(max_bids, floors):
        rate = np.count_nonzero(np.where(bq_bids < f, 0, bq_bids))
        rate = rate / N
        pr.append(rate)
    return pr


def get_delta_samples(max_bids, floors, gamma, num_samples):
    """
    Given bids data, get the bernoulli samples with the participation rate determined by the bids observed
    """
    pr_arr = get_participation_rate(max_bids, floors)
    temp = []
    for pr, g in zip(pr_arr, gamma):
        if 0 <= pr+g <= 1:
            temp.append(bernoulli.rvs(p=pr+g, size=num_samples))
        else:
            print("Optimisation not respecting the gamma constraint!")
            t = max(min(pr+g, 1),0)
            temp.append(bernoulli.rvs(p=0, size=num_samples))
    return np.array(temp)


def get_max_bids_samples(max_bids, floors, epsilon, num_samples):
    """
    Get max-bids samples using the observed bids data
    """
    temp = []
    for bq_bids, f, e in zip(max_bids, floors, epsilon):
        valid_bids = bq_bids[bq_bids >= f]
        valid_bids = np.add(valid_bids, e)
        temp.append(np.random.choice(valid_bids, size=num_samples))
    return np.array(temp)


def get_revenue(bids, X, ao, num_samples, truncate=True):
    """
    Get the revenue as a function of floor using the observed bids data. We need this into constraint equation
    """
    floors = X[:3]
    epsilon = X[3:6]
    gamma = X[6:9]
    max_bid_samples = get_max_bids_samples(bids, floors, epsilon, num_samples)
    delta_samples = get_delta_samples(bids, floors, gamma, num_samples)
    X = np.multiply(max_bid_samples, delta_samples)
    Y = np.amax(X, axis=0)
    Z = np.argmax(X, axis=0)
    R1 = np.sum(Y[np.argwhere(Z == 0)])
    R2 = np.sum(Y[np.argwhere(Z == 1)])
    R3 = np.sum(Y[np.argwhere(Z == 2)])
    if truncate:
        return [round(np.mean(Y), 3) * ao, round(R1 / num_samples, 3) * ao, round(R2 / num_samples, 3) * ao,
                round(R3 / num_samples, 3) * ao]
    else:
        return [np.mean(Y) * ao, R1 * ao / num_samples, R2 * ao / num_samples, R3 * ao / num_samples]


def get_fr(bids, X):
    """
    Get fill rate as a function of floor using the observed bids data. fr is part of optimisation objective function
    """
    floors = X[:3]
    gamma = X[6:9]
    prs = get_participation_rate(bids, floors)
    fr = 1
    for pr, g in zip(prs, gamma):
        if 0 <= pr + g <= 1:
            fr = fr * (1 - (pr+g))
        else:
            t = max(min(pr+g, 1),0)
            t=0
            fr = fr*(1-t)
    return 1-fr


def l2_norm(arr):
    out = 0
    for a in arr:
        out += a * a
    return out/len(arr)


def find_optimum_floors_shgo(bids, budget, num_samples):
    """
    We want to minimise out objective function which is -100 * fr + l1 * ||epsilon|| + l2 * ||gamma||, subject to revenue
    constraint:    0.95 * budget <= revenue(floor) <= budget
    and 0 <= participation rate + gamma <= 1
    together with bounds on floors, epsilon and gamma
    """
    # X: input vector for optimisation [floor1, floor2, floor3, epsilon1, epsilon2, epsilon3, gamma1, gamma2, gamma3]
    # bounds contains the bounds for X vector
    AO = len(bids[0])
    bounds = [(2.5, 3.5), (3.5, 4.5), (1.0, 2.0), (-0.15, 0.15), (-0.15, 0.15), (-0.15, 0.15), (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1)]

    # l1 l2 are hyperparameters for the objective function
    l1 = 50
    l2 = 500

    def objective_function(X):
        # objective function is the fill rate minus penalties on epsilon and gamma
        ep = X[3:6]
        gm = X[6:9]
        fr = get_fr(bids, X)
        return -100*fr + l1 * l2_norm(ep) + l2 * l2_norm(gm)

    # budget constraint: we want the revenue to remain between 0.95 of budget and budget
    def get_lower_budget_constraint(X, i):
        rev_arr = get_revenue(bids, X, AO, num_samples, False)
        return rev_arr[1:][i] - 0.95 * budget[i]

    def get_upper_budget_constraint(X, i):
        rev_arr = get_revenue(bids, X, AO, num_samples, False)
        return budget[i] - rev_arr[1:][i]

    # participation rate constraint: we want participation rate to remain between 0.1 and 0.8
    def get_lower_delta_constraint(X, i):
        flr = X[:3]
        gm = X[6:9]
        pr = np.add(get_participation_rate(bids, flr),gm)
        return pr[i]

    def get_upper_delta_constraint(X, i):
        flr = X[:3]
        gm = X[6:9]
        pr = np.add(get_participation_rate(bids, flr), gm)
        return 1-pr[i]

    constraints = []
    lower_delta_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_lower_delta_constraint(x, i)} for i in
                              range(3)]
    constraints.extend(lower_delta_constraint)
    upper_delta_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_upper_delta_constraint(x, i)} for i in
                              range(3)]
    constraints.extend(upper_delta_constraint)
    lower_budget_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_lower_budget_constraint(x,i)} for i in range(3)]
    constraints.extend(lower_budget_constraint)
    upper_budget_constraint = [{'type': 'ineq', 'fun': lambda x, i=i: get_upper_budget_constraint(x, i)} for i in range(3)]
    constraints.extend(upper_budget_constraint)

    res = shgo(func = objective_function, bounds = bounds, sampling_method='sobol', constraints=constraints)

    print(res.x)
    print(res.fun)

    if res.x is not None:
        print(get_revenue(bids, res.x, AO, num_samples, False))
        print(get_fr(bids, res.x))
        print(res.x[:3])


if __name__ == '__main__':
    ao = 10000
    budgets = np.array([25000, 35000, 2000])
    initial_floors = np.array([3, 4, 1.5])


    # do the simulate and collect the data
    s_rev, s_fr, bq_rev, bq_fr, all_bids = simulate_auction_collect_bids(budgets, initial_floors, ao)
    print(f"Supply revenue: {s_rev}")
    print(f"Supply fill rate: {s_fr}")
    print(f"Each BQ revenue: {bq_rev}")
    print(f"Each BQ fill rate: {bq_fr}")

    # one solution satisfying all constraints
    initial_X = np.append(initial_floors, [0, -0.1, -0.1, 0.0, 0.03, 0.018])
    print(f"One solution with X: {initial_X}")
    revenue_init_flr = get_revenue(all_bids, initial_X, ao, 100000, truncate=False)
    print(np.divide(np.array(revenue_init_flr[1:]), budgets) * 100)
    fr_init_flr = get_participation_rate(all_bids, initial_floors)
    print(fr_init_flr)

    #finding optimum via shgo
    print("Starting shgo optimisation, to find optimum floors")
    find_optimum_floors_shgo(all_bids, budgets, 100000)

一、收益评估变异性对SHGO的影响 #

是的，收益函数的蒙特卡洛变异性是SHGO失效的核心原因：

• SHGO是确定性全局优化器，依赖一致的函数值和约束判断来规划搜索路径。而你的get_revenue每次调用都会重新随机采样，返回值存在波动，会导致SHGO无法准确判断当前点是否满足约束、是否为更优解，最终因逻辑混乱返回None。
• COBYLA和BOBYQA属于局部优化器，对函数噪声的容忍度更高，迭代逻辑更依赖局部趋势而非全局一致的函数值，因此能正常运行。

二、让SHGO正常工作的具体建议 #

1. 消除收益函数的随机性 #

• 固定随机种子：在所有涉及随机采样的函数开头固定numpy种子，确保同一输入对应完全一致的输出：

  def get_revenue(bids, X, ao, num_samples, truncate=True):
      np.random.seed(42)  # 固定种子，消除随机性
      # 原有代码...
  

  注意：固定种子会牺牲采样的统计随机性，但这是SHGO能正常工作的必要前提。如果需要保留统计准确性，可提前对每个输入多次采样取均值，将结果固化为查询表，优化时直接查表。

• 预计算收益均值：提前对可行域内的关键点多次运行蒙特卡洛采样，计算收益的均值并存储，优化时直接调用预计算结果，避免实时采样的波动。

2. 调整SHGO参数提升稳定性 #

• 增加初始采样点数量：SHGO的全局搜索依赖初始采样覆盖可行域，9维问题可将n参数设为2000（默认是100*维度）：

  res = shgo(func=objective_function, bounds=bounds, sampling_method='sobol', constraints=constraints, n=2000)
  

• 放宽约束容忍度：