Оптимизировать минимизацию ограничений неравенства, несовместимых с SLSQP - PullRequest
Новая
       31

Оптимизировать минимизацию ограничений неравенства, несовместимых с SLSQP

1 голос
/ 08 июля 2019

Я пытаюсь выполнить некоторую подгонку кривой с помощью Оптимизации Сципи. Я не делаю это с 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) работает), и при этом оптимизация проходит успешно, но я бы хотел этого избежать (мне нужно инвертировать функцию в одной точке, и тогда это может стать грязным с возвращением к первоначальному значению). Кроме того, мне кажется, что отношения довольно простые, и значения между температурой и серой не настолько сильно различаются. Что бы это могло быть? Спасибо!

1 Ответ

1 голос
/ 08 июля 2019

Здесь у вас есть различные проблемы ограничения: во-первых, выбор решателя, который сильно влияет на то, какую оптимизацию вы можете выполнять (ограниченную, ограниченную и т. Д.).Кроме того, в вашем случае вы заинтересованы в параметрах и располагаете предопределенным набором (x, y), поэтому вам следует обрабатывать данные многомерным способом для вычислений, связанных с (x,y).Тем не менее, определение ограничений, которые вы использовали, насколько мне известно, подходит для одномерных выходов.Поэтому, как предполагает SO-вопрос , использование градиента является хорошей идеей.

К сожалению, при тестировании этого по вашему случаю результаты не были убедительными для меня, несмотря на безошибочный запуск кода.Немного подправив ваш код, мне удалось найти достойный обходной путь, но я не уверен, что он лучший.Моя идея состоит в том, чтобы использовать 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])

Сюжеты:

enter image description here

...