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

Размещено на http://www.allbest.ru

Министерству образования и науки Российской Федерации

Федеральное государственное автономное образовательное учреждение

высшего профессионального образования

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

(ВолГУ)

Институт математики и информационных технологий

Лабораторная работа №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 ,и одинаковых точностях (, ), при центральной РП и правой РП (ЦРП/ПРП).

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

n	h
4	49 / 49*	48 / 55*	48 / 172	48 / 51
5	31 / 71*	31 / 49*	31 / 34	31 / 26
6	24 / 37*	24 / 31*	24 / 28	24 / 26
7	37 / 71*	35 / 42*	35 / 40	35 / 49
8	40 / 49*	43 / 64*	43 / 59	43 / 55
9	30 / 46*	30 / 36*	30 / 48	30 / 37
10	42 / 31*	42 / 50*	42 / 45	42 / 42
11	36 / 34*	36 / 34*	36 / 37	36 / 41
12	38 / 49*	39 / 48*	39 / 39	39 / 37
13	37 / 53*	38 / 52*	38 / 43	38 / 40
14	37 / 43*	37 / 57*	37 / 43	37 / 51
15	41 / 46*	42 / 75*	42 / 48	42 / 43
16	35 / 49*	35 / 66*	35 / 37	35 / 35
17	41 / 52*	40 / 51*	40 / 55	40 / 41
18	41 / 37*	44 / 90*	42 / 45	43 / 44

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

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

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

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

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

n	h
5	22 / 21*	20 / 36*	20 / 24	20 / 21
6	13 / 25*	15 / 19*	15 / 26	15 / 18
7	20 / 29*	19 / 29*	19 / 23	19 / 19

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

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

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

#include <stdio.h>

#include <math.h>

#include <windows.h>

#define EPS 0.000000002

#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(" %s\n",ToRus("разностная схема нахождения градиента"));

printf("%s:\tx0 = (",ToRus("Начальная точка"));

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

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

}

printf("),\t f(x0) = %.6lf\n",function(x));

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

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

printf("%s:\th = %lf\n",ToRus("Шаг в разностной производной"),h);

g_old = gradient(x,h,ds);

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

printf("\n%s\t\t\t %s\t\t\t %s\n",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=%.6lf\n",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)<EPS){

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

printf("\t%s\n",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=%.6lf\n",yI);

k ++;

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

puts("---------------------------------------------------------------------------");

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

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

printf(") = %.6lf\n",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<nv;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<nv;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<nv;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<nv;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<nv;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<nv;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<nv;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<nv;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<nv;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;

}

Размещено на Allbest.ru