优化最小化与SLSQP不相容的不等式约束

Question

我试着用Scipy的Optimize运行一些曲线拟合。我不使用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.])

我想找到一条曲线来拟合这种单调递增的关系。在

我已经找到了一些代码，这就是我所拥有的：我计算多项式并将其用于目标函数，并且我约束一阶导数对于每个点都是正的。为了加快操作速度，我还使用polyfit值作为x0：

^{pr2}$

但是，优化总是返回：

     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）起作用），这样优化就成功地进行了，但是我想避免这样做（我必须在某一点上反转函数，然后回到原始值时会变得混乱）。而且，在我看来，这种关系是非常基本的，温度和硫之间的值在不同的尺度上并不是那么极端。可能是什么？谢谢！在

Answer 1

这里有各种各样的限制问题：首先是求解器的选择，这会极大地影响到你可以做什么样的优化（约束、有界等）。 除此之外，在您的例子中，您对参数感兴趣，并处理一个预定义的集(x, y)，因此您应该以多维的方式处理数据，以便进行与(x,y)相关的计算。但是，据我所知，您使用的约束定义适用于一维输出。因此，正如SO-question所示，使用梯度是个好主意。在

不幸的是，当在您的案例中进行测试时，结果并不能让我信服，尽管运行了没有错误的代码。在修改了你的代码之后，我设法找到了一个不错的解决方法，但我不确定是否是最好的。我的想法是使用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()

输出：

^{pr2}$

绘图：