Минимизация функции методом сопряженных градиентов

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

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

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

ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«ОРЕНБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ»

Факультет экономики и управления

Кафедра математических методов и моделей в экономике

ОТЧЕТ

по учебной-вычислительной практике

на базе кафедры математических методов и моделей в экономике

ГОУ ОГУ 080116.65.9011.06 П

Руководители от кафедры:

профессор, заведующий кафедры ММиМЭ А.Г. Реннер

ассистент кафедры ММиМЭ Р.М. Шаяхметова

ассистент кафедры ММиМЭ Т.А. Зеленина

Исполнители:

Студент группы 09 ММЭ Д.А. Евдокимов

Оренбург 2011

Содержание

1. Задание на практику

2. Безусловная минимизация функции многих переменных

2.1 Минимизация функции вдоль заданного направления

2.2 Градиентный метод наискорейшего спуска

2.3 Метод сопряженных градиентов

2.4 Блок-схема алгоритма минимизации

2.5 Листинг программы «Безусловная минимизация функции»

1. Задание на практику

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

2. Безусловная минимизация функции многих переменных

Метод сопряжённых градиентов является быстросходящимся методом, но он сходится только из начальных точек, достаточно близких к точке минимума.

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

2.1 Минимизация функции вдоль заданного направления

Численное решение задачи одномерной минимизации:

Имеем некоторую начальную точку и вектор направления . Требуется с заданной точностью е найти точку минимума функции в данном направлении и минимизирующее функцию значение б.

С помощью метода удвоения шага локализуем интервал, на котором находится бmin. Суть метода удвоения шага состоит в том, что мы выбираем некоторое начальное значение б (в данном случае выберем нулевое) и некоторую величину шага h. Следующую точку итерационной последовательности бi находим по формуле бi=бi-1+2i-1*h. Затем вычисляем хi=х0+бih. Итерации продолжаем до тех пор, пока не выполнится условие f(хi)>f(хi-1). При выполнении этого условия полагаем бmin=[бi-2; бi].

Затем уточним границы интервала (и, соответственно, приближённое значение точки минимума) с помощью одного из быстросходящихся методов. Например, с помощью метода золотого сечения.

Суть метода золотого сечения состоит в использовании одноимённой пропорции. Точка на отрезке является его «золотым сечением», если она делит отрезок так, что его большая часть соотносится с ним самим так же, как и меньшая часть с большей. На отрезке существует две таких точки. Обозначим их б(1) и б(2). Пусть a,b - левая и правая границы отрезка соответственно, тогда их координаты можно выразить, используя формулы . Находя х(1)=х0+б(1)h и х(2)=х0+б(2)h, сравниваем f(х(1)) и f(х(2)). Если f(х(1))>f(х(2)), то . Если f(х(1))<f(х(2)), то . Если f(х(1))=f(х(2)), то . При необходимости находим новые точки золотого сечения. Следует заметить, что точка б(1) является точкой золотого сечения для отрезка [a; б(2)], а точка б(2) - точкой золотого сечения для [б(1); b]. Абсолютный критерий останова имеет вид . Относительный критерий останова: . steptol обычно принимают равным 10-p, где p - желаемое число значащих цифр после запятой.

2.2 Градиентный метод наискорейшего спуска

Имеем оптимизационную задачу минимизации функции: , где . Каждое следующее приближение находим по формуле . Если вектор hk равен антиградиенту функции в точке хk, а бk выбирается из условия , то мы имеем метод наискорейшего спуска.

Метод наискорейшего спуска является методом первого порядка, так как использует информацию о значении первых производных функции. К группе таких методов относятся различные градиентные методы.

Градиентные методы представляют собой методы спуска, в которых в качестве направлений убывания выбирается направление антиградиента функции:-grad

В качестве абсолютного критерия останова обычно выбирают либо неравенство (критерий сильной сходимости) или ¦grad¦. В качестве относительного критерия останова обычно выбирают либо , либо .

С помощью метода минимизации функции вдоль заданного направления осуществляется поиск минимума из заданной точки. В методе наискорейшего спуска в качестве вектора направления используют вектор антиградиента.

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

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

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

Геометрическая иллюстрация метода представлена на рисунке 1.



Рисунок 1 - геометрическая иллюстрация метода наискорейшего спуска

2.3 Метод сопряженных градиентов

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

Рассмотрим следующую частную задачу оптимизации:

Здесь - симметричная положительно определённая матрица размера . Такая задача оптимизации называется квадратичной. Имеем последовательность точек хi=хi-1+бdi-1, где б выбирается из условия , а направление d - по правилу , где . Для квадратичной функции метод сопряжённых градиентов сходится за n шагов, где n - число переменных в функции.

В общем случае (неквадратичная функция) для задачи минимизации имеем , если и , если

В качестве абсолютного критерия останова обычно выбирают либо неравенство (критерий сильной сходимости) или ¦grad¦. В качестве относительного критерия останова обычно выбирают либо , либо .

Рисунок 2 - Геометрическая иллюстрация метода сопряженных градиентов

2.4 Блок-схема алгоритма минимизации



2.5 Листинг программы «Безусловная минимизация функции»

алгоритм сопряженный градиент минимизация

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, XPMan;

type

TForm1 = class(TForm)

Button1: TButton;

Label1: TLabel;

Edit2: TEdit;

XPManifest1: TXPManifest;

Label2: TLabel;

Edit3: TEdit;

Button4: TButton;

Button5: TButton;

Label3: TLabel;

Edit4: TEdit;

edt1: TEdit;

lst1: TListBox;

lbl1: TLabel;

edt2: TEdit;

lst2: TListBox;

btn1: TButton;

btn2: TButton;

lbl2: TLabel;

lst3: TListBox;

btn3: TButton;

btn4: TButton;

edt3: TEdit;

Label4: TLabel;

procedure Button4Click(Sender: TObject);

procedure Button5Click(Sender: TObject);

procedure Button1Click(Sender: TObject);

procedure btn1Click(Sender: TObject);

procedure btn3Click(Sender: TObject);

procedure btn2Click(Sender: TObject);

procedure btn4Click(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

const nvar=5;

var x,x0,x00,xleft,d:array [0..(nvar-1)] of real;

eps,h,lam0: real;

n:integer;

macht,tipich: array [0..nvar-1] of real; //Коэф-ты масштабирования и типичные значения х и ф-ции

Form1: TForm1;

implementation

{$R *.dfm}

function myfunc(x: array of real):real;

begin

myfunc:=sqr(x[0])+2*x[0]+2*sqr(x[1])+4*x[1]+5;

end;

procedure FindAlpha(x0,d:array of Real; h,eps:real; n:integer; var lam0:real);

var lam: array [1..3] of real;

i:integer;

Procedure newlam;

begin

lam[1]:=lam[2];

lam[2]:=lam[3];

lam[3]:=lam[3]+h;

h:=2*h;

end;

Procedure newx;

var i:integer;

begin

newlam;

for i:=0 to (n-1) do

x[i]:=x0[i]+d[i]*lam[3];

end;

procedure FindInterval;

var i:Integer;

begin

if myfunc(x)>myfunc(x0) then begin

h:=-h;

newx;

end;

while myfunc(x)<myfunc(x0) do begin

for i:=0 to (n-1) do begin

xleft[i]:=x0[i];

x0[i]:=x[i];

end;

newx;

end;

end;

procedure MakeItAccurate;

var i:integer;

lam1,lam2:real; //Золотое сечение

x1,x2: array [0..(nvar-1)] of real;

chosen: byte; //1 - интервал [a, lam2], 2 - интервал [lam1, b], 1 - интервал [lam1, lam2]

procedure choose(lam1,lam2:real); //выбираем новый отрезок

var i:integer;

u,v:real; //Чтобы постоянно функцию не вычислять

begin

for i:=0 to (n-1) do begin

x1[i]:=x0[i]+lam1*d[i];

x2[i]:=x0[i]+lam2*d[i];

end;

u:=myfunc(x1);

v:=myfunc(x2);

if u<v then chosen:=1;

if u>v then chosen:=2;

if u=v then chosen:=3;

end;

procedure zollam; //Считаем новые лямбда

begin

if chosen=1 then begin

lam[3]:=lam2;

lam2:=lam1;

lam1:=lam[3]-0.618*(lam[3]-lam[1]);

end;

if chosen=2 then begin

lam[1]:=lam1;

lam1:=lam2;

lam2:=lam[1]+0.618*(lam[3]-lam[1]);

end;

if chosen=3 then begin

lam[1]:=lam1;

lam[3]:=lam2;

lam1:=lam[3]-0.618*(lam[3]-lam[1]);

lam2:=lam[1]+0.618*(lam[3]-lam[1]);

end;

end;

function makestop:boolean;

var norm:real;

i:integer;

function relxnorm(x1,x2:array of real):real;

var relx:array [0..nvar-1] of real;

k:integer;

norme:real;

begin

for k:=0 to (n-1) do begin

relx[k]:=abs(x2[k]-x1[k]);

if abs(x2[k])>tipich[k] then relx[k]:=relx[k]/abs(x2[k])

else relx[k]:=relx[k]/tipich[k];

end;

norme:=0;

for k:=0 to (n-1) do

norme:=norme+sqr(relx[k]);

norme:=sqrt(norme);

relxnorm:=norme;

end;

begin

choose(lam1,lam2);

for i:=0 to (n-1) do begin

xleft[i]:=x0[i]+d[i]*lam[1];

x[i]:=x0[i]+d[i]*lam[3];

end;

if relxnorm(xleft,x)<eps then makestop:=true

else makestop:=false;

end;

begin

for i:=0 to (n-1) do x0[i]:=x00[i+1];

lam1:=lam[3]-0.618*(lam[3]-lam[1]);

lam2:=lam[1]+0.618*(lam[3]-lam[1]);

while not makestop do zollam;

end;

begin

lam[1]:=0;

lam[2]:=0;

lam[3]:=0;

newx;

FindInterval;

MakeItAccurate;

for i:=1 to n do

x[i]:=(xleft[i]+x[i])/2;

lam0:=(lam[3]+lam[1])/2;

end;

Procedure spusken(x0:array of real; eps,h:real; n:integer; var x1:array of real);

var xo,xn,agrad:array [0..(nvar-1)] of Real; //xold, xnew, antigrad.

alpha,gradnorm:Real;

j:integer;

procedure grads(x:array of Real; var agrad: array of Real);

var i:integer;

g:real;

gradtole: array [0..(nvar-1)] of Real;

begin

agrad[0]:=-(2*x[0]+2);

agrad[1]:=-(4*x[1]+4);

g:=myfunc(x);

if abs(g)>tipich[n] then g:=abs(g)

else g:=tipich[n];

for i:=0 to (n-1) do begin

gradtole[i]:=abs(-agrad[i]);

if abs(x[i])>tipich[i] then gradtole[i]:=gradtole[i]*abs(x[i])

else gradtole[i]:=gradtole[i]*tipich[i];

gradtole[i]:=gradtole[i]/g;

end;

gradnorm:=gradtole[0];

for i:=0 to (n-1) do

if gradnorm<gradtole[i] then gradnorm:=gradtole[i];

end;

procedure xnew(alpha:real; x0:array of Real; var x:array of Real);

var i:integer;

begin

for i:=0 to (n-1) do

x[i]:=x0[i]+alpha*agrad[i];

end;

function stop (norm,eps:Real):Boolean;

begin

if norm<eps then stop:=true

else stop:=false;

end;

begin

for j:=0 to (n-1) do xo[j]:=x0[j];

repeat

grads(xo,agrad);

FindAlpha(xo,agrad,h,eps,n,alpha);

xnew(alpha,xo,xn);

for j:=0 to (n-1) do xo[j]:=xn[j];

until stop(gradnorm,eps);

for j:=0 to (n-1) do x1[j]:=xn[j];

end;

procedure soprgrad(x0:array of Real; eps,h:real; n:Integer; var x:array of real);

var xold, xnew, nd, antigrad, predd: array [0..(nvar-1)] of real; //predd - предыдущее d, newd - новое d

beta,teknorm,prednorm,alpha:Real; //текушая и предыдущая нормы

k,counter:integer;

procedure grads(x:array of Real; var agrad: array of Real; gradnorm:real);

var g:real;

i:integer;

gradtole: array [0..(nvar-1)] of Real;

begin

agrad[0]:=-(2*x[0]+2);

agrad[1]:=-(4*x[1]+4);

g:=myfunc(x);

if abs(g)>tipich[n] then g:=abs(g)

else g:=tipich[n];

for i:=0 to (n-1) do begin

gradtole[i]:=abs(-agrad[i]);

if abs(x[i])>tipich[i] then gradtole[i]:=gradtole[i]*abs(x[i])

else gradtole[i]:=gradtole[i]*tipich[i];

gradtole[i]:=gradtole[i]/g;

end;

gradnorm:=gradtole[0];

for i:=1 to (n-1) do

if gradnorm<gradtole[0] then gradnorm:=gradtole[i];

end;

function newbeta(ti:Integer; norm,prednorm:real):Real; //ti - текущая итерация

begin

if ((ti mod n)=0) then newbeta:=0

else newbeta:=Sqr(norm/prednorm);

end;

procedure newd(agrad,d:array of Real; beta:Real; var nd:array of real);

var i:integer;

begin

for i:=0 to (n-1) do

nd[i]:=agrad[i]+beta*d[i];

end;

function theend(norm,eps:real):boolean;

begin

if norm<(eps/1000) then theend:=true

else theend:=false;

end;

procedure findx(alpha:real; d,x0:array of Real; var x:array of Real);

var i:integer;

begin

for i:=0 to (n-1) do

x[i]:=x0[i]+alpha*d[i];

end;

begin

for k:=0 to (n-1) do xold[k]:=x0[k];

spusken(xold,eps,h,n,xnew);

grads(xnew,antigrad,prednorm);

for k:=0 to (n-1) do xold[k]:=xnew[k];

counter:=0;

repeat

beta:=newbeta(counter,teknorm,prednorm);

newd(antigrad,predd,beta,nd);

FindAlpha(xold,nd,h,eps,n,alpha);

findx(alpha,nd,xold,xnew);

for k:=0 to (n-1) do xold[k]:=xnew[k];

prednorm:=teknorm;

counter:=counter+1;

until theend(teknorm,eps);

for k:=0 to (n-1) do x[k]:=xnew[k];

end;

procedure TForm1.Button4Click(Sender: TObject);

var g:real;

begin

Lst1.AddItem(Edit3.Text, Edit3);

try

g:=StrToFloat(Edit3.Text);

except

ShowMessage('Enter real number!');

Lst1.Items.Delete(Lst1.Items.Count-1);

end;

end;

procedure TForm1.Button5Click(Sender: TObject);

begin

Lst1.Items.Delete(Lst1.ItemIndex);

end;

procedure TForm1.Button1Click(Sender: TObject);

var i:integer;

g:real;

begin

try

eps:=StrToFloat(Edit4.Text);

h:=StrToFloat(Edt1.Text);

except

ShowMessage('Enter real number!');

Edit4.Text:='';

edt1.Text:='';

end;

n:=Lst1.Items.Count;

for i:=0 to (n-1) do begin

x0[i]:=StrToFloat(Lst1.Items[i]);

macht[i]:=StrToFloat(Lst2.Items[i]);

tipich[i]:=StrToFloat(Lst3.Items[i]);

end;

g:=myfunc(tipich);

tipich[n]:=g;

soprgrad(x0,eps,h,n,x);

Edit2.Text:='';

for i:=0 to (n-1) do

Edit2.Text:=Edit2.Text+FloatToStr(x[i])+'; ';

end;

procedure TForm1.btn1Click(Sender: TObject);

var g:real;

begin

Lst2.AddItem(Edt2.Text, Edt2);

try

g:=StrToFloat(Edt2.Text);

except

ShowMessage('Enter real number!');

Lst2.Items.Delete(Lst2.Items.Count-1);

end;

end;

procedure TForm1.btn3Click(Sender: TObject);

var g:real;

begin

Lst3.AddItem(Edt3.Text, Edit3);

try

g:=StrToFloat(Edt3.Text);

except

ShowMessage('Enter real number!');

Lst3.Items.Delete(Lst1.Items.Count-1);

end;

end;

procedure TForm1.btn2Click(Sender: TObject);

begin

Lst2.Items.Delete(Lst2.ItemIndex);

end;

procedure TForm1.btn4Click(Sender: TObject);

begin

Lst3.Items.Delete(Lst3.ItemIndex);

end;

end.

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