2

我正在尝试使用 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 以加快操作：

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)

但是，优化总是返回：

     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) 有效），通过这样做优化会成功进行，但我想避免这样做（我必须在某一点反转函数，然后它可以回到原始值会变得混乱）。此外，在我看来，这种关系非常基本，温度和硫之间的值在不同的尺度上并没有那么极端。会是什么呢？谢谢！

4

2 回答 2

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()

输出：

.
.
.
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])

情节：

在此处输入图像描述

于 2019-07-08T15:41:55.170 回答
0

问题在于边界。对于极小的浮点数，某些数字的二进制表示不存在。在内部，优化器正在比较例如 99.9999999 -100 >0 并确定它们不相等（边界不满足），如果您的约束是 XY==.0 。达到 maxEval 后，它得出的结论是没有满足约束的组合。为了避免这种情况，如果解决方案中没有问题，请将边界更改为 (0, 0.000001) 而不是 (0.,0.)

于 2021-06-10T20:33:38.253 回答