优化最小化与 SLSQP 不兼容的不等式约束
Optimize minimize Inequality constraints incompatible with SLSQP
我正在尝试使用 Scipy 的优化来 运行 一些曲线拟合。我不使用 polyfit 执行此操作,因为考虑到我的关系,我想确保曲线是单调的。所以假设我在硫和温度之间有以下关系:
sulphur = array([ 71., 82., 50., 113., 153., 177., 394., 1239., 2070., 2662., 3516., 4000., 4954., 6314.])
temperature = array([ 70., 90., 140., 165., 210., 235., 265., 330., 350., 390., 410., 435., 540., 580.])
我想找到一条曲线来拟合这种单调递增的关系。
我找到了一些代码,这就是我所拥有的:我计算多项式并在 objective 函数中使用它,并且我将每个点的一阶导数限制为正数。我还使用 polyfit 值作为 x0 以加快操作速度:
x = sul
y = temperature
initial = list(reversed(np.polyfit(sul, temperature, 3)))
Nfeval = 1
def polynomial(p, x):
return p[0]+p[1]*x+p[2]*x**2+p[3]*x**3
def constraint_1st_der(p, x):
return p[1]+2*p[2]*x+3*p[3]*x**2
def objective(p, x):
return ((polynomial(p, x) - y)**2).sum()
def f(p):
return objective(p, x)
def callback(p):
global Nfeval
print(Nfeval, p, constraint_1st_der(p, x))
Nfeval += 1
cons = {'type' : 'ineq', 'fun' : lambda p : constraint_1st_der(p, x)}
res = optimize.minimize(f, x0=np.array(initial), method='SLSQP', constraints=cons, callback = callback)
但是,始终优化 returns:
fun: 4.0156824919527855e+23
jac: array([0.00000000e+00, 0.00000000e+00, 7.02561542e+17, 3.62183986e+20])
message: 'Inequality constraints incompatible'
nfev: 6
nit: 1
njev: 1
status: 4
success: False
x: array([ -111.35802358, 1508.06894349, -2969.11149743, 2223.26354865])
我尝试规范化(比如 sul_norm = sul / max(sul) 有效),通过这样做优化成功进行,但我想避免这样做(我必须反转在某一点起作用,然后它可能会因返回原始值而变得混乱)。此外,在我看来,这种关系非常基本,温度和硫之间的值在不同的尺度上并没有那么极端。会是什么呢?谢谢!
您在这里遇到了各种限制问题:首先是求解器的选择,它会极大地影响您可以进行的优化类型(约束、有界等)。
最重要的是,在您的情况下,您对参数感兴趣并配置预定义集 (x, y)
,因此您应该以 multi-dimensional 方式处理数据以进行与 (x,y)
相关的计算。但是,据我所知,您使用的约束定义适用于一维输出。因此,正如 SO-question 所建议的那样,使用渐变是一个好主意。
不幸的是,在对你的案例进行测试时,结果对我来说并不令人信服,尽管没有错误 code-runs。稍微调整您的代码后,我设法找到了一个不错的解决方法,但我不确定这是否是最好的解决方法。我的想法是使用 Nelder-Mead Solver
并使用等式约束确保你的导数向量都是正的。代码如下:
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
np.set_printoptions(precision=3)
# init data
sulphur = np.array([ 71., 82., 50., 113., 153., 177., 394., 1239., 2070., 2662., 3516., 4000., 4954., 6314.])
temperature = np.array([ 70., 90., 140., 165., 210., 235., 265., 330., 350., 390., 410., 435., 540., 580.])
# init x and y
x = sulphur
y = temperature
# compute initial guess
initial = list(reversed(np.polyfit(x, y, 3)))
Nfeval = 1
# define functions
polynomial = lambda p, x: p[0] + p[1]*x + p[2]*x**2 + p[3]*x**3
derivative = lambda p, x: p[1] + 2*p[2]*x + 3*p[3]*x**2
def constraint(p):
if (derivative(p, x) > 0).all() : return 0
else : return -1
def callback(p):
global Nfeval
print("Evaluations nbr: %3s | p: %5s |cons: %3s" % (Nfeval,
p,
constraint(p)))
Nfeval += 1
func = lambda p: np.linalg.norm(y - polynomial(p, x))
cons = {'type' : 'eq', 'fun' : constraint}
res = optimize.minimize(func,
x0 = np.array(initial),
method ='Nelder-Mead',
constraints = cons,
callback = callback)
print('----------------------------------------------------------------------------------------')
print(res)
# plot results
f = plt.figure(figsize=(10,4))
ax1 = f.add_subplot(131)
ax2 = f.add_subplot(132)
ax3 = f.add_subplot(133)
ax1.plot(x, y,'ro', label = 'data')
ax1.plot(x, polynomial(res.x, x), label = 'fit using optimization', color="orange")
ax1.legend(loc=0)
ax2.plot(x, derivative(res.x, x), label ='derivative')
ax2.legend(loc=0)
ax3.plot(x, y,'ro', label = 'data')
ax3.plot(x, polynomial(initial, x), label = 'fit using polyfit', color="green")
ax3.legend(loc=0)
plt.show()
输出:
.
.
.
Evaluations nbr: 95 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
Evaluations nbr: 96 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
Evaluations nbr: 97 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
Evaluations nbr: 98 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
----------------------------------------------------------------------------------------
final_simplex: (array([[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.881e-09]]), array([159.565, 159.565, 159.565, 159.565, 159.565]))
fun: 159.5654373399882
message: 'Optimization terminated successfully.'
nfev: 168
nit: 99
status: 0
success: True
x: array([ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09])
地块:
问题出在Bounds上。对于极小的浮点数,某些数字的二进制表示不存在。在内部,优化器正在比较实例 99.9999999 -100 >0 并确定它们不相等(绑定不满足),如果您的约束是 X-Y==.0 。在达到 maxEval 后,它得出结论,没有满足约束的组合。为避免这种情况,如果解决方案中没有问题,请将边界更改为 (0, 0.000001) 而不是 (0.,0.)
我正在尝试使用 Scipy 的优化来 运行 一些曲线拟合。我不使用 polyfit 执行此操作,因为考虑到我的关系,我想确保曲线是单调的。所以假设我在硫和温度之间有以下关系:
sulphur = array([ 71., 82., 50., 113., 153., 177., 394., 1239., 2070., 2662., 3516., 4000., 4954., 6314.])
temperature = array([ 70., 90., 140., 165., 210., 235., 265., 330., 350., 390., 410., 435., 540., 580.])
我想找到一条曲线来拟合这种单调递增的关系。
我找到了一些代码,这就是我所拥有的:我计算多项式并在 objective 函数中使用它,并且我将每个点的一阶导数限制为正数。我还使用 polyfit 值作为 x0 以加快操作速度:
x = sul
y = temperature
initial = list(reversed(np.polyfit(sul, temperature, 3)))
Nfeval = 1
def polynomial(p, x):
return p[0]+p[1]*x+p[2]*x**2+p[3]*x**3
def constraint_1st_der(p, x):
return p[1]+2*p[2]*x+3*p[3]*x**2
def objective(p, x):
return ((polynomial(p, x) - y)**2).sum()
def f(p):
return objective(p, x)
def callback(p):
global Nfeval
print(Nfeval, p, constraint_1st_der(p, x))
Nfeval += 1
cons = {'type' : 'ineq', 'fun' : lambda p : constraint_1st_der(p, x)}
res = optimize.minimize(f, x0=np.array(initial), method='SLSQP', constraints=cons, callback = callback)
但是,始终优化 returns:
fun: 4.0156824919527855e+23
jac: array([0.00000000e+00, 0.00000000e+00, 7.02561542e+17, 3.62183986e+20])
message: 'Inequality constraints incompatible'
nfev: 6
nit: 1
njev: 1
status: 4
success: False
x: array([ -111.35802358, 1508.06894349, -2969.11149743, 2223.26354865])
我尝试规范化(比如 sul_norm = sul / max(sul) 有效),通过这样做优化成功进行,但我想避免这样做(我必须反转在某一点起作用,然后它可能会因返回原始值而变得混乱)。此外,在我看来,这种关系非常基本,温度和硫之间的值在不同的尺度上并没有那么极端。会是什么呢?谢谢!
您在这里遇到了各种限制问题:首先是求解器的选择,它会极大地影响您可以进行的优化类型(约束、有界等)。
最重要的是,在您的情况下,您对参数感兴趣并配置预定义集 (x, y)
,因此您应该以 multi-dimensional 方式处理数据以进行与 (x,y)
相关的计算。但是,据我所知,您使用的约束定义适用于一维输出。因此,正如 SO-question 所建议的那样,使用渐变是一个好主意。
不幸的是,在对你的案例进行测试时,结果对我来说并不令人信服,尽管没有错误 code-runs。稍微调整您的代码后,我设法找到了一个不错的解决方法,但我不确定这是否是最好的解决方法。我的想法是使用 Nelder-Mead Solver
并使用等式约束确保你的导数向量都是正的。代码如下:
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
np.set_printoptions(precision=3)
# init data
sulphur = np.array([ 71., 82., 50., 113., 153., 177., 394., 1239., 2070., 2662., 3516., 4000., 4954., 6314.])
temperature = np.array([ 70., 90., 140., 165., 210., 235., 265., 330., 350., 390., 410., 435., 540., 580.])
# init x and y
x = sulphur
y = temperature
# compute initial guess
initial = list(reversed(np.polyfit(x, y, 3)))
Nfeval = 1
# define functions
polynomial = lambda p, x: p[0] + p[1]*x + p[2]*x**2 + p[3]*x**3
derivative = lambda p, x: p[1] + 2*p[2]*x + 3*p[3]*x**2
def constraint(p):
if (derivative(p, x) > 0).all() : return 0
else : return -1
def callback(p):
global Nfeval
print("Evaluations nbr: %3s | p: %5s |cons: %3s" % (Nfeval,
p,
constraint(p)))
Nfeval += 1
func = lambda p: np.linalg.norm(y - polynomial(p, x))
cons = {'type' : 'eq', 'fun' : constraint}
res = optimize.minimize(func,
x0 = np.array(initial),
method ='Nelder-Mead',
constraints = cons,
callback = callback)
print('----------------------------------------------------------------------------------------')
print(res)
# plot results
f = plt.figure(figsize=(10,4))
ax1 = f.add_subplot(131)
ax2 = f.add_subplot(132)
ax3 = f.add_subplot(133)
ax1.plot(x, y,'ro', label = 'data')
ax1.plot(x, polynomial(res.x, x), label = 'fit using optimization', color="orange")
ax1.legend(loc=0)
ax2.plot(x, derivative(res.x, x), label ='derivative')
ax2.legend(loc=0)
ax3.plot(x, y,'ro', label = 'data')
ax3.plot(x, polynomial(initial, x), label = 'fit using polyfit', color="green")
ax3.legend(loc=0)
plt.show()
输出:
.
.
.
Evaluations nbr: 95 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
Evaluations nbr: 96 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
Evaluations nbr: 97 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
Evaluations nbr: 98 | p: [ 1.400e+02 1.830e-01 -4.203e-05 3.882e-09] |cons: 0
----------------------------------------------------------------------------------------
final_simplex: (array([[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09],
[ 1.400e+02, 1.830e-01, -4.203e-05, 3.881e-09]]), array([159.565, 159.565, 159.565, 159.565, 159.565]))
fun: 159.5654373399882
message: 'Optimization terminated successfully.'
nfev: 168
nit: 99
status: 0
success: True
x: array([ 1.400e+02, 1.830e-01, -4.203e-05, 3.882e-09])
地块:
问题出在Bounds上。对于极小的浮点数,某些数字的二进制表示不存在。在内部,优化器正在比较实例 99.9999999 -100 >0 并确定它们不相等(绑定不满足),如果您的约束是 X-Y==.0 。在达到 maxEval 后,它得出结论,没有满足约束的组合。为避免这种情况,如果解决方案中没有问题,请将边界更改为 (0, 0.000001) 而不是 (0.,0.)