Scipy优化失败：为何续保倾向模型无法收敛？

问题概述 #

执行利润最大化的客户保费优化任务，需满足整体续保率目标约束，但优化流程失败，运行result返回奇异矩阵错误，且最优保费与初始值完全一致。

错误输出 #

fun: -298.2633664479155
     jac: array([-0.82517242, -0.60799408, -0.44937897, -0.82517242, -0.72691727,
       -0.84170532, -0.74860764, -0.18024063, -0.74860764, -0.96907043])
 message: 'Singular matrix C in LSQ subproblem'
    nfev: 11
     nit: 1
    njev: 1
  status: 6
 success: False
       x: array([338.39, 409.57, 453.83, 337.21, 377.76, 330.7 , 366.62, 559.69,
       366.95, 156.04])

Jacobian矩阵符合预期（保费上升续保率下降，不同保单变化有差异）。

代码重现 #

import numpy as np
import pandas as pd
import scipy.optimize as opt

# 数据
cost_prices = np.array([293.2 , 371.08, 479.59, 328.2 , 295.23, 273.17, 303.83, 368.09,
        229.72, 208.64])
initial_premiums = ([338.39, 409.57, 453.83, 337.21, 377.76, 330.7 , 366.62, 559.69,
        366.95, 156.04])

# 目标函数
def objective(premiums, cost_prices):
    propensities = propensity_to_renew(premiums)
    profit = np.sum((premiums - cost_prices) * propensities)
    return -profit  # Scipy默认最小化，因此取负

# 分段续保倾向模型
Lower = [0, 215,    231.7,  248.4,  265.1,  281.8,  298.5,  315.2,  331.9,  348.6,  365.3,  382,    398.7,  415.4,  432.1,  448.8,  465.5,  482.2,  498.9,  515.6,  532.3,  549,    565.7,  582.4,  599.1,  615.8,  632.5,  649.2,  665.9,  682.6,  699.3,  716,    732.7,  749.4,  766.1,  782.8,  799.5,  816.2,  832.9,  849.6,  866.3,  883,    899.7,  916.4,  933.1,  949.8,  966.5,  983.2,  999.9,  1016.6, 1033.3, 1050]
Upper = [214.999,   231.7,  248.4,  265.1,  281.8,  298.5,  315.2,  331.9,  348.6,  365.3,  382,    398.7,  415.4,  432.1,  448.8,  465.5,  482.2,  498.9,  515.6,  532.3,  549,    565.7,  582.4,  599.1,  615.8,  632.5,  649.2,  665.9,  682.6,  699.3,  716,    732.7,  749.4,  766.1,  782.8,  799.5,  816.2,  832.9,  849.6,  866.3,  883,    899.7,  916.4,  933.1,  949.8,  966.5,  983.2,  999.9,  1016.6, 1033.3, 1050,   1000000]
COMBINED = [0,  -0.07492837,    -0.222324113125,    -0.366118206625,    -0.506502842875,    -0.64346552175, -0.77682655125, -0.906585931125,    -1.03292335375, -1.155851319,   -1.27517763475, -1.39108199325, -1.50357689425, -1.612457645875,    -1.717749248125,    -1.819798585,   -1.91807908025, -2.013117310375,    -2.104553891,   -2.192376322125,    -2.276801796,   -2.357805312375,    -2.4352071795,  -2.509174589125,    -2.579577849375,    -2.646521652,   -2.7098888055,  -2.769834001625,    -2.826369740125,    -2.879291329375,    -2.92881596125, -2.974713943625,    -3.017214968625,    -3.05629403625, -3.0917589545,  -3.123814415375,    -3.15225572675, -3.17731258075, -3.19875528525, -3.216776032625,    -3.231195130375,    -3.24217977075, -3.24977995375, -3.253765987375,    -3.254342563625,    -3.251304990375,    -3.24484545975, -3.234988971875,    -3.2215183345,  -3.204625739625,    -3.1841564955,  -3.17255777875]

elas = pd.DataFrame(
    {'Lower': Lower,
     'Upper': Upper,
     'COMBINED': COMBINED
    })

def propensity_to_renew(premium):
    prop_vec = []
    for prem in premium:
        val = prem
        lp = float(elas.query(' @val >= Lower and @val < Upper')['COMBINED'])
        prop_vec.append(1/(1+np.exp(-(-1.496+lp))))
    return prop_vec

# 续保率约束
def retention_constraint(premiums):
    propensities = propensity_to_renew(premiums)
    target_retention = 0.07  # 示例目标续保率
    return np.mean(propensities) - target_retention

# 约束与边界
constraints = [{'type': 'eq', 'fun': retention_constraint}]
bounds = [(0.7*premium,1.3*premium) for premium in initial_premiums]

# 执行优化
result = opt.minimize(objective, initial_premiums, args=(cost_prices), bounds=bounds, constraints=constraints)

optimal_premiums = result.x
optimal_profit = -result.fun

print("Optimal Premiums:", optimal_premiums)
print("Optimal Profit:", optimal_profit)

排查发现 #

替换为连续可导的logistic模型时，优化可正常完成：

def propensity_to_renew(premium):
    return 1 / (1 + np.exp(0.01 * (premium - 100)))

两种模型的续保倾向随保费变化趋势高度相似：均随保费上升单调下降，整体呈logistic曲线形态。

解决建议 #

1. 根本原因分析 #

原分段续保倾向模型是分段常数函数，在分段边界处不可导，而Scipy默认的SLSQP优化器依赖目标函数和约束函数的光滑性（连续可导）。非光滑点会导致雅可比矩阵计算异常，进而出现“奇异矩阵C”的错误，优化器无法找到有效搜索方向，直接返回初始值。

2. 具体修复方案 #

方案一：光滑化分段模型

将分段常数的COMBINED值插值为连续可导的函数，消除非光滑点：

from scipy.interpolate import CubicSpline

# 用区间中点作为插值节点
mid_points = [(l + u)/2 for l, u in zip(Lower[:-1], Upper[:-1])]
# 三次样条插值生成连续可导的lp函数
cs = CubicSpline(mid_points, COMBINED[:-1])

def propensity_to_renew_smoothed(premium):
    premium = np.asarray(premium)
    lp = cs(premium)
    # 处理超出插值范围的保费，使用边界值
    lp[premium < mid_points[0]] = COMBINED[0]
    lp[premium > mid_points[-1]] = COMBINED[-1]
    return 1/(1+np.exp(-(-1.496 + lp)))

替换原propensity_to_renew函数后，SLSQP优化器可正常运行。

方案二：更换支持非光滑优化的算法

使用COBYLA算法（无需计算雅可比矩阵），需将边界条件转化为不等式约束：

# 将边界转化为COBYLA兼容的不等式约束
bounds_constraints = []
for idx, (low_bound, high_bound) in enumerate(bounds):
    # 下限约束：x[idx] >= low_bound
    bounds_constraints.append({
        'type': 'ineq',
        'fun': lambda x, idx=idx, low=low_bound: x[idx] - low
    })
    # 上限约束：x[idx] <= high_bound
    bounds_constraints.append({
        'type': 'ineq',
        'fun': lambda x, idx=idx, high=high_bound: high - x[idx]
    })

# 合并续保率约束与边界约束
total_constraints = constraints + bounds_constraints

# 使用COBYLA执行优化
result = opt.minimize(
    objective, 
    initial_premiums, 
    args=(cost_prices), 
    method='COBYLA', 
    constraints=total_constraints,
    options={'maxiter': 1000}  # 可适当调整迭代次数
)

方案三：微调初始点避免边界问题

检查初始保费是否刚好落在Lower/Upper的分界点上，数值精度问题可能导致查询异常。可对初始值进行微小调整（如±0.01），避免落在分段边界。

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