Список работ

Поиск минимума функции методом Нелдера-Мида

Содержание

Введение

Приводится Фортран-реализация поиска минимума функции методом Нелдера-Мида. В качестве тестовой берется функция Розенброка

f(x1, x2) = (1 – x1)2 + 100 * (x2 – x12)2

Используется следующий код функции:

real(8) function Rosenbrock(X, n)
 real(8) X(n)
 Rosenbrock = (1.0_8 - X(1))**2 + 100.0_8 * (X(2) - X(1) * X(1) )**2
end function Rosenbrock

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

Симплекс и операции с симплексом в методе Нелдера-Мида

В начале работы алгоритма Нелдера-Мида строится регулярный симплекс: в пространстве Rn - это правильный многогранник, образованный n + 1 равноотстоящими друг от друга вершинами. Так, в R2 – это равносторонний треугольник. В пространстве Rn координаты вершин регулярного симплекса, в котором первая вершина
X = (x1, x2, ..., xn),
можно задать следующей матрицей:

Матрица с координатами вершин регулярного симплекса

в которой в столбце i находятся координаты i-й вершины симплекса и

Формулы для формирования матрицы с координатами вершин регулярного симплекса

где l - длина ребра симплекса.
Так, в R2 из вершины X = (0, 0) регулярный симплекс (рис. 1) задается следующей матрицей:

R2-регулярный симплекс.
Изображение R2-регулярного симплекса.

Рис. 1. R2-регулярный симплекс из вершины X = (0, 0)

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

! Создает из точки X регулярный симплекс с длиной ребра L и с n + 1 вершиной
! Формирует массив FN значений оптимизируемой функции F в вершинах симплекса
subroutine makeSimplex(X, L, n)
 real(8) X(:), L
 integer n
 real(8) qn, q2, r1, r2
 qn = sqrt(1.0_8 + n) - 1.0_8
 q2 = L / sqrt(2.0_8) * n
 r1 = q2 * (qn + n)
 r2 = q2 * qn
 simplex(:, 1) = X
 do i = 2, n + 1
  simplex(:, i) = X + r2
 end do
 do i = 2, n + 1
  simplex(i - 1, i) = simplex(i - 1, i) - r2 + r1
 end do
 do i = 1, n + 1
  ! Значения функции в вершинах симплекса
  FN(i) = F(simplex(:, i), n)
 end do
end subroutine makeSimplex

В алгоритме Нелдера-Мида используются следующие симплекс-операции:

Показано отражение R2-симплекса.

Рис. 2. Отражение R2-симплекса

Показана редукция R2-симплекса.

Рис. 3. Редукция R2-симплекса

Показано сжатие R2-симплекса.

Рис. 4. Сжатие R2-симплекса

Показано растяжение R2-симплекса.

Рис. 5. Растяжение R2-симплекса

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

Формула вычисления центра тяжести симплекса.

где k - номер отражаемой вершины, r - номер текущего шага алгоритма.
Координаты отраженной вершины вычисляются по следующей формуле:

Отражение симплекса.

Обычно cR = 1.0. Координаты прочих вершин симплекса при выполнении операции отражения не меняются.
Центр тяжести симплекса возвращается следующей функцией:

function center_of_gravity(simplex, k, n)
 real(8) center_of_gravity(n)
 real(8), intent(in) :: simplex(n, n + 1)
 integer, intent(in) :: k, n
 integer j
 do j = 1, n
  center_of_gravity(j) = sum(simplex(j, :))
 end do
 center_of_gravity = (center_of_gravity - simplex(:, k)) / n
end function center_of_gravity

Отражение вершины с номером k относительно центра тяжести симплекса обеспечивает следующая подпрограмма:

subroutine reflection(simplex, k, cR, n)
 real(8) simplex(n, n + 1)
 ! cR – коэффициент отражения
 real(8), intent(in) :: cR
 integer, intent(in) :: k, n
 simplex(:, k) = (1 +cR) * center_of_gravity(simplex, k, n) - simplex(:, k)
end subroutine reflection

Редукция симплекса, – уменьшение длины всех ребер симплекса, выполняется к выбранной в алгоритме вершине Xk симплекса. Новые координаты редуцированного симплекса определяются по следующей формуле:

Редукция симплекса.

где γ - коэффициент редукции (положительное число, меньшее 1.0; часто берется значение 0.5), r - номер текущего шага алгоритма.
За редукцию симплекса отвечает следующая подпрограмма:

subroutine reduction(simplex, k, gamma, n)
 real(8) simplex(n, n + 1)
 real(8), intent(in) :: gamma
 integer, intent(in) :: k, n
 real(8) xk(n)
 integer j
 xk = simplex(:, k)
 do j = 1, n
  simplex(:, j) = xk + gamma * (simplex(:, j) - xk)
 end do
 simplex(:, k) = xk
end subroutine reduction

Сжатие симплекса выполняется в направлении Xk - Xc, где Xk - выбранная в алгоритме вершина симплекса, а Xc - центр его тяжести. В результате сжатия изменяются координаты вершины Xk:

Сжатие симплекса.

где β - коэффициент сжатия (положительное число, меньшее 1.0; нередко берется значение 0.5), r - номер текущего шага алгоритма.
Координаты прочих вершин симплекса не изменяются.
Растяжение симплекса выполняется в направлении Xk - Xc, где Xk - выбранная в алгоритме вершина, а Xc - центр его тяжести. В результате растяжения изменяются координаты вершины Xk:

Растяжение симплекса.

где α - коэффициент растяжения (положительное число, большее 1.0; можно употребить, например, значение 2.0), r - номер текущего шага алгоритма.
Координаты прочих вершин симплекса не изменяются.
Сжатие или растяжение симплекса выполняет следующая подпрограмма:

subroutine shrinking_expansion(simplex, k, alpha_beta, n)
 real(8) simplex(n, n + 1)
 real(8), intent(in) :: alpha_beta
 integer, intent(in) :: k, n
 real(8) xc(n)
 xc = center_of_gravity(simplex, k, n)
 simplex(:, k) = xc + alpha_beta * (simplex(:, k) - xc)
end subroutine shrinking_expansion

В случае сильно овражных функций может происходить вырождение (сплющивание) симплекса, что приводит к необходимости восстановления симплекса, которое заключается в следующем:

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

! Восстанавливает симплекс
subroutine simplexRestore()
 real(8) fmi
 real(8) X2(n)
 integer imi(1), imi2(1)
 fmi = minval(fn)
 imi = minloc(fn)
 imi2 = minloc(fn, fn /= fmi)
 X2 = simplex(:, imi(1)) - simplex(:, imi2(1))
 call makeSimplex(simplex(:, imi(1)), sqrt(dot_product(X2, X2)), n)
end subroutine simplexRestore

Алгоритм метода Нелдера-Мида

  1. Начало.
  2. Задать начальную точку X = (x1, x2, …, xn) и n – размер вектора X.
  3. Задать L и L_thres соответственно начальное и наименьшее значения длины ребра симплекса.
  4. Задать cR, alpha, beta и gamma соответственно коэффициенты отражения, растяжения, сжатия и редукции симплекса.
  5. Простроить из точки X регулярный симплекса с длиной ребра L и числом вершин n + 1.
  6. Вычислить значения функции в вершинах начального симплекса.
  7. Пока в симплексе есть ребро, длина которого больше L_thres
      Найти в симплексе вершину k, в которой значение функции максимально (Fma).
      Найти Fmi – минимальное значение функции в текущем симплексе.
      Выполнить отражение вершины k и найти F_R – значение функции в отраженной вершине.
      Если F_R > Fma Тогда
        Выполнить, перемещая вершину k, сжатие симплекса и найти F_S – значение функции в вершине k после сжатия.
       Если F_R > Fma Тогда
        Вернуться к симплексу до отражения и выполнить его редукцию, сохраняя положение вершины k.
        Вычислить значения функции в вершинах редуцированного симплекса.
       Иначе
        Продолжить работу со сжатым симплексом.
       Конец Если
      Иначе Если F_R < Fmi Тогда
       Выполнить растяжение симплекса, перемещая вершину k.
       Вычислить F_E – значение функции в вершине k после растяжения.
       Если F_E > Fmi Тогда
        Отказаться от растяжения и продолжить работу с симплексом, полученным после отражения.
       Иначе
        Продолжить работу с симплексом, полученным после растяжения.
       Конец Если
      Конец Если
     Конец Цикла
  8. Напечатать в качестве результата координаты вершины симплекса с минимальным значением функции.
  9. Останов.

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

В главной программе задаются параметры алгоритма и начальный симплекс (регулярный). Оптимизируемая функция объявляется в главной программе с атрибутом External и поэтому может быть использована в качестве параметра процедуры.

program tst
 integer, parameter :: n = 2
 real(8) X(n), L, L_thres, cR, alpha, beta, gamma
 real(8), external :: Rosenbrock
 L = 0.4_8 ! Начальная длина ребра симплекса
 L_thres = 1.0e-5_8 ! Предельное значение длины ребра симплекса
 cR = 1.0_8 ! Коэффициент отражения симплекса
 alpha = 2.0_8 ! Коэффициент растяжения симплекса
 beta = 0.5_8 ! Коэффициент сжатия симплекса
 gamma = 0.5_8 ! Коэффициент редукции симплекса
 ! Первая вершина начального симплекса (начальная точка)
 X(1) = 5.0_8
 X(2) = -5.0_8
 print *, 'Start point = ', X
 !
 ! Поиск минимума функции Розенброка
 call nelMead(Rosenbrock, X, n, L, L_thres, cR, alpha, beta, gamma)
 !
 ! Результат
 print *, 'Result = ', X
end program tst
!
! Выполняет поиск экстремума (минимума) функции F
subroutine nelMead(F, X, n, L, L_thres, cR, alpha, beta, gamma)
 real(8) F, X(n), L, L_thres, cR, alpha, beta, gamma
 integer n, j, jMx, i
 real(8) X_R(n), simplex(n, n + 1), FN(n + 1)
 real(8) Fmi, Fma, F_R, F_S, F_E
 integer ima(1), k, kr
 ! kr_todo - число шагов алгоритма, после выполнения которых симплекс восстанавливается
 integer, parameter :: kr_todo = 60
 j = 0
 ! Предельное число шагов алгоритма (убрать после отладки)
 jMx = 10000
 print *, 'L = ', L
 print *, 'L_thres = ', L_thres
 print *, 'cR = ', cR
 print *, 'alpha = ', alpha
 print *, 'beta = ', beta
 print *, 'gamma = ', gamma
 call makeSimplex(X, L, n)
 kr = 0
 k = 1
 do while (notStopYet() .and. j < jMx)
  j = j + 1
  kr = kr + 1
  if (kr == kr_todo) then
   kr = 0
   ! Восстановление симплекса
   call simplexRestore()
  end if
  Fmi = minval(FN)
  Fma = maxval(FN)
  ima = maxloc(FN)
  ! Номер отражаемой вершины
   k = ima(1)
   X = simplex(:, k)
  ! Отражение
  call reflection(simplex, k, cR, n)
  X_R = simplex(:, k)
  ! Значение функции в вершине k симплекса после отражения
  F_R = F(simplex(:, k), n)
  if (F_R > Fma) then
  ! Сжатие
   call shrinking_expansion(simplex, k, beta, n)
   ! Значение функции в вершине k симплекса после его сжатия
   F_S = F(simplex(:, k), n)
   if (F_S > Fma) then
    simplex(:, k) = X
  ! Редукция
    call reduction(simplex, k, gamma, n)
    do i = 1, n + 1
     if (i == k) cycle
     ! Значения функций в вершинах симплекса после редукции
      ! В вершине k значение функции сохраняется
     FN(i) = F(simplex(:, i), n)
    end do
   else
    FN(k) = F_S
   end if
  else if (F_R < Fmi) then
   ! Растяжение
   call shrinking_expansion(simplex, k, alpha, n)
   ! Значение функции в вершине k симплекса после его растяжения
   F_E = F(simplex(:, k), n)
   if (F_E > Fmi) then
    simplex(:, k) = X_R
    FN(k) = F_R
   else
    FN(k) = F_E
    end if
   else
    FN(k) = F_R
  end if
  end do
  ! Результат
  print *, 'Number of iterations j = ', j
contains
 !
 ! Возвращает .true., если длина хотя бы одного ребра симплекса превышает L_thres,
 ! или .false. - в противном случае
 logical function notStopYet()
  integer i, j
  real(8) X(n), X2(n)
  notStopYet = .false.
  do i = 1, n
   X = simplex(:, i)
   do j = i + 1, n + 1
     X2 = X - simplex(:, j)
     if (sqrt(dot_product(X2, X2)) > L_thres) then
      notStopYet = .true.
      return
     end if
    end do
  end do
 end function notStopYet
 !
 ! Второй вариант определения точки останова алгоритма
 ! Возвращает .true., если хотя бы одна разница значений функций в вершинах симплекса превышает L_thres,
 ! или .false. - в противном случае
 logical function notStopYet2()
  integer i, j
  real(8) fv
  notStopYet2 = .false.
  do i = 1, n
   fv = FN(i)
   do j = i + 1, n + 1
     if (abs(fv - FN(j)) > L_thres) then
      notStopYet2 = .true.
      return
     end if
    end do
  end do
 end function notStopYet2
 !
 ! Восстанавливает симплекс
 subroutine simplexRestore()
  real(8) fmi
  real(8) X2(n)
  integer imi(1), imi2(1)
  fmi = minval(fn)
  imi = minloc(fn)
  imi2 = minloc(fn, fn /= fmi)
  X2 = simplex(:, imi(1)) - simplex(:, imi2(1))
  call makeSimplex(simplex(:, imi(1)), sqrt(dot_product(X2, X2)), n)
 end subroutine simplexRestore
 !
 ! Создает из точки X регулярный симплекс с длиной ребра L и с n + 1 вершиной
 subroutine makeSimplex(X, L, n)
  real(8) X(:), L
  integer n
  real(8) qn, q2, r1, r2
  qn = sqrt(1.0_8 + n) - 1.0_8
  q2 = L / sqrt(2.0_8) * n
  r1 = q2 * (qn + n)
  r2 = q2 * qn
  simplex(:, 1) = X
  do i = 2, n + 1
   simplex(:, i) = X + r2
  end do
  do i = 2, n + 1
   simplex(i - 1, i) = simplex(i - 1, i) - r2 + r1
  end do
  do i = 1, n + 1
   ! Значения функции в вершинах созданного симплекса
   FN(i) = F(simplex(:, i), n)
  end do
 end subroutine makeSimplex
 !
 ! Вычисляет центр тяжести симплекса
 function center_of_gravity(simplex, k, n)
  real(8) center_of_gravity(n)
  real(8), intent(in) :: simplex(n, n + 1)
  integer, intent(in) :: k, n
  integer j
  do j = 1, n
   center_of_gravity(j) = sum(simplex(j, :))
  end do
  center_of_gravity = (center_of_gravity - simplex(:, k)) / n
 end function center_of_gravity
 !
 ! Выполняет операцию отражения
 subroutine reflection(simplex, k, cR, n)
  real(8) simplex(n, n + 1)
  ! cR - коэффициент отражения
  real(8), intent(in) :: cR
  integer, intent(in) :: k, n
  simplex(:, k) = (1 +cR) * center_of_gravity(simplex, k, n) - simplex(:, k)
 end subroutine reflection
 !
 ! Обеспечивает редукцию симплекса
 subroutine reduction(simplex, k, gamma, n)
  real(8) simplex(n, n + 1)
  real(8), intent(in) :: gamma
  integer, intent(in) :: k, n
  real(8) xk(n)
  integer j
  xk = simplex(:, k)
  do j = 1, n
   simplex(:, j) = xk + gamma * (simplex(:, j) - xk)
  end do
  simplex(:, k) = xk
 end subroutine reduction
 !
 ! Растяжение или сжатие симплекса
 subroutine shrinking_expansion(simplex, k, alpha_beta, n)
  real(8) simplex(n, n + 1)
  real(8), intent(in) :: alpha_beta
  integer, intent(in) :: k, n
  real(8) xc(n)
  xc = center_of_gravity(simplex, k, n)
  simplex(:, k) = xc + alpha_beta * (simplex(:, k) - xc)
 end subroutine shrinking_expansion
end subroutine nelMead
!
real(8) function Rosenbrock(X, n)
 real(8) X(n)
 integer n
 Rosenbrock = (1.0_8 - X(1))**2 + 100.0_8 * (X(2) - X(1) * X(1) )**2
end function Rosenbrock

Второй вариант определения точки останова алгоритма реализуется функцией notStopYet2, которая возвращает .true., если хотя бы одна разница значений функций в различных вершинах симплекса превышает L_thres. В противном случае функция вернет .false.

Заключение

При испытании программы поиск минимума функции Розенброка начинался с различных точек:

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

Литература

  1. Nichtrestringierte Optimierung. Numerik III, 21. Januar 2013.
  2. О. В. Бартеньев. Современный Фортран. - М.: ДИАЛОГ-МИФИ, 2000. - 613 с.

Список работ

Рейтинг@Mail.ru