Scipy.optimize.minimize非线性方程组SLSQP：不等式约束仍得负解

问题背景 #

求解含11个变量的11个非线性联立方程组，要求所有变量为正。已设置非负不等式约束，但求解结果中仍出现接近0的负解（如Pz=-1.38777878e-17）；更换初始猜测值后，解恢复非负。需明确该现象的原因，并寻求解决方法，同时推荐适用于带约束场景的其他优化算法（如Nelder-Mead）。

求解代码 #

import numpy as np
from scipy.optimize import minimize
import sympy as sym
import math as math

### 定义参数值
Cpz= 2 # 参数
intensity=  0.04996599987353118 # 参数
length = 10 # 参数
Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz, Dx, Dy, Dz = 0, length, 0, length, length,0,  0, 0, 0, length, 0, 0
ni, nm= 1.1, 1 # 光学油脂和闪烁晶体的折射率 # 参数
thetacrit= sym.asin(nm/ni) # 弧度制
print(f"thetacrit in radians is {thetacrit}")

def my_fun(param):
    thetaC1= param[0]
    thetaC2= param[1]
    Cpx= param[2]
    Px=param[3]
    Pz=param[4]
    Cphorizdist= param[5]
    Cpvertdist= param[6]
    Cpdist3D= param[7]
    Chorizdist= param[8]
    Cvertdist= param[9]
    Cdist3D= param[10]
    f= np.zeros(11)
    
    f[0]= sym.asin((nm/ni)*sym.sin(thetaC2))- thetaC1
    f[1]= np.absolute(Px-Cpx)- Cphorizdist
    f[2]= np.absolute(Pz-Cpz)- Cpvertdist
    f[3]= ( (Cphorizdist)**2 + (Cpvertdist)**2 )**(1/2)-Cpdist3D
    f[4]= np.absolute(Cpx-Cx)- Chorizdist
    f[5]= np.absolute(Cpz-Cz)-Cvertdist
    f[6]= (Chorizdist**2 + Cvertdist**2)**(1/2)- Cdist3D 
    f[7]= Cphorizdist/Cpdist3D-sym.sin(thetaC1)
    f[8]= Cpx/Cdist3D- sym.sin(thetaC2)
    f[9]= Cphorizdist/Cpvertdist- sym.tan(thetaC1)
    f[10]= 1/((Cpdist3D+Cdist3D)**2)-intensity
    
    return np.dot(f,f) # 最小化残差平方和

def my_cons(param):
    thetaC1= param[0]
    thetaC2= param[1]
    Cpx= param[2]
    Px=param[3]
    Pz=param[4]
    Cphorizdist= param[5]
    Cpvertdist= param[6]
    Cpdist3D= param[7]
    Chorizdist= param[8]
    Cvertdist= param[9]
    Cdist3D= param[10]

    f = np.zeros(13)
  
    # SLSQP的不等式约束要求返回值 >=0
    f[0]= thetaC1
    f[1]= math.pi/2-thetaC1
    f[2]= thetaC2
    f[3]= math.pi/2-thetaC2
    f[4]= Cpx
    f[5]= Px
    f[6]= Pz
    f[7]= Cphorizdist
    f[8]= Cpvertdist
    f[9]= Cpdist3D
    f[10]= Chorizdist
    f[11] = Cvertdist
    f[12] = Cdist3D
    return f

cons = {'type' : 'ineq', 'fun': my_cons}
res = minimize(my_fun, (0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5), method='SLSQP',\
           constraints=cons,options={"maxiter": 5000})
print(f"求解结果：{res}")

原因分析 #

1. 浮点数值精度误差：得到的-1.38777878e-17本质是接近0的计算噪声，远小于常规数值精度阈值（如1e-10）。SLSQP算法默认的约束容忍度（ftol）为1e-6，允许约束存在微小违反，因此会出现这种边界附近的数值偏差。
2. 初始猜测的影响：优化算法的收敛路径依赖初始点，部分初始点可能引导算法进入约束边界的数值不稳定区域，而合适的初始点会让算法收敛到严格满足约束的解。
3. 非线性方程组的特性：目标函数在约束边界处的梯度变化可能不明显，导致算法难以精确停留在非负区域，出现微小的数值越界。

解决方法 #

1. 调整优化器精度参数：在minimize的options中设置更严格的容忍度，比如：

   options={"maxiter": 5000, "ftol": 1e-10, "eps": 1e-12}
   

   缩小ftol（约束容忍度）和eps（有限差分的步长），让算法对约束违反更敏感。
2. 变量替换实现严格非负：将所有非负变量替换为正数的平方，例如Pz = pz_sq**2，其中pz_sq为实数变量。这样从数学上保证变量非负，彻底避免数值越界问题。修改后的目标函数需要对应替换变量，并重新计算残差。
3. 添加约束惩罚项：在目标函数中对违反非负约束的情况施加惩罚，例如：

   penalty = 0
   for var in param:
       if var < 0:
           penalty += 1e10 * abs(var) # 对负变量施加极大惩罚
   return np.dot(f,f) + penalty
   

   强制算法远离负区域。
4. 改用支持边界约束的算法：比如L-BFGS-B，无需定义复杂的约束函数，直接通过bounds参数设置变量下限为0：

   bounds = [(0, math.pi/2), (0, math.pi/2)] + [(0, None)]*9 # 前两个变量有上下限，其余仅下限0
   res = minimize(my_fun, (0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5), method='L-BFGS-B', bounds=bounds, options={"maxiter":5000})
   

推荐的带约束优化算法 #

• L-BFGS-B：最适合单纯的边界约束（如非负）场景，属于拟牛顿法，效率高，无需手动定义约束函数，直接设置bounds即可。
• Nelder-Mead：无导数优化算法，本身不支持约束，但可通过可行域内初始化+目标函数惩罚项实现约束优化。例如，当变量为负时，目标函数返回极大值，迫使算法在可行域内搜索。
• Trust-Region-Reflective：适用于带边界约束的优化，需要目标函数提供梯度（可通过scipy.optimize.approx_fprime计算数值梯度），收敛稳定性较好。

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