Найти экстремумы функции python

Тестирование метода Хука-Дживса на функции Розенброка показало его чувствительность к выбору начальной точки поиска экстремума функции.
В этой работе, опять же на функции Розенброка, реализуется и тестируется метод Нелдера-Мида [1].
Приводится Python-реализации поиска минимума функции методом Нелдера-Мида.
При желании можно познакомиться и с Фортран-, и C#-реализациями метода Нелдера-Мида.
Замечание. Для поиска максимума нужно умножить на -1 результат, получаемый при вычислении значения оптимизируемой функции.
В качестве тестовой берется функция Розенброка:

Глобальный минимум функции равен 0.0 и находится в точке (x1, x2) = (1.0, 1.0).
При работе с функцией Розенброка для проверки используемого метода в качестве начальной часто берут точку
(-5, 10)
или точку
(-2.048, 2.048).
Приводимая программа может работать с произвольной оптимизируемой функцией.
Код, реализующий метод Нелдера-Мида, предваряет Python-программа построения графика функции Розенброка.

Вывод графика функции Розенброка

Выводится изометрическая проекция 3d-поверхности. Используется библиотека Matplotlib [2].
График показан на рис. 1.

Рис. 1. График функции Розенброка

Для его вывода употреблен следующий код:

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

# Формирование сетки
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1.5, 3, 0.1)
X, Y = np.meshgrid(X, Y)
# Функция Розенброка
Z = (1.0 — X)**2 + 100.0 * (Y — X * X)**2
#
fig = plt.figure()
# Будем выводить 3d-проекцию графика функции
ax = fig.gca(projection = ‘3d’)
# Вывод поверхности
surf = ax.plot_surface(X, Y, Z, cmap = cm.Spectral, linewidth = 0, antialiased = False)
# Метки осей координат
ax.set_xlabel(‘X’)
ax.set_ylabel(‘Y’)
ax.set_zlabel(‘Z’)
# Настройка оси X
for label in ax.xaxis.get_ticklabels():
label.set_color(‘black’)
label.set_rotation(-45)
label.set_fontsize(9)
# Настройка оси Y
for label in ax.yaxis.get_ticklabels():
label.set_fontsize(9)
# Настройка оси Z
for label in ax.zaxis.get_ticklabels():
label.set_fontsize(9)
# Изометрия
ax.view_init(elev = 30, azim = 45)
# Шкала цветов
fig.colorbar(surf, shrink = 0.5, aspect = 5)
# Отображение результата (рис. 1)
plt.show()

Python-реализация метода Нелдера-Мида

Используется библиотека SciPy [3]. Реализуются два варианта вызова библиотечной процедуры оптимизации:

Задается начальная точка поиска минимума функции.
Задается начальный симплекс поиска минимума функции.

# Вариант 1
import numpy as np
import math # Для sqrt()
import scipy.optimize as opt
# Функция Розенброка
def Rosenbrock(X):
return (1.0 — X[0])**2 + 100.0_8 * (X[1] — X[0] * X[0] )**2
#
n = 2
x0 = np.zeros(2, dtype = float) # Вектор с двумя элементами типа float
# Начальная точка поиска минимума функции
x0[0] = -5.0
x0[1] = 10.0
xtol = 1.0e-5 # Точность поиска экстремума
# Находим минимум функции
res = opt.minimize(Rosenbrock, x0, method = ‘Nelder-Mead’, options = )
print(res)

Результат:
final_simplex: (array([[1.00000132, 1.0000028 ],
[1.00000014, 0.99999997],
[0.99999681, 0.99999355]]), array([4.38559817e-12, 9.00569749e-12, 1.05977059e-11]))
fun: 4.385598172677925e-12
message: ‘Optimization terminated successfully.’
nfev: 269 (число оценок функции)
nit: 143 (число итераций)
status: 0
success: True
x: array([1.00000132, 1.0000028 ])

# Вариант 2
import numpy as np
import math # Для sqrt()
import scipy.optimize as opt
# Функция Розенброка
def Rosenbrock(X):
return (1.0 — X[0])**2 + 100.0_8 * (X[1] — X[0] * X[0] )**2
#
# Процедура формирования начального симплекса
def makeInitialSimplex(X, L, n, initialSimplex):
qn = math.sqrt(1.0 + n) — 1.0
q2 = L / math.sqrt(2.0) * n
r1 = q2 * (qn + n)
r2 = q2 * qn
initialSimplex[0, :] = X
for j in range(n):
initialSimplex[j + 1, :] = X + r2
for i in range(n):
initialSimplex[i + 1, i] += (r1 — r2)
#
n = 2
x0 = np.zeros(2, dtype = float) # Вектор с двумя элементами типа float
# Начальная точка поиска минимума функции
x0[0] = -5.0
x0[1] = 10.0
xtol = 1.0e-5 # Точность поиска экстремума
# Начальная симплекс поиска минимума функции
initialSimplex = np.zeros((n + 1, n), dtype = float)
L = 0.4 # Длина ребра начального симплекса
# Формируем начальный симплекс
makeInitialSimplex(x0, L, n, initialSimplex)
# Находим минимум функции
res = opt.minimize(Rosenbrock, x0, method = ‘Nelder-Mead’, options = )
print(res)

Результат:
final_simplex: (array([[1.00000168, 1.00000305],
[0.99999702, 0.99999376],
[0.99999679, 0.99999416]]), array([1.18711392e-11, 1.70297440e-11, 4.41663303e-11]))
fun: 1.1871139195622139e-11
message: ‘Optimization terminated successfully.’
nfev: 255 (число оценок функции)
nit: 138 (число итераций)
status: 0
success: True
x: array([1.00000168, 1.00000305])

Литература

Nichtrestringierte Optimierung. Numerik III, 21. Januar 2013.
Matplotlib: Научная графика в Python. [Электронный ресурс] URL: https://pythonworld.ru/novosti-mira-python/scientific-graphics-in-python.html.
SciPy Tutorial. Optimization. [Электронный ресурс] URL: https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html.

Источник

Пример использования библиотеки sympy для поиска экстремума функции

найдем значение x, которому соответсвует минимум или максимум функции.

 # numpy всегда пригодится import numpy as np # Для построения графиков import matplotlib.pyplot as plt # Научный питон from scipy import optimize # Символьный питон from sympy import *

Объявляем используемые в функции символьные параметры

a, b, h, k, x, v1 = symbols('a, b, h, k, x, v1')

f = (sqrt(x**2+a**2) + k*sqrt(b**2+(h-x)**2))/v1

Найдем её производную при помощи функции библиотеки sympy diff. Первый аргумент функции diff – дифференцируемое выражение, второй – переменная, по которой необходимо найти производную:

В результате переменная df будет содержать следующее выражение

>> df (k*(-h + x)/sqrt(b**2 + (h - x)**2) + x/sqrt(a**2 + x**2))/v1

В полученном выражении для производной заменим символы (параметры) a, b, h, k, v1 их значениями (a=10, b=10, h=10, k=5, v1=5). Для этого создаем словарь

params = a:10, b:10, h:10.0, k:5, v1: 5>

который подставим в найденную производную, используюя метод subs

Результатом будет выражение, которое зависит только от x:

>> df_par x/(5*sqrt(x**2 + 100)) + (x - 10.0)/sqrt((-x + 10.0)**2 + 100)

Создадим на основе символьного выращения f_par лямбда-функцию от x

df_num = lambda xnum: df_par.subs( x: xnum> )

Численным методом найдем значение x, при котором производная обращается в 0. Для этого используем функцию root модуля scipy.optimize, передав этой функции имя лямбда-функции и начальное приближение \(x_0 = 5\):

sol = optimize.root(df_num, 5.0)

Результатом будет следующее значение \(x\):

2023

2022

2021

2020

2019

Решение нестационарной задачи теплопроводности в MATLAB 10 Dec
Модель физического маятника в Simulink 01 Aug
Движение механической системы с двумя степенями свободы 29 Jul
Модель движения наноспутника 03 Jul
Плоская модель сети 30 Jun
Пример использования библиотеки sympy для поиска экстремума функции 09 Jun
Способ очистки орбит от объектов космического мусора 22 Feb
Движение спускаемого аппарата в атмосфере 05 Jan
Истечение газа из ёмкости постоянного объёма 02 Jan

2018

Источник

Генетический алгоритм на Python для поиска глобальных экстремумов

Недавно в рамках учебной деятельности мне понадобилось использовать старый добрый генетический алгоритм в целях поиска минимума и максимума функции от двух переменных. Однако, к моему удивлению, на просторах интернета не нашлось подобной реализации на питоне, а в википедии в статье по генетическому алгоритму именно этот раздел не освещен.

И поэтому я решил написать свой небольшой пакет на Python с наглядной визуализацией алгоритма, по которой будет удобно настраивать этот алгоритм и искать тонкости выбранной модели.

В этой небольшой статье я хотел бы поделиться процессом, наблюдениями и итогами.

Принцип работы алгоритма

Я не стану рассказывать о глобальном принципе работы генетических алгоритмов, но если вы ещё не слышали о таком, то ознакомиться с ним вы можете в википедии.

На данный момент в пакете реализован только один ГА, который параметризуется входными данными через простую гуишку. Расскажу кратко о выбранных генетических функциях и основных алгоритмических решениях.

Однохромосомная особь несёт в каждом своем гене информацию о соответствующей координате х или у. Популяция определяется множеством особей, однако популяция сегментируется по 4 особи. Данное решение, конечно же, обусловлено попыткой избежать сходимость к локальному оптимуму, так как стоит задача о поиске глобального экстремума. Такое разбиение, как показала практика, во многих случаях не дает доминировать одному генотипу во всей популяции, а, наоборот, придает «эволюции» бóльшую динамику. Для каждой такой части популяции применяется следующий алгоритм:

Селекция происходит схоже с метод ранжирования. Выбирается 3 особи с наилучшим показателями фитнесс-функции (т.е. производится сортировка особей в порядке возрастания/убывания заданной пользователем функции, которая и выступает в роли функции приспособления).
Далее применяется функция скрещивания таким образом, что новое поколение (а точнее новый сегмент популяции в 4 особи) получает 2 пары немутировавших генов от особи с лучшим показанием фитнесс-функции и по паре мутировавших генов от двух других особей. Подробнее о составлении функции мутации будет написано в следующем параграфе.

Первичные тесты и наблюдения

Итак, протестируем данный алгоритм на двух простых примерах:

После тестирования и изучения работы алгоритма методом пристального взгляда и случайного тыка выявилось несколько гипотез-закономерностей:

Погрешность алгоритма прямо пропорциональна количеству особей, но в среднем погрешность исчисляется на сотые, хотя при неудачных параметрах ошибка может достигать и десятых
На данный момент мутация происходит прибавлением к гену рандомного числа из полуинтервала из-за чего особи не могут удержаться в окрестности экстремума и происходит относительно сильное колебание в расстоянии между особями разных поколений. В некоторых случаях такие колебания могут оказаться полезными для выхода из локального экстремума, например, периодической функции, однако экстремумы могут находиться на большом расстоянии друг от друга (расстояние в евклидовом понимании)
Во многих случаях точка экстремума достигается особями за 5-15 поколений, а остальные поколения бесполезно «прыгают» в окрестности этого экстремума
Нулевое поколение заполняется случайными числами только в квадрате
и возможны случаи, когда данный квадрат покрывает локальный экстремум, что нам не подходит

Экстремумы функции g будут находится в точках

Этот пример подтверждает и иллюстрирует все выше сказанные наблюдения.

Апгрейд генетического алгоритма

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

к мутировавшему гену. Такая инвариантность мутации иногда мешает корректной работе алгоритма, но есть эффективный способ исправить этот недочёт.

Введём новый параметр, который назовём «сила мутации» (mutation range) и который будет показывать, в каком полуинтервале мутирует ген. Сделаем так, чтобы данный коэффициент мутации был обратно пропорциональна номеру поколения. Т.е. чем больше номер поколения, тем слабее мутируют гены. Это решение позволяет настраивать стартовую область и улучшать точность вычислений при необходимости.

Как видно из примера, теперь, с каждым поколением популяция всё больше сходится в точку экстремума и вычисляет наиболее точные значения за счёт слабых колебаний.

А что же с проблемой локальных экстремумов? Рассмотрим уже знакомый пример.

Видим, что теперь идея с разделением популяции на части работает так, как было задумано. Без этой сегментации особи ранних поколений могли бы выявить ложный доминирующий генотип в локальном экстремуме, что привело бы к неправильному ответу в поставленной задаче. Также заметно качественное улучшение в точности ответа благодаря введенной зависимости мутации от номера поколения.

Итоги

Резюмирую полученный результат:

Правильный набор параметров позволяет достаточно точно найти глобальный экстремум функции от двух переменных
Разбиение популяции на части позволяет во многих случаях избежать сходимости к локальному оптимуму
Введение силы мутации позволяет также избежать сходимости к локальному оптимуму и в разы увеличивает точность результата
Визуализация может помочь в отладке и улучшении алгоритма, а так же помогает понимать ход и принципы самого генетического алгоритма
повертеть графики, попроверять свои гипотезы и посмотреть код можно тут

В следующих сериях.

Реализация других функций селекции, скрещивания и мутации. Например, интересно посмотреть, как в поставленной задаче будет работать турнирный метод селекции aka естественный отбор
Сравнение по времени и эффективности подобных алгоритмов со стандартными методами оптимизации, например, градиентным спуском
Новые функции пакета (вероятно)

Источник