Методы многомерной безусловной минимизации. 
Сравнение правой РП и центральной РП на примере минимизации функции нескольких аргументов методом сопряженных г

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

«Волгоградский государственный университет»

(ВолГУ) Институт математики и информационных технологий Лабораторная работа № 3

по курсу «Численные методы оптимизации»

Методы многомерной безусловной минимизации. Сравнение правой РП и центральной РП на примере минимизации функции нескольких аргументов методом сопряженных градиентов Выполнили:

Студенты 3-го курса группы ПМ-101

Самородов Е. А.

Гусынин О. С.

Болотин А. В.

Принял:

Яновский Т.А.

Волгоград 2013

Задание:

Сравнить правую и центральную разностную производную на примере минимизации функции нескольких аргументов методом сопряженных градиентов. В качестве функции взята функция Пауэлла.

Метод:

Метод сопряженных градиентов (метод Флетчера-Ривса) В основе метода лежит построение направлений поиска минимума, являющихся линейными комбинациями градиента предыдущих направлений поиска. При этом весовые коэффициенты выбираются так, чтобы сделать направления сопряженными относительно матрицы Гессе. Для повышения скорости сходимости метода, в случае не квадратичной используется рестарт: через каждые n циклов направление поиска заменяется на Наряду с начальной точкой и вектором положительных приращений координат, алгоритм метода требует априорного задания параметра точности поиска е > 0 .

Алгоритм:

Функция:

Результаты работы программы:

Был проведен ряд расчетов с целью определения отличий центральной РП и правой РП. Для этой цели в качестве переменных параметров были взяты eps1 и h. Где eps1 — максимальная величина нормы градиента при которой продолжается расчет, h — шаг который используется в РП для нахождения производной функции.

В качестве метода одномерной минимизации взят метод золотого сечения.

При минимизации отрезка использован постоянный шаг 0.01. Для цели работы он не важен, как и точность с которой будет найден alfa методом золотого сечения.

Пример работы программы:

Используем для расчета Центральную РП.

Здесь под МСГ понимается количество итераций пройденных методом сопряженных градиентов, под ОМ — количество итераций при нахождении отрезка минимизации, под ЗС — количество итерации в методе золотого сечения.

Теперь используем Правую РП.

На первый взгляд видно, что разница между ЦРП и ПРП только в количестве итераций. С помощью ЦРП требуемый результат найден чуть быстрее. Теперь изменим eps1 и h. Возьмем eps1 = 0.1, h =0.001. С центральной РП:

С правой РП:

Очевидно, что использование ЦРП, при более грубом шаге h и низкой точности eps1, способствует быстрейшему и более точнейшему нахождению минимума функции в отличие от ПРП.

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

Ниже приведена таблица количеств итерации при различных n и h, и одинаковых точностях (,), при центральной РП и правой РП (ЦРП/ПРП).

многомерный минимизация производная функция

Как видно из таблицы правая РП более чувствительна к изменению шага h чем центральная. Из этого следует, что ЦРП имеет меньшею погрешность нахождения производной, чем ПРЦ.

К тому же, при больших значения h применение ПРП, в данном примере, привило к очень малому изменению аргументов, т. е. к зацикливанию.

В таблице (*) помечено то количество итераций, где применение ПРП привело к зацикливанию, что не позволило достичь заданной точности (из-за этого значения некоторых чисел у ПРП меньше чем у ЦРП).

Можно увидеть, что при h равном и использование ЦРП эффективнее чем ПРП.

Для сравнения ПРП и ЦРП при большом шаге h точность для останова сделаем более грубой. Тогда увидим, что:

Видно, что при h равном и использование ЦРП по-прежнему эффективнее, чем ПРП. А при h равном и, даже при не достижении требуемой точности у ПРП, количество итераций при использовании ЦРП меньше чем при использовании ПРП.

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

Код программы на Си:

#include

#include

#include

#define EPS 0.2

#define f function

#define nv the_number_of_variables

#define CENTRAL 0

#define RIGHT 1

#define LEFT 2

const int the_number_of_variables = 4;

double function (double * x){

return pow (x[0]+10*x[1], 2)+5*pow (x[2]-x[3], 2)+pow (x[1]-2*x[2], 4)+10*pow (x[0]-x[3], 4);

}

char * ToRus (const char text []);

double * gradient (double * x, double h, int ds);

double CentralDS (double * x, double h, int i);

double RightDS (double * x, double h, int i);

double LeftDS (double * x, double h, int i);

double (*DS[3])(double * x, double h, int i) = {CentralDS, RightDS, LeftDS};

double NormaV (double * x, int n);

double * minimi (double * x, double * p, double dx);

double zoloto (double * x, double * p, double a, double b, double eps);

double argminf (double * x, double * g, double eps);

int main (){

int i, k, n, ds;

double *g, *g_old, *x, *p, yI, yII;

double eps1, eps2, h, Sc, Sz, betta, alfa;

x = (double *)malloc (nv*sizeof (double));

p = (double *)malloc (nv*sizeof (double));

k = 1;

n = 6;//цикличность (через сколько итерации произойдет рестарт)

eps1 = 0.001;

eps2 = 1e-6;

h =1e-6;//шаг при вычисления градиента в РС

ds = CENTRAL;//RIGHT

printf («%s «, ToRus («При расчетах использована»));

switch (ds+1){

case 1: printf («%s», ToRus («центральная»));break;

case 2: printf («%s», ToRus («правая»));break;

case 3: printf («%s», ToRus («левая»));break;

}

printf («%sn», ToRus («разностная схема нахождения градиента»));

printf («%s:tx0 = („, ToRus (“ Начальная точка»));

for (i = 0; i < nv; i ++){

x[i] = 1;printf («%.0lf,», x[i]);

}

printf («), t f (x0) = %.6lfn», function (x));

printf («%s:teps1 = %.3lfn», ToRus («Максимальная величина нормы градиента при останове»), eps1);

printf («%s:teps2 = %lfn», ToRus («Точность вычисления аргумента alfa на каждой итерации»), eps2);

printf («%s:th = %lfn», ToRus («Шаг в разностной производной»), h);

g_old = gradient (x, h, ds);

for (i = 0; i < nv; i ++)p[i] = - g_old[i];

printf («n%sttt %sttt %sn», ToRus («МСГ [ОМ/ЗС]»), ToRus («Аргументы»), ToRus («Функция»));

puts («—————————————————————————————————————-»);

printf («%2d «, k);

yI = f (x);

alfa = argminf (x, p, eps2);

for (i = 0; i < nv; i ++){x[i] += alfa*p[i]; }

for (i = 0; i < nv; i ++){printf («x%d=%.4lf «, i+1,x[i]);}

printf («f=%.6lfn», f (x));

k ++;yII = f (x);

do{

g = gradient (x, h, ds);

Sc = 0;

Sz = 0;

for (i = 0; i < nv; i ++){

Sc += g[i]*g[i];

Sz += g_old[i]*g_old[i];

}

betta = Sc/Sz;

if (k%n==1){

betta = 0;

if (fabs (yI-yII)

eps1 = NormaV (g, nv)+1.;

printf («t%sn», ToRus («Остановка по причине черезвычайно малого изменения y (x)»));

}

}

for (i = 0; i < nv; i ++){

g_old[i] = g[i];

p[i] = betta*p[i] - g[i];

}

printf («%2d «, k);

alfa = argminf (x, p, eps2);

for (i = 0; i < nv; i ++){x[i] += alfa*p[i]; }

for (i = 0; i < nv; i ++){printf («x%d=%.4lf «, i+1,x[i]);}

yII = yI;

yI = f (x);

printf («f=%.6lfn», yI);

k ++;

}while (NormaV (g, nv) > eps1);

puts («—————————————————————————————————————-»);

printf («%s: f („, ToRus (“ Итог расчетов»));

for (i = 0; i < nv; i ++){printf («%.4lf,», x[i]);}

printf («) = %.6lfn», f (x));

free (x);

free (g);

free (g_old);

free (p);

return 0;

}

double CentralDS (double * x, double h, int i){

int j;

double * y;

double df;

y = (double *)malloc (nv*sizeof (double));

for (j=0;j

y[j] = x[j];

}

y[i] += 0.5*h;

df = f (y);

y[i] -= h;

df -= f (y);

df /= h;

free (y);

return df;

}

double RightDS (double * x, double h, int i){

int j;

double * y;

double df;

y = (double *)malloc (nv*sizeof (double));

for (j=0;j

y[j] = x[j];

}

y[i] += h;

df = (f (y) — f (x))/h;

free (y);

return df;

}

double LeftDS (double * x, double h, int i){

int j;

double * y;

double df;

y = (double *)malloc (nv*sizeof (double));

for (j=0;j

y[j] = x[j];

}

y[i] -= h;

df = (f (x) — f (y))/h;

free (y);

return df;

}

double * gradient (double * x, double h, int ds){

int i;

if (ds > 2 || ds < 0) ds = 0;

double * g;

g = (double *)malloc (nv*sizeof (double));

for (i=0;i

g[i] = DS[ds](x, h, i);

}

return g;

}

double NormaV (double * x, int n){

int i;

double Sk = 0;

for (i = 0; i < n; i++){

Sk += x[i]*x[i];

}

return sqrt (Sk);

}

double argminf (double * x, double * p, double eps){

double alfa;

double * ab = minimi (x, p,0.01);//0.01 — оптимальный шаг (подобран экспериментально)

alfa = zoloto (x, p, ab[0], ab[1], eps);

free (ab);

return alfa;

}

double * minimi (double * x, double * p, double dx){

int i;

double m=0, a, b, k=1, alfa = 0;

double * xda = (double*)malloc (nv*sizeof (double));

double * xdb = (double*)malloc (nv*sizeof (double));

for (i=0;i

xda[i] = x[i] + alfa*p[i];

xdb[i] = x[i] + (alfa+dx)*p[i];

}

if (f (xda) < f (xdb))m =- 1;

else m = 1;

dx *= m;

a = alfa+dx;

b = alfa+3*dx;

int j=1;

do{

for (i=0;i

xdb[i] = x[i] + b*p[i];

xda[i] = x[i] + a*p[i];

}

k*=2;

b += k*dx;

a = b-2*k*dx;

j++;

}while (f (xda) > f (xdb));

printf («[%2d/», j);

free (xda);

free (xdb);

b -= k*dx;

a = b-2*k*dx;

double *s;

s = (double*)malloc (2*sizeof (double));

if (m==1){ s[0] = a; s[1] = b;}

else {s[0] = b; s[1] = a;}

return s;

}

double zoloto (double * x, double * p, double a, double b, double eps){

int i;

double x1, x2, tau, A, B;

double * xd1 = (double*)malloc (nv*sizeof (double));

double * xd2 = (double*)malloc (nv*sizeof (double));

tau = (sqrt (5.)+1.)/2.;

x1 = a+(b-a)/(tau*tau);

x2 = a+(b-a)/tau;

for (i=0;i

xd1[i] = x[i] + x1*p[i];

xd2[i] = x[i] + x2*p[i];

}

B = f (xd2);

A = f (xd1);

int j=0;

do{

if (A <= B){

b = x2;

x2 = x1;

x1 = a + b — x2;

B = A;

for (i=0;i

xd1[i] = x[i] + x1*p[i];

}

A = f (xd1);

}

else{

a = x1;

x1 = x2;

x2 = a + b — x1;

A = B;

for (i=0;i

xd2[i] = x[i] + x2*p[i];

}

B = f (xd2);

}

j++;

}while (fabs (b-a) > eps);

printf («%2d]t», j);

free (xd1);

free (xd2);

return (b+a)/2.;

}

char * ToRus (const char text []){

char * buf = (char*)malloc (strlen (text)+1);

CharToOem (text, buf);

return buf;

}