Алгоритм метода Ньютона
1. Задать размерность задачи оптимизации п, координаты начальной точки , точность поиска .
2. Положить счетчик числа итераций .
3. Определить направление вектора градиента целевой
функции
в точке . Для вычисления координат вектора градиента использовать разностную формулу (2.3) .
4. Проверить условие окончания поиска
Если условие выполнено, то расчет окончен , иначе перейти к пункту 5.
5. Сформировать матрицу Гессе , используя разностные формулы вычисления вторых (2.5) и смешанных производных (2.6).
6. Проверить положительную определенность матрицы
Гессе . Если матрица положительно определена , то перейти к пункту 7, иначе — к пункту 8.
7. Определить координаты точки и перейти к пункту 10.
8. Вычислить шаг по формуле (2.4), используя резуль
таты вычислений пункта 3 и разностные формулы (2.5), (2.6).
9. Определить координаты точки по методу наискорейшего градиентного спуска.
10. Положить и перейти к пункту 3.
Программы оптимизации
Для подтверждения работоспособности программы минимизации функции, из учебника [6], решим пример № 2.22. В данном примере, необходимо найти минимум целевой функции, методом Ньютона, с точностью
Для разработки программы была использована среда PascalABC.
Текст программы
Programlol;
labelstop;
constdx=0.0001;
typearray2D=array[1..4,1..4] ofreal;
varx,df:array ofreal;
d2f:array2D;
eps,h,grad,s:real;
i,j,n,q:integer;
functionf(x:array ofreal):real;
Begin
f:=2*sqr(x[1])+4*sqr(x[2])+8*sqr(x[3])+2*x[1]*x[2]-x[1]*x[3]+2*x[2]*x[3]+6*x[1]-7*x[3]
End;
Procedureinvert(n,q:integer;matr1:array2D;varmatr:array2D);
labelM1,M2;
vara:array[1..4,1..8] ofreal;
i,j,k,m:integer;
t:real;
Begin
m:=2*n; q:=0;
fori:=1 ton do
forj:=1 tom do
ifj<=n thena[i,j]:=matr1[i,j] else
ifj=n+i thena[i,j]:=1.0 elsea[i,j]:=0;
fori:=1 ton do
begink:=i;
M1:ifa[k,i]=0 then
beginq:=1;
ifk<n thenk:=k+1 else gotoM2;
gotoM1;
end;
ifq=1 then
forj:=1 tom do
begint:=a[k,j]; a[k,j]:=a[i,j]; a[i,j]:=t
end;
forj:=m downtoi doa[i,j]:=a[i,j]/a[i,i];
fork:=1 ton do
ifk<>i then
forj:=m downto1 do
a[k,j]:=a[k,j]-a[i,j]*a[k,i];
end;
q:=0;
fori:=1 ton do
forj:=1 ton domatr[i,j]:=a[i,j+n];
M2: end;
{Процедура вычисления координат вектора градиента и фомирования матрицы вторых производных}
procedureFor_Mat_d2f(n:integer; vargrad:real);
Var
i,j:integer;
s,f0:real;
Begin
//расчет производных
f0:=f(x);s:=0;
fori:=1 ton do
Begin
x[i]:=x[i]+dx;
df[i]:=(f(x)-f0)/dx;
s:=sqr(df[i]);
x[i]:=x[i]-dx;
end;
grad:=sqrt(s);
//расчет вторых производных
fori:=1 ton do
Begin
s:=-2*f(x);
x[i]:=x[i]+dx;
s:=s+f(x);
x[i]:=x[i]-2*dx;
s:=s+f(x);
x[i]:=x[i]+dx;
d2f[i,i]:=s/sqr(dx);
end;
//расчет смешанных производных
fori:=1 ton-1 do
forj:=i+1 ton do
Begin
s:=f(x); {1}
x[i]:=x[i]-dx; x[j]:=x[j]-dx;
s:=s+f(x); {4}
x[j]:=x[j]+dx;
s:=s-f(x);{2}
x[i]:=x[i]+dx; x[j]:=x[j]-dx;
s:=s-f(x);{3}
x[j]:=x[j]+dx;
d2f[i,j]:=s/sqr(dx);
d2f[j,i]:=d2f[i,j];
end;
end; {For_Mat_d2f}
//*****************
Begin
writeln('исходныеданные');
writeln;
write('введите размерность задачи оптимизации n=');
readln(n);
SetLength(x,n+1);
SetLength(df,n+1);
write('введите точность вычислений eps=');
readln(eps);
writeln('введите начальные значения переменных');
fori:=1 ton do
Begin
write('x[',i,']=');
readln(x[i]);
end;
Repeat
For_Mat_d2f(n,grad);
invert(n,q,d2f,d2f);
ifq=1 then beginwriteln('определительравеннулю');
gotostop end;
fori:=1 ton do
Begin
s:=0;
forj:=1 ton do
s:=s+d2f[i,j]*df[j];
x[i]:=x[i]-s;
end;
untilgrad<eps;
writeln;
writeln('результатыоптимизации:'); writeln;
fori:=1 ton do
writeln('x[',i,']=',x[i]);
writeln;
writeln('значение функции цели = ',f(x));
stop:readln;
end.
Результаты расчета
Полученные значения совпадают с результатами минимизации из учебника [6].