У меня довольно простая ограниченная задача оптимизации, но я получаю разные ответы в зависимости от того, как я это делаю.Давайте сначала разберемся с функцией импорта и красивой печати:
import numpy as np
from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint, SR1
def print_res( res, label ):
print("\n\n ***** ", label, " ***** \n")
print(res.message)
print("obj func value at solution", obj_func(res.x))
print("starting values: ", x0)
print("ending values: ", res.x.astype(int) )
print("% diff", (100.*(res.x-x0)/x0).astype(int) )
print("target achieved?",target,res.x.sum())
Образцы данных очень просты:
n = 5
x0 = np.arange(1,6) * 10_000
target = x0.sum() + 5_000 # increase sum from 15,000 to 20,000
Вот ограниченная оптимизация (включая якобианов).Словом, целевая функция, которую я хочу минимизировать, - это просто сумма квадратов процентных изменений от начальных значений до конечных значений.Линейное ограничение равенство просто требует, чтобы x.sum() равнялась константе.
def obj_func(x):
return ( ( ( x - x0 ) / x0 ) ** 2 ).sum()
def obj_jac(x):
return 2. * ( x - x0 ) / x0 ** 2
def constr_func(x):
return x.sum() - target
def constr_jac(x):
return np.ones(n)
И для сравнения я переформулировал как неограниченную минимизацию с использованием ограничения равенства длязамените x[0] на функцию x[1:].Обратите внимание, что неограниченная функция передается x0[1:], тогда как ограниченная функция передается x0.
def unconstr_func(x):
x_one = target - x.sum()
first_term = ( ( x_one - x0[0] ) / x0[0] ) ** 2
second_term = ( ( ( x - x0[1:] ) / x0[1:] ) ** 2 ).sum()
return first_term + second_term
Затем я пытаюсь минимизировать тремя способами:
- Без ограничений с помощью 'Nelder-Mead '
- Ограничено' trust-constr '(с / без якобиана)
- Ограничено с' SLSQP '(с / без якобиана)
Код:
##### (1) unconstrained
res0 = minimize( unconstr_func, x0[1:], method='Nelder-Mead') # OK, but weird note
res0.x = np.hstack( [target - res0.x.sum(), res0.x] )
print_res( res0, 'unconstrained' )
##### (2a) constrained -- trust-constr w/ jacobian
nonlin_con = NonlinearConstraint( constr_func, 0., 0., constr_jac )
resTCjac = minimize( obj_func, x0, method='trust-constr',
jac='2-point', hess=SR1(), constraints = nonlin_con )
print_res( resTCjac, 'trust-const w/ jacobian' )
##### (2b) constrained -- trust-constr w/o jacobian
nonlin_con = NonlinearConstraint( constr_func, 0., 0. )
resTC = minimize( obj_func, x0, method='trust-constr',
jac='2-point', hess=SR1(), constraints = nonlin_con )
print_res( resTC, 'trust-const w/o jacobian' )
##### (3a) constrained -- SLSQP w/ jacobian
eq_cons = { 'type': 'eq', 'fun' : constr_func, 'jac' : constr_jac }
resSQjac = minimize( obj_func, x0, method='SLSQP',
jac = obj_jac, constraints = eq_cons )
print_res( resSQjac, 'SLSQP w/ jacobian' )
##### (3b) constrained -- SLSQP w/o jacobian
eq_cons = { 'type': 'eq', 'fun' : constr_func }
resSQ = minimize( obj_func, x0, method='SLSQP',
jac = obj_jac, constraints = eq_cons )
print_res( resSQ, 'SLSQP w/o jacobian' )
Вот несколько упрощенных выводов (и, конечно, вы можете запустить код для получения полного вывода):
starting values: [10000 20000 30000 40000 50000]
***** (1) unconstrained *****
Optimization terminated successfully.
obj func value at solution 0.0045454545454545305
ending values: [10090 20363 30818 41454 52272]
***** (2a) trust-const w/ jacobian *****
The maximum number of function evaluations is exceeded.
obj func value at solution 0.014635854609684874
ending values: [10999 21000 31000 41000 51000]
***** (2b) trust-const w/o jacobian *****
`gtol` termination condition is satisfied.
obj func value at solution 0.0045454545462939935
ending values: [10090 20363 30818 41454 52272]
***** (3a) SLSQP w/ jacobian *****
Optimization terminated successfully.
obj func value at solution 0.014636111111111114
ending values: [11000 21000 31000 41000 51000]
***** (3b) SLSQP w/o jacobian *****
Optimization terminated successfully.
obj func value at solution 0.014636111111111114
ending values: [11000 21000 31000 41000 51000]
Примечания:
(1) и (2b) являются правдоподобными решениями в том смысле, что они достигают значительно более низких значений целевой функции, и интуитивно мы ожидаем, что переменные с большими начальными значениями будут двигаться больше (как в абсолютном, так и в процентном отношении)термины), чем меньшие.
Добавление якобиана к «trust-const» приводит к тому, что он получает неправильный ответ (или, по крайней мере, худший ответ), а также превышает максимальные итерации.Может быть, якобиан ошибается, но функция настолько проста, что я почти уверен, что она правильная (?)
«SLSQP» не работает без или безЯкобиан поставил, но работает очень быстро и утверждает, что успешно завершен.Это кажется очень тревожным из-за того, что получение неправильного ответа и утверждение, что он был успешно завершен, является в значительной степени худшим из возможных результатов.
Первоначально я использовал очень маленькие начальные значения и цели (всего 1/1000).из того, что я имею выше), и в этом случае все 5 подходов выше работают нормально и дают одинаковые ответы.Мои выборочные данные все еще очень малы, и кажется, что они обрабатывают 1,2,..,5, но не 1000,2000,..5000.
FWIW. Обратите внимание, что все 3 неверных результата попали в цельдобавляя 1000 к каждому начальному значению - это удовлетворяет ограничению, но далеко не минимизирует целевую функцию (переменные b / c с более высокими начальными значениями должны быть увеличены больше, чем более низкие, чтобы минимизировать сумму квадратов процентных разниц).
Итак, мой вопрос на самом деле заключается в том, что здесь происходит, и почему только (1) и (2b), кажется, работают?
В целом, я хотел бы найтихороший основанный на Python подход к этой и аналогичным задачам оптимизации, и он рассмотрит ответы с использованием других пакетов, кроме scipy, хотя лучший ответ в идеале также будет касаться того, что здесь происходит с scipy (например, это ошибка пользователя или ошибка, которую я должен публиковать на github?).