Оптимизация с помощью Python (scipy.optimize)

Question

Я пытаюсь максимизировать следующую функцию, используя scipy.optimize от Python. Однако после многих попыток он не работает. Функция и мой код вставляются ниже. Спасибо за помощь!

Проблема

Maximize [sum (x_i/y_i)**gamma]**(1/gamma) 
subject to the constraint sum x_i = 1; x_i is in the interval (0,1). 

x является вектором выбора переменных; y - вектор параметров; gamma - параметр. Сумма x s должна быть равна сумме. И каждый x должен находиться в интервале (0,1).

Код

def objective_function(x, y): 
    sum_contributions = 0 
    gamma = 0.2 

    for count in xrange(len(x)): 
     sum_contributions += (x[count]/y[count]) ** gamma 
    value = math.pow(sum_contributions, 1/gamma) 
    return -value 

cons = ({'type': 'eq', 'fun': lambda x: np.array([sum(x) - 1])}) 

y = [0.5, 0.3, 0.2] 
initial_x = [0.2, 0.3, 0.5] 

opt = minimize(objective_function, initial_x, args=(y,), method='SLSQP', 
constraints=cons,bounds=[(0, 1)] * len(x)) 

Answer 1

Иногда числовая оптимизатор не работает по какой-либо причине. Мы можем параметризовать проблему несколько иначе, и она будет работать. (И может работать быстрее)

К примеру, за пределы (0,1), мы можем иметь преобразования функции, такие, что значения в (-inf, +inf), после того, как трансформируются, будет в конечном итоге в (0,1)

Мы можем сделать подобный трюк с ограничения равенства. Например, мы можем уменьшить размер от 3 до 2, так как последний элемент в x должен быть 1-sum(x).

Если он по-прежнему не работает, мы можем переключиться на оптимизатор, который не требует информации от производной, такой как Nelder Mead.

А также есть Lagrange multiplier.

In [111]: 

def trans_x(x): 
    x1 = x**2/(1+x**2) 
    z = np.hstack((x1, 1-sum(x1))) 
    return z 

def F(x, y, gamma = 0.2): 
    z = trans_x(x) 
    return -(((z/y)**gamma).sum())**(1./gamma) 
In [112]: 

opt = minimize(F, np.array([0., 1.]), args=(np.array(y),), 
       method='Nelder-Mead') 
opt 
Out[112]: 
    status: 0 
    nfev: 96 
success: True 
    fun: -265.27701747828007 
     x: array([ 0.6463264, 0.7094782]) 
message: 'Optimization terminated successfully.' 
    nit: 52 

Результат:

In [113]: 

trans_x(opt.x) 
Out[113]: 
array([ 0.29465097, 0.33482303, 0.37052601]) 

И мы можем визуализировать его с:

In [114]: 

x1 = np.linspace(0,1) 
y1 = np.linspace(0,1) 
X,Y = np.meshgrid(x1,y1) 
Z = np.array([F(item, y) for item 
       in np.vstack((X.ravel(), Y.ravel())).T]).reshape((len(x1), -1), order='F') 
Z = np.fliplr(Z) 
Z = np.flipud(Z) 
plt.contourf(X, Y, Z, 50) 
plt.colorbar() 

Answer 2

Даже трудные эти вопросы немного устарели, я хотел бы добавить альтернативное решение, которое может быть полезно для других, которые наткнулись на этот вопрос в будущем ,

Оказывается, наша проблема разрешима аналитически. Вы можете начать записывать лагранжиан (равенство стесненной) задача оптимизации:

L = \sum_i (x_i/y_i)^\gamma - \lambda (\sum x_i - 1) 

Оптимальное решение найдено, установив первую производную от этого лагранжиана к нулю:

0 = \partial L/\partial x_i = \gamma x_i^{\gamma-1}/\y_i - \lambda 
=> x_i \propto y_i^{\gamma/(\gamma - 1)} 

Используя это представление проблема оптимизации может быть решена просто и эффективно:

In [4]: 
def analytical(y, gamma=0.2): 
    x = y**(gamma/(gamma-1.0)) 
    x /= np.sum(x) 
    return x 
xanalytical = analytical(y) 
xanalytical, objective_function(xanalytical, y) 
Out [4]: 
(array([ 0.29466774, 0.33480719, 0.37052507]), -265.27701765929692) 

решение КТ Чжа элегантно, но это может нарушить позитивности ограничения на третьем ворковать rdinate.Для gamma = 0.2 это не кажется проблемой на практике, но и для различных гамм вы легко не столкнулись с проблемами:

In [5]: 
y = [0.2, 0.1, 0.8] 
opt = minimize(F, np.array([0., 1.]), args=(np.array(y), 2.0), 
       method='Nelder-Mead') 
trans_x(opt.x), opt.fun 
Out [5]: 
(array([ 1., 1., -1.]), -11.249999999999998) 

Для других задач оптимизации с той же вероятностью симплексных ограничений как ваша проблема, но для которых не существует аналитическое решение, возможно, стоит рассмотреть проецируемые градиентные методы или аналогичные. Эти методы используют тот факт, что существует быстрый алгоритм проектирования произвольной точки на этом множестве, см. https://en.wikipedia.org/wiki/Simplex#Projection_onto_the_standard_simplex.

(Чтобы увидеть полный код и лучшее представление уравнений посмотрим на Jupyter ноутбук http://nbviewer.jupyter.org/github/andim/pysnippets/blob/master/optimization-simplex-constraints.ipynb)