Приводится Фортран-реализация поиска минимума функции методом Хука-Дживса.
Замечание. Для поиска максимума нужно умножить на -1 результат, получаемый при вычислении значения оптимизируемой функции.
В качестве тестовой берется функция Розенброка:
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-разрядном процессоре с двойной точностью. Программа может работать с произвольной оптимизируемой функцией.
В этом и следующем алгоритмах F – это оптимизируемая функция, переменная X сначала содержит начальную точку поиска, затем координаты текущей точки, а после завершения цикла – координаты точки-результата.
Функция Erkunden, выполняющая пробный шаг и возвращающая найденную точку, основана на следующем алгоритме:
В главной программе задаются параметры алгоритма и начальная точка поиска экстремума. Оптимизируемая функция объявляется в главной программе с атрибутом External и поэтому может быть использована в качестве параметра процедуры.
program tst
integer, parameter :: n = 2
real(8) X(n), delta, delta_thres, rho
real(8), external :: Rosenbrock
delta = 0.4_8
delta_thres = 1.0e-5_8
rho = 0.5_8
! Начальная точка
X(1) = 5.0_8
X(2) = -5.0_8
! Поиск минимума функции Розенброка
call hjeeves(Rosenbrock, X, n, delta, delta_thres, rho)
!
! Результат
print *, 'X = ', X
end program tst
!
! Выполняет поиск экстремума (минимума) функции F
subroutine hjeeves(F, X, n, delta, delta_thres, rho)
real(8) :: F, X(n), delta, delta_thres, rho
integer n
real(8) :: XNxt(n), XNxt2(n)
integer k, kMx, j
k = 0
! Предельное число итераций (убрать после отладки)
kMx = 100000
print *, 'Initial point = ', X
print *, 'rho = ', rho
print *, 'delta_thres = ', delta_thres
do while (delta > delta_thres .and. k < kMx)
xNxt = Erkunden(X)
! Функция ql сравнивает два вектора на предмет их равенства
if (ql(X, XNxt)) then
delta = rho * delta
else
XNxt2 = 2.0_8 * XNxt - X
XNxt2= Erkunden(XNxt2)
if (F(XNxt2, n) < F(XNxt, n)) then
X = XNxt2
else
X = XNxt
end if
end if
k = k + 1
end do
print *, 'Number of iterations k = ', k
!
contains
!
! Выполняет покоординатное сравнение двух векторов
logical function ql(v, v2)
real(8) v(n), v2(n)
integer j
ql = .true.
do j = 1, n
ql = ql .and. (v(j) == v2(j))
end do
end function ql
!
! Поочередно делает шаг по каждой из координатных осей
! Выбирается шаг, дающий лучший результат
function Erkunden(X)
real(8), dimension(n) :: Erkunden, X
real(8) xNxt(n), fvX, fvXNxt
integer j
fvX = F(X, n)
Erkunden = X
do j = 1, n
XNxt(j) = X(j) + delta
fvXNxt = F(XNxt, n)
if (fvXNxt < fvX) then
Erkunden = XNxt
fvX = fvXNxt
else
XNxt(j) = X(j) - delta
fvXNxt = F(XNxt, n)
if (fvXNxt < fvX) then
Erkunden = XNxt
fvX = fvXNxt
end if
end if
end do
end function Erkunden
!
end subroutine hjeeves
!
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
Приведенная реализация метода Хука-Дживса находит минимум тестовой функции Розенброка, но успех достигается не всегда. Так, если в качестве начальной точки взять, например, точку
(5.0 -5.0),
то минимум будет найден (потребуется 10674 итерации). Если же начать с точки
(-5.0 5.0),
то при заданных в программе параметрах алгоритма вычисления будут закончены после 28799 итераций, но минимум так и не будет достигнут.