Scipy minimize不能正确地解决约束优化(尽管显示成功)

vx6bjr1n  于 2023-08-05  发布在  其他
关注(0)|答案(1)|浏览(230)

下面是一个简单的优化问题。这是一个凸优化,但cvxpy不能识别它,所以我使用SciPy最小化来解决它。问题如下:

min_{x_1, x_2} (x_1 + x_2 - c)^2

subject to: x_1 / (sig - x_1^2) =  x_2 / (sig - x_2^2)

x_1 , x_2 >= 0

字符串
密码是

import numpy as np
import cvxpy as cp
from scipy.linalg import sqrtm
from scipy.optimize import minimize, Bounds, LinearConstraint, 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)

# initial guess for the x values
x0 = np.ones(n) 

# Constraints
bounds = Bounds([0, 0])  # [min x0, min x1]
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")

const_result = [x[i] / (sigma_eigs[i] - x[i]**2) for i in range(n)]


使用'trust-constr'的结果是

[-2.75911573e-08  1.22830289e+07]


解决方案是

x = [-2.75911573e-08  1.22830289e+07]


错误:某些x值为非正残差为:150872747356720.12错误:解不准确。残差很高
优化消息是
消息:xtol终止条件满足。
显然,在这种方法下,变量之一是发散的。
使用其他方法也会给出不满足约束的随机结果。例如,在默认求解器result = minimize(equations, x0 ,constraints=constr1, bounds= bounds)中运行会给出另一个错误答案:[0. 0.]解为x = [0. 0.]所有x值均为正残差为:4.401707373400404错误:解不准确。残差很高
任何帮助都是非常感谢的。

c9qzyr3d

c9qzyr3d1#

我简化了你的代码(见下文),并没有得到你的“正确”的解决方案。仔细观察,当我尝试你的“正确”解决方案时,我看到约束没有得到满足,所以它实际上是不正确的。唯一满足约束的点(在0和100之间)是x=[0,0](我通过蛮力测试了这个点)。

import numpy as np
from scipy.optimize import minimize

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

def equations(x):
    x1, x2 = x
    eq = x1 + x2 - c
    return eq**2

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

result = minimize(equations, np.ones(n),
                  constraints={"type": "eq", "fun": non_linear_eq},
                  bounds=[(0, None), (0, None)])

x = result.x
fun = result.fun

print(f"The solution is {x = }")

if np.all(x >= 0):
    print("All x values are positive.")
else:
    print("Error: Some of the x values are non-positive.")

print(f"The objective function at the solution: {fun}")

constraint_result = non_linear_eq(x)
if np.isclose(constraint_result, 0):
    print("The constraint is satisfied.")
else:
    print(f"Error: The constraint is not satisfied. {constraint_result = }")

x1 = np.linspace(0, 100, 1000)
x2 = np.linspace(0, 10, 1000)
X1, X2 = np.meshgrid(x1, x2)
x_2d = [x1[:, None], x2[None, :]]
F = equations(x_2d)
const_satisfied = np.where(np.isclose(non_linear_eq(x_2d), 0))
print(f"Points that satisfy the constraint: {np.stack(const_satisfied).T}")

字符串
输出量:

The solution is x = array([0., 0.])
All x values are positive.
The objective function at the solution: 4.401707373400404
The constraint is satisfied.
Points that satisfy the constraint: [[0 0]]

相关问题