使用scipy.optimize.minimize求解优化问题的异常解排查

问题背景 #

尝试用Scipy实现某优化问题（原问题对应图示），已知收益率数组：

r = [0.96366965, 0.93341242, 0.90676733, 0.88186071, 0.8582291,  0.83472442, 0.80977363, 0.783683, 0.75743122, 0.73259614]

当前实现代码如下：

import numpy as np
from scipy.optimize import minimize, Bounds

def objective_function(x):
    e_t = np.zeros(10)
    e_t[0] = 1 - 1 / (1 + r[0])

    for t in range(1, 9):
        sub = 0
        for ti in range(0, t):
            sub += np.sum(x[: (ti + 1)]) / (1 + r[ti]) ** (ti + 1)

        e_t[t] = 1 - 1 / (1 + r[t]) ** (t + 1) - sub

    e_t[9] = 0

    return np.sum(e_t**2)

def constraint_function(t):
    def calculate_et(x):
        sub = 0
        for ti in range(0, t):
            sub += np.sum(x[: (ti + 1)]) / (1 + r[ti]) ** (ti + 1)

        e_t = 1 - 1 / (1 + r[t]) ** (t + 1) - sub
        return e_t

    return calculate_et

x0 = np.zeros(10)
cons = (
    {"type": "ineq", "fun": constraint_function(1)},
    {"type": "ineq", "fun": constraint_function(2)},
    {"type": "ineq", "fun": constraint_function(3)},
    {"type": "ineq", "fun": constraint_function(4)},
    {"type": "ineq", "fun": constraint_function(5)},
    {"type": "ineq", "fun": constraint_function(6)},
    {"type": "ineq", "fun": constraint_function(7)},
    {"type": "ineq", "fun": constraint_function(8)},
    {"type": "ineq", "fun": constraint_function(9)},
)

result = minimize(
    objective_function,
    x0,
    bounds=Bounds(lb=0),
    constraints=cons,
)
alpha = result.x

c = np.zeros(10)
c[0] = alpha[0]
for i in range(1, 10):
   c[i] = alpha[i-1] + alpha[i]

运行后得到的c数组为：

[9.02212140e-01 9.02212140e-01 2.39808173e-13 1.06359366e-13
4.57411886e-14 1.75415238e-14 6.16173779e-15 2.17534324e-15
4.67095103e-16 1.25975522e-17]

预期的c数组更接近：

[0.0318, 0.0318, 0.0318, 0.0318, 0.0318, 0.0318, 0.0320, 0.0325, 0.0336, 0.0353]

当前代码的核心问题 #

1. 目标函数逻辑错误
   • 第一个误差项e_t[0]的计算完全独立于优化变量x，是固定值，无法被优化，导致目标函数的优化方向偏离原问题要求。
   • 循环计算e_t时的嵌套求和逻辑错误，未正确关联优化变量x与每期的定价误差，导致误差项的计算不符合原问题的数学定义。
2. 约束条件与目标函数不一致
   • 约束函数的计算式和目标函数中的e_t一致，但约束类型设为ineq（要求返回值≥0），结合平方和最小化的目标，优化器会倾向于让后续e_t尽可能为0，最终得到不符合预期的极端解。
3. 变量映射关系错误
   • 从alpha到c的转换逻辑c[i] = alpha[i-1] + alpha[i]不符合原问题中变量的定义，导致最终输出的c数组完全偏离预期。
4. 初始值与优化算法选择不佳
   • 初始值x0 = np.zeros(10)过于极端，容易让优化器陷入局部最优解；未指定优化算法，默认的SLSQP在处理非凸约束问题时，可能因初始值问题收敛到错误解。

可行解决思路 #

1. 修正目标函数与约束的逻辑 #

• 重新梳理原优化问题的数学定义，确保所有误差项e_t的计算都依赖优化变量x。例如，原问题的e_t应为债券定价误差，正确的计算逻辑应为：

  def objective_function(x):
      e_t = np.zeros(10)
      cum_x = np.cumsum(x)  # 计算累积和，对应每期的现金流总和
      discount = 1 / (1 + r) ** np.arange(1, 11)  # 贴现因子数组
      for t in range(10):
          sum_term = np.sum(cum_x[:t+1] * discount[:t+1])
          e_t[t] = 1 - sum_term
      e_t[9] = 0  # 按原问题要求保留最后一项为0
      return np.sum(e_t ** 2)
  

• 修正约束条件：确保约束的计算式与原问题要求一致，若约束为e_t ≥ 0，则约束函数应返回正确的e_t值，且与目标函数的误差定义统一。

2. 修正变量映射关系 #

• 直接将优化变量x作为c数组，或根据原问题的数学定义建立正确的映射关系，删除错误的相邻项相加逻辑。

3. 优化初始值与算法选择 #

• 使用更合理的初始值：比如用预期解的近似值x0 = np.full(10, 0.03)作为初始点，避免优化器陷入零值附近的局部最优。
• 指定更适合的优化算法：例如使用trust-constr算法，它在处理复杂约束和非凸问题时稳定性更强；或调整SLSQP的参数（如增大maxiter、减小ftol）。

4. 简化计算逻辑，避免嵌套循环 #

• 用NumPy向量化操作替代嵌套循环，提升计算效率的同时减少逻辑错误，如上面目标函数示例中的cum_x和discount数组的使用。

5. 验证逻辑正确性 #

• 在运行优化前，手动代入预期的c数组计算目标函数值和约束值，验证逻辑是否符合原问题要求；同时检查初始值对应的输出，确认计算逻辑无误。

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