Double f, f1, f2, f3, x0, x1, x2, x3, L, t1, t2,

tolerance =0.00001;

t 1 = 0.3819660113; // (1 – tau ) – левая точка внутри

t 2 = (1.0 – t 1); // ( tau )    – правая точка внутри

// Значения функции на границах и внутри интервала

x0 = *a;           f0 = Fx(x0);

x1 = *a + t1*(*b - *a); f1 = Fx(x1);

x2 = *a + t2*(*b - *a); f2 = Fx(x2);

x 3 = * b ;       f 3 = Fx ( x 3);

// Основной цикл поиска

do {

if ( f 1 <= f 2) { //Исключаем отрезок [ x ₂ , x ₃ ]

L = x2 – x0;

x3 = x2; f3 = f2;

x2 = x1; f2 = f1;

x 1 = x 0 + t 1* L ; f = f 1 = Fx ( x 1);

} else {  // Исключаем отрезок [ x ₀ , x ₁ ]

L = x3 – x1;

x0 = x1; f1 = f0;

x1 = x2; f1 = f2;

x2 = x0 + t2*L; f = f2 = Fx(x2);

}

} while (fabs(L) > tolerance);

*a = x0; *b = x3;

return (f);

}

 

Пример 5 . Метод золотого сечения

Рассмотрим задачу из примера 4, в которой требуется минимизировать f ( x )=(100—х)² в интервале 60 ≤ х ≤ 150. Для того чтобы перейти к интервалу единичной длины, проведём замену переменной, положив w = (x - 60) / 90. Таким образом, задача принимает следующий вид: минимизировать f (w)=(40 – 90w)² при ограничении 0 ≤ w ≤ 1.

И т е р а ц и я  1.

I ₁ = (0, 1); L ₁ = 1. Проведём вычисления значений функции:

w ₁ = τ = 0.618;   f ( w ₁ ) = 244.0.

w ₂ =1— τ = τ² = 0.382; f ( w ₂ ) = 31.6.

Так как f ( w ₂ ) < f ( w ₁ ) и w ₂ < w ₁ интервал w ³ w ₁ исключается.

И т е р а ц и я  2.

I ₂ = (0, 0.618); L ₂ = 0.618 = τ. Следующее вычисление значения функции проводится в точке

w ₃ = τ - τ ²= τ (1 - τ) = τ ³ = 0.236; f ( w ₃ ) = 352.

Так как f ( w 3) > f ( w ₂ ) и w ₃ < w ₂, интервал w ≤ w ₃ исключается.

И т е р а ц и я  3.

I₃ = (0.236, 0.618); L₃ = 0.382 = τ².

Следующее вычисление значения функции проводится в точке, расположенной на расстоянии τ ´ (длину полученного интервала) от левой граничной точки интервала, или на расстоянии (1— τ) ´ (длину интервала) от правой граничной точки. Таким образом,

w ₄ = 0.618 - (1 - τ) L ₃ = 0.618 - τ² L ₃  = 0.618 - τ²(τ²) = 0,618 - τ⁴ = 0.472;

f ( w ₄ ) = 6.15.

Так как f ( w ₄ ) < f ( w ₂ ) и w ₄ > w ₂, интервал w ≤ w ₂ исключается.

В результате получен следующий интервал неопределённости: 0.382 ≤ w ≤ 0.618 для переменной w , или 94.4 < x ≤ 115.6 для переменной х.

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

τ ^{N -1} = τ⁵ = 0.09;

что соответствует интервалу длиной 8.1 для переменной х. Для сравнения напомним, что в аналогичной ситуации метод деления интервала пополам привёл к получению интервала длиной 11.5.

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

3.2.3. Квадратичная аппроксимация. Метод Пауэлла.

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

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

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

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

Если задана последовательность точек х₁, х₂, х₃ и известны соответствующие этим точкам значения функции f ₁ , f ₂ , и f ₃, то можно определить постоянные величины a ₀ , a ₁ и а₂ таким образом, что значения квадратичной функции

q ( x ) = a ₀ + a ₁ ( x—х₁) + а₂ (х—х₁) (х—х₂)            (16)

совпадут со значениями f ( x ) в трёх указанных точках. Перейдём к вычислению q ( x ) в каждой из трёх заданных точек. Прежде всего, так как

f ₁ = f ( x ₁ ) = q ( x ₁ ) = a ₀ ;                                           (17)

имеем a ₀ = f₁. Далее, поскольку

f₂ = f (x₂) = q (x₂) = f₁ + a₁ (x₂ – x₁);                       (18)

получаем a ₁ = ( f ₂ — f ₁ ) / ( x ₂ — x ₁ ). Наконец, при х = х₃:

f ₃ = f ( x ₃ ) = q ( x ₃ ) = f ₁ +( f ₂ — f ₁ )/( x ₂ — x ₁ )( x ₃ — x ₁ )+ a ₂ ( x ₃ — x ₁ )(х₃—х₂). (19)

Разрешая последнее уравнение относительно a ₂ , получаем

                 (20)

Таким образом, по трём заданным точкам и соответствующим значениям функции можно оценить параметры а₀, а₁ и а₂ аппроксимирующего квадратичного полинома с помощью приведённых формул.

Если точность аппроксимации исследуемой функции в интервале от x ₁ до х₃ с помощью квадратичного полинома оказывается достаточно высокой, то в соответствии с предложенной стратегией поиска построенный полином можно использовать для оценивания координаты точки оптимума. Напомним, что стационарные точки функции одной переменной определяются путём приравнивания к нулю её первой производной и последующего нахождения корней полученного таким образом уравнения, В данном случае из уравнения

можно получить

Поскольку функция f ( x ) на рассматриваемом интервале обладает свойством унимодальности, а аппроксимирующий квадратичный полином также является унимодальной функцией, то можно ожидать, что величина  окажется приемлемой оценкой координаты точки истинного оптимума х*.

Пример 6. Использование квадратичной аппроксимации

Покажем процедуру оценки координаты точки минимума функции f ( x )=2 x ² +(16/ x ) в интервале 1 ≤ х ≤ 5. Пусть х₁ = 1, х₃ = 5,а х₂ есть средняя точка интервала, то есть х₂ = 3. Вычислим соответствующие значения функции:

f ₁ = 18; f ₂ = 23.33; f ₃ = 53.2.

Для того чтобы найти оценку , необходимо найти значения параметров a ₁ и а₂ аппроксимирующей функции. Имеем:

Подстановка этих значений в формулу для х позволяет получить

х = (3+1)/2—(8/3)/[2 (46/15)] = 1.565.

Точный минимум достигается при х* = 1.5874.

Пауэллом [2] разработан метод, основанный на последовательном применении процедуры оценивания с использованием квадратичной аппроксимации. Схему алгоритма Пауэлла можно описать следующим образом. Пусть х₁ - начальная точка, Δx - выбранная величина шага по оси х.

Шаг 1. Вычислить x ₂ = x ₁ + Δx .

Шаг 2. Вычислить f ( x ₁ ) и f ( x ₂ ).

Шаг 3. Если f ( x ₁ ) > f ( x ₂ ) положить х₃ = x ₁ + 2 Δx . Если f ( x ₁ ) ≤ f ( x ₂ ), тоположить x ₃ = x ₁ - Δx .

Шаг 4. Вычислить f (х₃) по формулам (17)-(20) и найти

F_min = min (f₁, f₂, f₃);

x_min = точка x_i которая соответствует F_min .

Шаг 5. По трем точкам х₁, x ₂ , x ₃ вычислить , используя фор­мулу для оценивания с помощью квадратичной аппроксимации.

Шаг 6. Проверка на окончание поиска.

(а) Является ли разность (F_min - f ( )) достаточно малой?

(б) Является ли разность (x_min - ) достаточно малой?

Если оба условия выполняются, то закончить поиск. В противном случае перейти к шагу 7.

Шаг 7. Выбрать «наилучшую» точку из x_min и  и две точки по обе стороны от неё. Обозначить эти точки в естественном порядке х₁, x ₂ , x ₃ и перейти к шагу 4.

Заметим, что при первой реализации шага 5 границы интервала, содержащего точку минимума, не обязательно оказываются внутри найденного интервала неопределённости. При этом полученная точка  может находиться за точкой х₃. Для того чтобы исключить возможность слишком большого экстраполяционного перемещения, следует провести после шага 5 дополнительную проверку и в случае, когда точка  находится слишком далеко от х₃, заменить  точкой, координата которой вычисляется с учётом заранее установленной длины шага.

Пример 7. Метод Пауэлла

Рассмотрим задачу из примера 6: минимизировать f (х) = 2х² + (16/х).

Пусть начальная точка x ₁ = 1 и длина шага Δx = 1. Для проверки на окончание поиска используются следующие параметры сходимости:

Итерация 1.

Шаг1. х₂ = х₁ + Δx = 2.

Шаг 2. f ( x₁ ) = 18; f ( x ₂ ) = 16.

Шаг 3. f ( x ₁ ) > f ( x ₂ ), следовательно, положить х₃ = 1 + 2 = 3.

Шаг 4. f (х₃) = 23.33; F_min = 16; x_min = х₂.

Шаг 5.

Шаг 6. Проверка на окончание поиска:

(а)

Следовательно, продолжаем поиск.

Шаг 7. Выбираем  как «наилучшую» точку, a x ₁ = 1и х₂ = 2 - как точки, которые её окружают. Обозначаем эти точки в естественном порядке и переходим к итерации 2, которая начинается с шага 4.

Итерация 2

Шаг 4. x ₁ = 1; f ₁ = 18, х₂ = 1.714, f ₂ = 15.210 = F_min; x_min = x ₂; x ₃=2;  f ₃ = 16.

Шаг 5.

Шаг 6. Проверка на окончание поиска:

(а)

(б)

Следовательно, поиск закончен.

Программа, реализующая метод Пауэлла выглядит следующим образом.

#include <math.h>

double Powell(double *a, double *b)

{

// Программа поиска минимума на интервале [ a , b ]

// методом квадратичной аппроксимации Пауэлла.

// Значение функции f ( x ) вычисляется в подпрограмме Fx ( x )

// В результате интервал неопределённости уменьшается и

// его границы доступны в вызывающей программе.

// Также «наверх» передаётся минимальное значение функции.

extern double Fx(double);

double dx, x[4], f[4], xt, ft, dn, nm,

tolerance =0.001;

int j,k, s1, s2, s3;

dx = 0.05; // Шаг

// Начать процесс с первых трёх точек

x[0] = *a;       f[0] = Fx(x[0]);

x[1] = *a + dx;  f[1] = Fx(x[1]);

if (f[0] < f[1]) x[2] = *a - dx; else x[2] = *a + 2.0*dx;

f[2] = Fx(x[2]);

//Первый интерполированный минимум

 dn = ( x [1 ]-x[2])*f[0];

dn += (x[2]-x[0])*f[1] + (x[0]-x[1])*f[2];

 nm = (x[1]*x[1]-x[2]*x[2])*f[0];

nm += (x[2]*x[2]-x[0]*x[0])*f[1];

nm += (x[0]*x[0]-x[1]*x[1])*f[2];

x[3] = nm /(2.0*dn); f[3] = Fx(x[3]);

while(1) {// Бесконечный цикл интерполяции

for (j = 0; j < 3; j++)

for (k = 0;k < 4; k++) {

   if ( f [ j ] > f [ k ] {//Упорядочить x и f

       xt = x[j]; x[j] = x[k]; x[k] = xt;

       ft = f[j]; f[j] = f[k]; f[k] = ft;

   }//Конец оператора if ( )

} // Конец циклов по j и k

if (fabs(x[0]-x[1]) < tolerance) {

---

Дата добавления: 2019-09-13; просмотров: 205; Мы поможем в написании вашей работы!

Поделиться с друзьями:

---

---

Мы поможем в написании ваших работ!