Scipy minimize未正确求解约束优化问题（显示成功但结果异常）

问题定义 #

优化问题如下（其中sig₁=2.56008479，sig₂=0.24800215，c=2.098024636032762）：

$$
\min_{x_1, x_2} \ (x_1 + x_2 - c)^2 \
\text{subject to:} \
\frac{x_1}{\text{sig}_1 - x_1^2} = \frac{x_2}{\text{sig}_2 - x_2^2} \
x_1, x_2 \geq 0
$$

当前实现代码 #

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

n = 2
c = 2.098024636032762
sigma_eigs = np.array([2.56008479, 0.24800215])

def equations(x):
    eq = np.sum(x) - c
    return np.sum(eq**2)

# 初始猜测
x0 = np.ones(n) 

# 约束
bounds = Bounds([0, 0])  # 仅设置下限，未设置上限
non_linear_eq = lambda x: x[0]/(sigma_eigs[0]-x[0]**2) - x[1]/(sigma_eigs[1]-x[1]**2)
constr1 = NonlinearConstraint(non_linear_eq, 0, 0)

result = minimize(equations, x0, method='trust-constr', constraints=constr1, bounds=bounds)

x = result.x
print(x)

const_result = [x[i]/(sigma_eigs[i]-x[i]**2) for i in range(n)]
print("The solution is x = {}".format(x))
if np.any(x < 0):
    print("Error: some of the x values are non-positive")
else:
    print("All x values are positive")

residuals = equations(x)
print("Residuals are: {}".format(residuals))
if np.any(np.abs(residuals) > 1e-6):
    print("Error: The solution is not accurate. The residuals are high")
else:
    print("All residuals are within the tolerance")

求解问题 #

使用trust-constr方法得到的结果：

[-2.75911573e-08  1.22830289e+07]

输出信息：

x = [-2.75911573e-08  1.22830289e+07]
Error: some of the x values are non-positive
Residuals are: 150872747356720.12
Error: The solution is not accurate. The residuals are high

优化提示：

message: xtol termination condition is satisfied.

当前问题表现为变量发散、不满足非负约束，且残差极大；使用其他求解器会得到如[0. 0.]这类不满足精度要求的结果。

解决建议 #

1. 修正变量的可行域 #

原约束中的分母sigma_eigs[i] - x[i]^2不能为0或负数，否则分式无意义或出现数值奇异。因此需要为每个变量添加上限：x[i] < sqrt(sigma_eigs[i])，结合非负约束，Bounds应设置为：

upper_bounds = np.sqrt(sigma_eigs)
bounds = Bounds([0, 0], upper_bounds)

2. 选择合理的初始点 #

原初始点x0 = [1,1]中，x0[1]=1已经大于sqrt(0.24800215)≈0.498，超出可行域，导致求解器初始状态数值不稳定。建议选择可行域内的初始点，例如：

x0 = np.array([0.5, 0.4])

3. 提供约束的雅可比矩阵 #

手动计算约束的梯度（雅可比矩阵），避免求解器使用数值微分带来的误差，提升稳定性：

def non_linear_eq_jac(x):
    # 计算约束函数对x0的导数
    denom0 = sigma_eigs[0] - x[0]**2
    jac0 = (denom0 + 2*x[0]**2) / (denom0 ** 2)
    # 计算约束函数对x1的导数
    denom1 = sigma_eigs[1] - x[1]**2
    jac1 = -(denom1 + 2*x[1]**2) / (denom1 ** 2)
    return np.array([jac0, jac1])

constr1 = NonlinearConstraint(non_linear_eq, 0, 0, jac=non_linear_eq_jac)

4. 调整求解器参数 #

设置更高的精度要求和迭代次数，确保求解收敛到可行解：

options = {
    'xtol': 1e-8,
    'gtol': 1e-8,
    'maxiter': 1000,
    'verbose': 1  # 可选，查看求解过程细节
}
result = minimize(equations, x0, method='trust-constr', constraints=constr1, bounds=bounds, options=options)

修正后的完整代码 #

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

n = 2
c = 2.098024636032762
sigma_eigs = np.array([2.56008479, 0.24800215])

def objective(x):
    return (x[0] + x[1] - c) ** 2

def non_linear_eq(x):
    return x[0]/(sigma_eigs[0]-x[0]**2) - x[1]/(sigma_eigs[1]-x[1]**2)

def non_linear_eq_jac(x):
    denom0 = sigma_eigs[0] - x[0]**2
    jac0 = (denom0 + 2*x[0]**2) / (denom0 ** 2)
    denom1 = sigma_eigs[1] - x[1]**2
    jac1 = -(denom1 + 2*x[1]**2) / (denom1 ** 2)
    return np.array([jac0, jac1])

# 设置可行域
upper_bounds = np.sqrt(sigma_eigs)
bounds = Bounds([0, 0], upper_bounds)

# 可行域内的初始点
x0 = np.array([0.5, 0.4])

# 带雅可比的非线性约束
constr1 = NonlinearConstraint(non_linear_eq, 0, 0, jac=non_linear_eq_jac)

# 求解器参数
options = {
    'xtol': 1e-8,
    'gtol': 1e-8,
    'maxiter': 1000
}

result = minimize(objective, x0, method='trust-constr', constraints=constr1, bounds=bounds, options=options)

x = result.x
print("优化结果x =", x)
print("约束满足情况：", non_linear_eq(x))
print("目标函数值：", objective(x))
print("是否满足非负约束：", np.all(x >= 0))
print("是否满足分母约束：", np.all(x**2 < sigma_eigs))

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