Разработка алгоритма

Из предыдущего раздела можно выделить следующие этапы, которые должна выполнять программа:

1 Определить функцию;

2 Ввести приближение и точность;

3 Найти градиент;

4 Найти на направлении против градиента точку, соответствующую минимуму функции и лежащую на градиенте;

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

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

Градиент – это производная скалярной функции, определенной на векторном пространстве. Для численного вычисления градиента будем пользоваться следующей формулой [1, стр.72]

, (2.5)

где h – некоторое положительное число.

Следует подходить ответственно к выбору числа h, так как слишком большие и слишком маленькие его значения сильно искажают результат. Для функции (2.1) можно брать шаг 0,9 или 1. Данные результаты получены опытным путём.

Формула (2.4) очень сложна, поэтому поиск минимума на направлении антиградиента лучше всего проводить численно. В данной реализации осуществлять это будем по следующему алгоритму:

1 Задаем начальную опорную точку;

2 Делаем какой-то шаг относительно опорной точки в сторону антиградиента;

3 Решаем задачу одномерной оптимизации (поиск минимума) в промежутке от опорной точки, до точки, в которую мы попали после отступа.

4 Если в результате оптимизации минимумом оказалась граница данного промежутка, то делаем опорной точкой эту границу и возвращаемся в пункт 2, в противном случае результат – это и есть минимум.

Одним из простых алгоритмов одномерной оптимизации является троичный поиск [5]. При троичном поиске некоторый отрезок делится на три равные части, при этом определены четыре точки: исходные границы (внешние точки) и две точки внутри границ. Затем вычисляется значение функции во внутренних точках. В зависимости от результата устанавливается новый промежуток, одной границей которого будет внешняя точка, а другой какая-то внутренняя точка. Затем все действия повторяются с новыми границами. Общий принцип троичного поиска показан на рисунке 2.2. Например, новый интервал на рисунке 2.2 будет образован из точки a₀ и x₂, так как минимум лежит между этими точками. Процесс повторяется, пока интервал не станет требуемой длины, задаваемой точностью.

Рисунок 2.2 Метод троичного поиска

Из всего выше изложенного можно представить алгоритм графически (рисунок 2.3).

Рисунок 2.3 Метод наискорейшего спуска

Алгоритм троичного поиска представим в виде структурограммы (рисунок 2.4).

Ввод X₁{x1_i}_0…_n, X₂{x2_i}_0…_n, eps

leftTh{a_i}_0…n=(X₁×2+X₂)/3

rightTh{b_i}_0…n=(X₁+2×X₂)/3

f( leftTh )<f( righTh )

да

нет

X₂=rightTh

X₁=leftTh

length=0

для i от 0

length=length+(x2_i-x1_i)²

до n

length=

пока length>eps

Вывод (X₁+X₂)/2

Рисунок 2.4 Метод троичного поиска

Программирование

Ниже представлен листинг программы, реализующей метод наискорейшего спуска. Так как программа содержит очень большое число мелких функций, то все они были вынесены в отдельный файл «optim.cpp», претендующий на звание «библиотеки». Данный файл вставляется в листинг основной программы через свой заголовочный файл «optim.h». Исходный код основной программы хранится в «minimize.cpp».

// Файл: «minimize.cpp»

//------------------------------------------------------

#include <vcl.h>

#include <stdio.h>

#include <conio.h>

#pragma hdrstop

// prototypes

double func(int k,...);

double f(double *x);

#include "optim.h"

//------------------------------------------------------

double func(int k,...)

{

va_list par;

double *pp;

int i;

if(_var==NULL)

{

_inp=-2;

return 0.;

}

va_start(par,k);

pp=va_arg(par,double*);

for(i=0;i<k;i++)

{

_var[i]=*pp;

pp++;

}

va_end(par);

return f(_var);

}

double f(double *x)

{

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

}

#pragma argsused

int main(int argc, char* argv[])

{

int m=2,i;

double *dot=NULL, eps;

if((dot=set_array(m))==NULL)

{

printf("Failed to allocate memory for operational array\n");

return -1;

}

printf("Input the initital estimate\n Size of the estimate is %d\n",m);

for(i=0;i<m;i++)

{

printf("Variable %d: ",i+1);

scanf("%lf",&dot[i]);

}

printf("Input the precision: ");

scanf("%lf",&eps);

gradientLowering(func,m,dot,eps);

printf("Minimum of the function searched by gradient lowering method\n");

printf("[");

for(i=0;i<m;i++)

{

printf("%8.4lf",dot[i]);

if(i!=m-1)

printf(";");

}

printf("]\n");

free(dot);

getch();

return 0;

}

//------------------------------------------------------

// Файл-заголовок «optim.h»

//------------------------------------------------------

#include <math.h>

#include <stdlib.h>

#include <stdarg.h>

double *_lt=NULL;

double *_rt=NULL;

double *_var=NULL;

int _inp=0;

int setTempMas(int m);

void unsetTempMas();

double * set_array(int m);

int shEqu(int k,double *dot1,double *dot2);

int getMin(double (*func)(int k,...),int m, double *left,double *right,double *res,double eps);

double derOfFunc(double (*func)(int k,...),int m,double *dot,int p,double h);

int walk(double (*func)(int k,...),int m,double *sDot, double *eDot, double *tDot, double *h, double eps);

void gradientLowering(double (*func)(int k,...),int m, double *sDot, double eps);

void alarm();

#include "optim.cpp"

//------------------------------------------------------

// Файл «optim.cpp»

//------------------------------------------------------

void alarm()

{

switch(_inp)

{

case -1:printf("\nError 1: failed to allocate memory for operational arrays\n");

break;

case -2:printf("\nError 2: there was an error reading the arguments of the function f(x1,x2...xn)\n");

break;

case -4:printf("\nError 3: failed to allocate memory in gradientLowering()\n");

break;

case -6:printf("\nError 4: happened breakpoint in derOfFunc()\n");

break;

default:break;

}

unsetTempMas();

}

void gradientLowering(double (*func)(int k,...),int m, double *sDot, double eps)

{

double *chkPoint=NULL, *eDot=NULL, *tDot=NULL, *grad=NULL;

int i,check;

chkPoint=set_array(m);

eDot=set_array(m);

tDot=set_array(m);

grad=set_array(m);

if(chkPoint==NULL || eDot==NULL || tDot==NULL || grad==NULL)

_inp=-4;

if(setTempMas(m)!=-1 && _inp>0)

{

for(i=0;i<m;i++)

chkPoint[i]=sDot[i];

while(1)

{

check=0;

for(i=0;i<m;i++)

{

grad[i]=derOfFunc(func,m,sDot,i+1,0.9);

grad[i]=-grad[i];

}

walk(func,m,sDot,eDot,tDot,grad,eps);

for(i=0;i<m;i++)

{

sDot[i]=tDot[i];

if(fabs(sDot[i]-chkPoint[i])<eps)

check=1;

else

check=0;

}

if(check)

break;

for(i=0;i<m;i++)

chkPoint[i]=tDot[i];

}

free(chkPoint);

free(eDot);

free(tDot);

free(grad);

alarm();

}

int walk(double (*func)(int k,...),int m,double *sDot, double *eDot, double *tDot, double *h, double eps)

{

int i,t;

if(sDot==NULL || eDot==NULL || tDot==NULL || h==NULL)

return -1;

while(1)

{

for(i=0;i<m;i++)

eDot[i]=sDot[i]+h[i];

getMin(func,m,sDot,eDot,tDot,eps);

if((t=shEqu(m,sDot,tDot)) && t!=-1)

for(i=0;i<m;i++)

{

sDot[i]=tDot[i];

h[i]=-h[i];

}

else

if((t=shEqu(m,eDot,tDot)) && t!=-1)

for(i=0;i<m;i++)

sDot[i]=eDot[i];

else

break;

for(i=0;i<m;i++)

h[i]*=2;

}

if(t==-1)

return -1;

else

return 0;

}

int setTempMas(int m)

{

_lt=set_array(m);

_rt=set_array(m);

_var=set_array(m);

if(_lt==NULL || _rt==NULL || _var==NULL)

{

_inp=-1;

unsetTempMas();

return -1;

}

_inp=1;

return 0;

}

void unsetTempMas()

{

if(_lt!=NULL)

free(_lt);

if(_rt!=NULL)

free(_rt);

if(_var!=NULL)

free(_var);

_lt=_rt=_var=NULL;

_inp=0;

}

double *set_array(int m)

{

double *res=NULL;

res=(double *)calloc(m,sizeof(double));

if(res==NULL)

return NULL;

return res;

}

int shEqu(int k,double *dot1,double *dot2)

{

int i;

if(dot1==NULL || dot2==NULL)

return -1;

for(i=0;i<k;i++)

if(dot1[i]!=dot2[i])

return 0;

return 1;

}

int getMin(double (*func)(int k,...),int m, double *left,double *right,double *res,double eps)

{

double length=0;

int i;

for(i=0;i<m;i++)

length+=(right[i]-left[i])*(right[i]-left[i]);

length=sqrt(fabs(length));

if(length<eps)

{

for(i=0;i<m;i++)

res[i]=(left[i]+right[i])/2;

return 0;

}

for(i=0;i<m;i++)

{

_lt[i]=(left[i]*2+right[i])/3;

_rt[i]=(left[i]+2*right[i])/3;

}

if(func(m,_lt)<func(m,_rt))

return getMin(func,m,left,_rt,res,eps);

else

return getMin(func,m,_lt,right,res,eps);

}

double derOfFunc(double (*func)(int k,...),int m,double *dot,int p,double h)

{

int i;

if(_inp==0)

setTempMas(m);

if(dot==NULL)

_inp=-6;

if(p<=0 || p>m)

_inp=-6;

if(_inp<0)

return -1.;

if(h<=0.)

h=0.0001;

for(i=0;i<m;i++)

{

_lt[i]=dot[i];

_rt[i]=dot[i];

}

_lt[p-1]=_lt[p-1]-h;

_rt[p-1]=_rt[p-1]+h;

return (func(m,_rt)-func(m,_lt))/2*h;

}

//------------------------------------------------------

Основная программа представляет собой вызов функции gradientLowering, которая выполняет наискорейший спуск. Тело функции хранится в «optim.cpp». Функции необходимо знать математическую функцию, для которой нужно найти минимум, число аргументов в ней, массив с начальным приближением и точность вычисления. Размер массива с приближением должен соответствовать (или, по крайней мере, быть не меньше) числа переменных в функции, так как весь упор делается на это число. В случае удачи в массив с приближением запишется минимум. Функция ничего не возвращает, однако, в структуре используется аппарат по отлову всевозможных ошибок, кроме ошибок, связанных с вычислением функции.

Для создания массива с приближением используется функция set_array, тело которой хранится в «optim.cpp». Функции необходимо передавать желаемый размер массива. Функция возвращает указатель на массив в случае удачи и NULL значение в противном случае. Данная функция используется в некоторых функциях файла «optim.cpp» и может использоваться сама по себе.

Для объявления арифметической функции используется функция func. Для простоты, в данной реализации все арифметические операции функции (2.1) помещены в отдельную функцию f, который нужно передавать указатель на массив с аргументами. Функция func в данной реализации фактически генерирует массив, который передаёт потом f. Благодаря особенности func, исходная арифметическая формула может не иметь ограничений на число аргументов, главное знать их число. В основе func лежит механизм макросредств, обеспечивающих доступ к переменному списку параметров [2, стр. 249]. Данной функции необходимо передавать число аргументов (это обязательный параметр) и указатель на вектор с аргументами, который должен соответствовать числу аргументов (или, по крайней мере, быть не меньше этого числа). Функция перехватывает ошибку, когда ей передают NULL значение в качестве указателя.

Перейдем к непосредственному описанию функций файла «optim.cpp». В данном файле хранятся следующие функции:

- int setTempMas(int m);

- void unsetTempMas();

- double * set_array(int m);

- int shEqu(int k,double *dot1,double *dot2);

- int getMin(double (*func)(int k,...),int m, double *left,double *right,double *res,double eps);

- double derOfFunc(double (*func)(int k,...),int m,double *dot,int p,double h);

- int walk(double (*func)(int k,...),int m,double *sDot, double *eDot, double *tDot, double *h, double eps);

- void gradientLowering(double (*func)(int k,...),int m, double *sDot, double eps);

- void alarm().

Функция setTempMas инициализирует внутренние временные массивы, которыми пользуются некоторые функции файла «optim.cpp». Имена этих массивов зарезервированы и не могут использоваться, если подключен данный файл. Данные массивы имеют имена _lt, _rt, _var. Массивами _lt и _rt пользуются в своих целях функции getMin и derOfFunc; массивом _var пользуется func.

Функция unsetTempMas освобождает память от временных массивов.

Функция shEqu сравнивает две точки и возвращает 1, если это одна и та же точка, в противном случае 0. Используется в связке с функцией walk.

Функция getMin выполняет троичный поиск по алгоритму рисунка 2.4. В данной функции реализована рекурсия. Возвращает 0, если минимум найден, при этом результат записывается в массив res. Используется в связке с функцией walk.

Функция derOfFunc возвращает значение производной в точке. Функции необходимо передавать дифференцируемую функцию, число переменных, точку, в которой рассчитывается производная, номер переменной, по которой нужно дифференцировать, шаг. Здесь подразумевается, что переменные нумеруются от единицы, поэтому необходимо иметь это в виду при передаче номера. Если функции передан отрицательный или нулевой шаг, то для вычисления производной используется шаг 0,0001, рекомендуемый в некоторых источниках [4].

Функция walk фактически реализует поиск минимума на антиградиентном направлении. Функции необходимо передавать арифметическую функцию, число переменных, начальную точку, два временных массива (во второй в случае удачи записывается результат), шаговый массив (здесь это градиент) и точность.

Функция alarm служит для «отлова» ошибок, которые могут произойти в процессе работы программы, кроме ошибок, связанных с вычислением функции. Для оценки ситуации она пользуется глобальной переменной _inp, чьё имя зарезервировано и не может быть использовано, если подключен «optim.h». Все ошибки имеют свои коды, по которым alarm их вычисляет и выдаёт предупреждающие сообщения. В задачи alarm также входит очистка временных массивов функции setTempMas.

В данной реализации могут быть обнаружены ошибки указанные в таблице 2.

Таблица 2

Код ошибки	Сообщение для пользователя	Возможные причины
	Не удалось выделить память под временные массивы.	Передан ненатуральный размер при вызове функции setTempMas, либо недостаточно памяти.
	Не удалось прочитать аргументы функции.	Функции func передан NULL указатель.
	Не удалось выделить память в функции gradientLowering().	Передан ненатуральный размер числа переменных, либо недостаточно памяти.
	Аварийный останов при вычислении производной.	Функции derOfFunc передан NULL указатель вместо реальной точки, либо передан номер несуществующей переменной.

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

Анализ результатов

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

Теперь испытаем программу при пяти приближениях, данных в таблице 3. Для наглядности работы алгоритма, приведем рисунок, похожий на рисунок 2.1, для этих приближений.

Правильность получаемого ответа установим с помощью Mathcad, для чего воспользуемся внутренней функцией Minimize.

На рисунке 2.7 показаны две изоклины (линии в виде эллипсов) и пути следования точки при выполнении алгоритма рисунка 2.3. Сравнивая ответ, данный Mathcad с ответами таблицы 3, можно убедиться, что программа выдаёт достаточно правильный ответ. Также можно заметить, что при некоторых приближения метод сходится к ответу достаточно быстро, а при некоторых несколько медленнее.

Рисунок 2.5 Результат работы программы оптимизации

Таблица 3

№ п/п	Приближение	Выводимый ответ	Число затраченных итераций	Точность
x₁	x₂	x₁	x₂
			7,7987	-2,3996	0,0001
		-500,81	7,8006	-2,4001
	300,58	481,26	7,8010	-2.4003
			7,8006	-2,4001
	-50	-1	7,7997	-2,4000

Рисунок 2.6 Результат работы внутренней функции Mathcad

Рисунок 2.7 Пути следования точки при различных приближениях

ЗАКЛЮЧЕНИЕ

Выше было показано, как использование ЭВМ позволяет сэкономить время и силы для решения различного рода тривиальных задач, которые раньше приходилось решать человеку.

Представленные в курсовой работе коды не являются окончательными и могут быть многократно модифицированы.

Дата добавления: 2015-12-17; просмотров: 20; Мы поможем в написании вашей работы!

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

⇐ Предыдущая 1 2 3 45

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