Автор cайта:
Владимир
Потемкин

fortran-90@yandex.ru

Язык программирования Фортран

Вычисление минимума функции

FMIN

Программа Fmin вычисляет минимум функции f(x) в заданном интервале [a,b] комбинацией методов золотого сечения и параболической интерполяции.
Подробное описание программы "FMIN" можно найти в книге [Б3].

 

Вызов программы

Fmin(Fun, AX, BX, Eps)

Параметры программы

Fun, AX, BX, Eps - входные параметры;
Real Fmin - возвращаемое значение;
Real Fun(Real X) - исследуемая функция;
Real AX, BX - границы интервала, AX<BX;
Real Eps - погрешность решения.

 

Пример

!Поиск минимума функции одной переменной
program TestFMIN
    use NML
    implicit none
    integer:: Count=0
    real:: A, B, R, FR, Eps
    data A/-1.0/, B/1.0/, Eps/1.0E-5/
    !begin
      R=FMIN(Fun, A, B, Eps)
      FR=Fun(R); Count=Count-1
      print 10, R, Count
      print 11, FR
10    format(/'Minimum =',F11.7,'   Count =',I3)
11    format('The value of the function'/'at the point of minimum ',F12.7)
    !end

    contains

    real function Fun(x)
    real, intent(in):: x
      Fun=x*x*x*x+5.0
      Count=Count+1
      return
    end function Fun

end program TestFMIN

Minimum =  0.0131775   Count = 24
The value of the function
at the point of minimum    5.0000000

 

Вернуться к оглавлению    Скачать FMIN

FM25

Программа FM25 вычисляет минимум функции нескольких переменных f(x1,x2,...,xn) методом Хука-Дживса.
В методе Хука-Дживса используются только значения самой функции, без использований частных производных (градиента) функции.

Для начала вычислений нужно указать начальное приближение (базисную точку) b1 с координатами b0, b1,…, bn и начальный шаг hi (i=1,…,n) по каждой координате. Когда начальные параметры заданы, вычисляем значение функции в базисной точке b1 и в точке, смещённой по оси x1 на шаг h1. Если значение функции в смещённой точке меньше, чем в базисной, то принимаем смещённую точку за новую базисную. В противном случае (функция не стала меньше) делаем смещение в противоположную сторону на шаг -h1 и вычисляем функцию в этой точке. Если значение функции меньше, переносим базисную точку в эту сторону, иначе оставляем её на прежнем месте. Далее аналогично поступаем с координатой x2 и со всеми остальными до xn.

В результате такого исследования наша базисная точка могла сместиться на новое место по n-мерной поверхности, но могла и остаться на месте. Если смещения не произошло, уменьшаем величину шага по каждой координате и повторяем исследование вокруг этой же точки с меньшим шагом. Практика показывает, что эффективно уменьшать шаг в 10 раз.

Если базисная точка сместилась, далее выполняем поиск по образцу. Образец (т.е. новая базисная точка) выбирается на линии, соединяющей начальную и смещённую базисные точки и в направлении от b1 к b2. Расстояние увеличиваем в два раза. В векторной форме это будет иметь вид:
b3 = b1 + 2(b2 - b1).
Или в общем случае:
bi+2 = bi + 2(bi+1 - bi),
где bi+2 новая базисная точка (образец), вокруг которой снова будем выполнять исследование.

sample

Этот процесс продолжается, пока шаг h не станет меньше значения h < µxmin+ε, где xmin – минимальная из координат последней базисной точки, µ – квадратный корень из машинного эпсилон, ε – желаемая точность результата, которая задаётся при вызове.
Было бы полезно указывать свою величину шага для каждой координаты точки bi, но в данной программе, ради упрощения, шаг одинаков по всем координатам.

[Б12]. Метод Хука-Дживса

 

Вызов программы

call FM25(Fun, X, H, Eps)

Параметры программы

Fun, H, Eps - входные параметры;
Real X[N] - входной и выходной параметр;
Real Fun(Real X(N)) - исследуемая функция;
Real H - начальный шаг;
Real Eps - относительная точность;
Integer N - количество переменных.

 

Пример

!Поиск минимума функции нескольких переменных
!f(x1,x2,...,xn) методом Хука-Дживса
program TestFM25
    use NML
    implicit none
    integer, parameter:: N=3
    real X(N), H, Eps
    data X/4.0, -2.0, 3.0/, H/1.0/, Eps/1.0E-6/
    integer:: Count=0
    integer i
    real FR
    !begin
      call FM25(Fun, X, H, Eps)
      FR=Fun(X); Count=Count-1
      print 10, (X(i), i=1, N)
      print 11, FR
      print 12, Count
10    format(/'The point of a minimum'/<N>(4X,F8.5))
11    format('The value of the function'/'at the point of minimum ',F13.7)
12    format(/'The number of function evaluations ',I3)
    !end

    contains

    real function Fun(X)
    real, intent(in):: x(:)
      Fun=(X(1)-2.0)**2+(X(2)-5.0)**2+(X(3)+2.0)**4
      Count=Count+1
      return
    end function Fun

end program TestFM25

!=========== Вывод программы
 The minimum
  X[1]= 2.000000  X[2]= 5.000000  X[3]=-2.000000
 The minimum of the function is  0.0000000E+00
 The number of iterations          30

The point of a minimum
     2.00000     5.00000    -2.00000
The value of the function
at the point of minimum     0.0000000

The number of function evaluations  91

 

Вернуться к оглавлению    Скачать FM25

FM28

Программа FM28 вычисляет минимум функции нескольких переменных f(x1,x2,...,xn) методом Нелдера-Мида.
В методе Нелдера-Мида используются только значения самой функции, без использований частных производных (градиента) функции.

Регулярным симплексом в n-мерном пространстве называется множество (n+1) равноудалённых точек. Так, для двухмерного пространства это будет равнобедренный треугольник, для трёхмерного – правильный тетраэдр. Идея метода заключается в следующем. Исходя из некоторой начальной точки и начального шага строится (n+1)-мерный симплекс, произвольной формы, не обязательно регулярный. В вершинах симплекса вычисляются значения функции f1= f(x1), f2= f(x2), …, fn= f(xn). Находим наименьшее значение функции fl, наибольшее значение fh, значение, предшествующее наибольшему fg и соответствующие им точки xl, xh, xg. Выбираем числа α, β, γ (альфа, бета, гамма) как коэффициенты для операций отражения, сжатия, растяжения соответственно. Дальнейшие наши действия представлены в виде следующего алгоритма.

schema

1. Вычисляем центр тяжести xw как среднее арифметическое всех точек, за исключением точки xh.

summa

2. Выполняем отражение. Точку xr вычисляем как отражение точки xh относительно точки xw.
xr = - αxh + (1+α)xw
Вычисляем fr= f(xr) и в зависимости от значения fr выбираем один из вариантов:
 ─ Если fr < fl, то мы двигаемся в правильном направлении и выполняем растяжение. Продвигаем xr далее на расстояние γ. xe = γxr + (1 – γ)xw и вычисляем fe= f(xe). Если fe < fl то полагаем xh=xe, иначе (fe > fl) мы продвинулись слишком далеко и xh= xr. Переходим к пункту 5.
 ─ Если fr находится между fl и fg (fl < fr < fg) полагаем xh= xr. Переходим к пункту 5.
 ─ Если fr находится между fg и fh (fg < fr < fh) то меняем местами точки xr и xh. Переходим к пункту 3.

3. fr > fh. Выполняем сжатие. Строим точку xc = βxh + (1 - β)xw и вычисляем fc= f(xc). В зависимости от значения fc выбираем один из вариантов:
 ─ Если fc < fh то присваиваем xh= xc и идём к пункту 5.
 ─ Если fc > fh то это означает, что исходные точки были самыми удачными. Переходим к пункту 4.

4. Совершаем сжатие симплекса. Делим пополам расстояния от каждой точки симплекса до точки с наименьшим значением функции xl, т.е. все xi заменяются на xi + ½(xi - xl), i≠l

5. Завершающий шаг – проверка критерия завершения. Если критерий завершения не выполняется, то совершаем нужные переобозначения точек и значений функции в этих точках и возвращаемся к пункту 1.

В качестве критерия завершения используем среднеквадратичное отклонение значений функции в точках симплекса σ от их среднеарифметического значения.
Вычисления прекращаем, если σ окажется меньше заданного значения ε.

sigma

Коэффициенты для операций отражения, сжатия, растяжения Недлер и Мид рекомендуют брать α=1, β=0.5, γ=2.

[Б12]. Метод Нелдера-Мида

 

Вызов программы

call FM28(Fun, X, H, Alpha, Beta, Gamma, Eps)

Параметры программы

Fun, H, Alpha, Beta, Gamma, Eps - входные параметры;
Real X[N] - входной и выходной параметр;
Real Fun(Real X(N)) - исследуемая функция;
Real H - начальный шаг;
Real Alpha, Beta, Gamma - коэффициенты отражения, сжатия и растяжения;
Real Eps - абсолютная точность;
Integer N - количество переменных.

 

Пример

!Поиск минимума функции нескольких переменных
!f(x1,x2,...,xn) симплексным методом Нелдера-Мида
program TestFM28
    use NML
    implicit none
    integer, parameter:: N=2
    real X(N), H, Eps
    data X/1.5, 2.0/, H/0.5/, Eps/1.0E-5/
    real Alpha, Beta, Gamma  !коэффициенты сжатия
    data Alpha/1.0/, Beta/0.5/, Gamma/2.0/
    integer:: Count=0
    integer i
    real FR
    !begin
      call FM28(Fun, X, H, Alpha, Beta, Gamma, Eps)
      FR=Fun(X); Count=Count-1
      print 10, (X(i), i=1, N)
      print 11, FR
      print 12, Count
10    format(/'The point of a minimum'/<N>(4X,F8.5))
11    format('The value of the function'/'at the point of minimum ',F13.7)
12    format(/'The number of function evaluations ',I3)
    !end

    contains

    real function Fun(X)
    real, intent(in):: x(:)
      Fun=100.0*(X(2)-X(1)*X(1))**2+(1.0-X(1))**2
      Count=Count+1
      return
    end function Fun

end program TestFM28

!=========== Вывод программы
 The minimum
  X[1]= 1.000634  X[2]= 1.001357
 The minimum of the function is  1.1944242E-06
 The number of iterations          37

The point of a minimum
     1.00063     1.00136
The value of the function
at the point of minimum     0.0000012

The number of function evaluations 108

 

Вернуться к оглавлению    Скачать FM28

FM31

Программа FM31 вычисляет минимум функции нескольких переменных f(x) = f(x1,x2,...,xn) методом Дэвидона-Флетчера-Пауэлла.
В методе Давидона-Флетчера-Пауэлла в ходе вычислений используются значения самой функции и её частных производных (градиента).

Любую функцию f(x) в случае её непрерывности и непрерывности её производных можно разложить в ряд Тейлора в окрестности точки x0:
f(x - x0) - f(x0) = (x - x0)Tg(x0) + ½(x - x0)TG(x0)(x - x0) + ...
где g(x0) = ∇f(x0) - градиент функции, G(x0) - её матрица Якоби (матрица Гессе), градиент и матрица Якоби вычислены в точке x0.

В качестве разумной аппроксимации функции f(x) в окрестности точки минимума можно принять квадратичную функцию φ(x):
φ(x) = f(x0) + (x - x0)T∇(x0) + ½(x - x0)TG(x0)(x - x0)
Если её минимум находится в точке xm, то
∇φ(xm) = ∇f(x0) + G(x0)(xmx0) = 0
откуда
xm = x0G-1(x0)∇f(x0) = x0G-1(x0)g(x0)
Таким образом, для нахождения точки минимума можно предложить следующий итерационный процесс:
xi+1 = xiG-1(xi)g(xi)
или, в более общем виде
xi+1 = xi – λiG-1(xi)g(xi)
где длина шага λi определяется одномерной минимизацией в направлении G-1(xi)g(xi)

Метод Ньютона-Рафсона поиска минимума функции f(x) основан на последнем уравнении и требует вычисления и обращения матрицы Якоби на каждом шаге (т.е. при каждой итерации). Этот процесс весьма трудоёмкий в случае, если количество переменных n больше двух. В методе Давидона-Флетчера-Пауэлла на i-ом шаге вычисляется одномерный минимум в направлении –Higi = -Hig(xi), где Hi – симметричная положительно-определённая матрица, которая в конце вычислений становится равной -G-1(xm). Т.е. в этом методе удаётся избежать как вычисления матрицы Якоби, так и её обращения.

Для начала вычислений выбираем точку x0, матрицу H полагаем равной единичной матрице. Далее выполняем следующий итерационный процесс.
1. Выбираем направление поиска di = - Higi
2. Находим коэффициенты λi, минимизирующие функцию f(xi + λidi), что эквивалентно одномерному поиску минимума вдоль прямой xi + λidi
3. Полагаем vi = λidi и xi+1 = xi + vi
4. Находим f(xi+1) и gi+1. Если величины |gi+1| и |vi| достаточно малы, то завершаем процесс.
5. Полагаем ui = gi+1gi
6. Обновляем матрицу H:
Hi+1 = Hi + Ai + Bi
где Ai = viviT/(viTui), Bi = -HiuiuiTHi/(uiTHiui)
7. Увеличиваем i на 1 и возвращаемся к пункту 1.

Для минимизации функции ψ(xi + λdi) по параметру λ (т.е. одномерного поиска минимума f(x) в направлении xi + λdi) применяется кубическая интерполяция по значениям функции и её производной в точках p = xi и q = xi + λdi

Минимум достигается в точке r и равен [Б12]

где z = 3(ψp – ψq)/q + Gp + Gp и w = (z2 – GpGp)½

[Б12]. Метод Давидона-Флетчера-Пауэлла

 

Вызов программы

call FM31(Fun, Grad, X, Eps, Delta)

Параметры программы

Fun, Grad, X, Eps, Delta - входные параметры;
Real X[N] - входной и выходной параметр;
Real Fun(Real X[N]) - исследуемая функция;
Real Grad(Real X[N], Real G[N]) - функция вычисления градиента;
Real Eps - абсолютная точность вычисления градиента;
Real Delta - абсолютная точность вычисления аргумента;
Integer N - количество переменных.

 

Пример

!Поиск минимума функции нескольких переменных
!f(x1,x2,...,xn) методом Давидона-Флетчера-Пауэлла
program TestFMIN
    use NML
    implicit none
    integer, parameter:: N=4
    real X(N), Eps, Delta
    data X/3.0, -1.0, 0.0, 1.0/, Eps/5.0E-5/, Delta/1.0E-5/
    integer:: Count=0
    integer i
    real FR
    !begin
      call FM31(Fun, Grad, X, Eps, Delta)
      FR=Fun(X); Count=Count-1
      print 10, (i, X(i), i=1, N)
      print 11, FR
      print 12, Count
10    format(/'The point of a minimum'/<N>(2X,'X[',I1,']=',F10.7))
11    format('The value of the function'/'at the point of minimum ',E14.7)
12    format(/'The number of function evaluations ',I3)
    !end

    contains

    real function Fun(X)
    real, intent(in):: x(:)
    real Z
      Z=(x(1)+10.0*x(2))**2+5.0*(x(3)-x(4))**2
      Z=Z+(x(2)-2.0*x(3))**4+10.0*(x(1)-x(4))**4
      Fun=Z
      Count=Count+1
      return
    end function Fun

    real function Grad(X, G)
    real, intent(in):: X(:)
    real:: G(:)
    real GG
    integer i, n
      n=size(X); GG=0.0
      G(1)=2.0*(X(1)+10.0*X(2))+40.0*(X(1)-X(4))**3
      G(2)=20.0*(X(1)+10.0*X(2))+4.0*(x(2)-2.0*x(3))**3
      G(3)=10.0*(X(3)-X(4))-8.0*(x(2)-2.0*x(3))**3
      G(4)=-10.0*(x(3)-x(4))-40.0*(X(1)-X(4))**3
      do i=1,n; GG=GG+G(i)*G(i); enddo
      Grad=sqrt(GG)
      return
    end function Grad

end program TestFMIN

!=========== Вывод программы
The point of a minimum
  X[1]= 0.0022939  X[2]=-0.0002294  X[3]= 0.0013408  X[4]= 0.0013402
The value of the function
at the point of minimum  0.8188328E-10

The number of function evaluations 129

 

Вернуться к оглавлению    Скачать FM31

FM34

Программа FM34 вычисляет минимум функции нескольких переменных f(x) = f(x1,x2,...,xn) методом Флетчера-Ривса

Направление градиента это направление, в котором функция возрастает наиболее быстро. Следовательно, в противоположное направление будет показывать наибольшее убывание функции. Возьмем произвольную точку x0 и совершим перемещение из этой точки в направлении –∇f(x0) на расстояние d. Выполним теперь одномерный поиск и найдём то расстояние d, на котором наша функция f(x) примет минимальное значение. Это будет следующая точка x1 в приближении к минимуму функции f(x). Таким образом, можно построить итерационный процесс
xi+1 = xi – di∇f(xi)
где di – параметр d (расстояние), на шаге i минимизирующий функцию
φ(d) = f(xi – di∇f(xi))
Мы получили метод наискорейшего градиентного спуска.

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

В окрестностях минимума разложение в ряд Тейлора любой функции можно аппроксимировать квадратичной функцией (если её вторые производные не равны нулю):
F(x) = a + xTb + ½xTGx
где a – константа, b – постоянный вектор, G – симметричная положительно-определённая матрица. Минимум F(x) в точке xm достигается при условии
∇F(xm) = b + Gxm = 0
откуда xm = - G-1b

Два направления p и q в n-мерном пространстве называются сопряженными, если их скалярное произведение относительно симметричной положительно-определённой матрицы G (порядка n×n)
pTGq = 0
В n-мерном пространстве n взаимно-сопряженных направлений являются линейно-независимыми и могут выступать в качестве базиса этого пространства. Можно показать [Б12], что если поиск минимума квадратичной функции F(x) производить по взаимно-сопряженным направлениям, то он будет найден не более чем через n шагов.

Алгоритм поиска можно составить так. В качестве первого направления из начальной точки x0 возьмём направление наискорейшего спуска
d1 = - ∇F(x0) = - g1
и найдём значение λ1, минимизрующую функцию φ(λ) = f(x1 + λd1). Положим
x2 = x1 + λ1d1
и проведём минимизирующий поиск в направлении d2, сопряжённом с d1 и найдём
x3 = x2 + λ2d2
и т.д. до xn. В случае квадратичной функции F(x) минимум должен быть найден. Данный метод применим и к произвольной функции, если поиск производится вблизи точки минимума, правда, в этом случае процесс может не завершиться на n-ом шаге. Флетчер и Ривс полагают, что в этой ситуации после построения всех сопряжённых направлений должен быть проведён рестарт и каждое n-ое направление поиска должно быть направлением наискорейшего спуска.

Для начала вычислений выбираем точку x1. Далее выполняем следующий итерационный процесс.
1. Определяем направление наискорейшего спуска di = - gi
2. Находим значение λ1, минимизирующее функцию f(xi + λdi)
3. Полагаем xi+1 = xi + λ1di
4. Если x+1i точка минимума, то заканчиваем вычисления, иначе переходим к пункту 5.
5. Если i = n, то полагаем i = 1 и переходим к пункту 1, иначе к пункту 6.
6. Вычисляем k = |gi+1|2/|gi|2; полагаем di+1 = - kdi
7. Увеличиваем i на 1 и возвращаемся к пункту 2.

Для вычисления минимума функции f(xi + λdi) применяем кубическую интерполяцию по значениям функции и её градиента в двух точках (см. FM31).

[Б12]. Метод Флетчера-Ривса

 

Вызов программы

call FM34(Fun, Grad, X, Eps)

Параметры программы

Fun, Grad, X, Eps, Delta - входные параметры;
Real X[N] - входной и выходной параметр;
Real Fun(Real X[N]) - исследуемая функция;
Real Grad(Real X[N], Real G[N]) - функция вычисления градиента;
Real Eps - абсолютная точность вычисления аргумента;
Integer N - количество переменных.

 

Пример

!Поиск минимума функции нескольких переменных
!f(x1,x2,...,xn) методом Флетчера-Ривса
program TestFMIN
    use NML
    implicit none
    integer, parameter:: N=3
    real X(N), Eps
    data X/9.0, -7.0, 11.0/, Eps/5.0E-5/
    integer:: Count=0
    integer i
    real FR
    !begin
      call FM34(Fun, Grad, X, Eps)
      FR=Fun(X); Count=Count-1
      print 10, (i, X(i), i=1, N)
      print 11, FR
      print 12, Count
10    format(/'The point of a minimum'/<N>(2X,'X[',I1,']=',F9.6))
11    format('The value of the function'/'at the point of minimum ',F13.7)
12    format(/'The number of function evaluations ',I3)
    !end

    contains

    real function Fun(X)
    real, intent(in):: x(:)
    real Z
      Z=3.0*(X(1)-1.0)**2+2.0*(X(2)-2.0)**2+(X(3)-3.0)**2
      Fun=Z
      Count=Count+1
      return
    end function Fun

    real function Grad(X, G)
    real, intent(in):: X(:)
    real:: G(:)
    real GG
    integer i, n
      n=size(X); GG=0.0
      G(1)=6.0*(X(1)-1.0)
      G(2)=4.0*(X(2)-2.0)
      G(3)=2.0*(X(3)-3.0)
      do i=1,n; GG=GG+G(i)*G(i); enddo
      Grad=sqrt(GG)
      return
    end function Grad

end program TestFMIN

!=========== Вывод программы
The point of a minimum
  X[1]= 1.000000  X[2]= 2.000000  X[3]= 3.000000
The value of the function
at the point of minimum     0.0000000

The number of function evaluations   7

 

Вернуться к оглавлению    Скачать FM34