Методы оптимизации
Методы многомерной оптимизации
Работы ЛЮБОЙ сложности в MATLAB, для заказа -- контакты в шапке сайта.
Цель работы
Изучение среды МАТЛАБ, создание программы, реализующей алгоритмы наискорейшего спуска и Ньютона, для поиска минимума заданных сферической (f1) и эллипсоидной (f2) функций.
1. Входные данные:
В качестве области поиска выбирается [-5; 5]D,
где D=3 – размерность пространства.
2 Теоретический материал
2.1. Метод наискорейшего спуска
Новое приближение к точке оптимума x(i,j,k) находится из простого соотношения:
(3)
где частные производные находятся по координатам точки m-ного приближения к искомой точке оптимума;
α=[a,b] – шаговый множитель, который находится из условия:
;
- функция, зависящая в данной точке приближения только от альфа..
Для этого. для каждой функции и каждого метода вводится своя область поиска значений α ∊ [a,b] и число дискретных точек n в этой области, для каждой из которых находятся значения F(α), соответственно выбранной для анализа функции (сферическая или эллипсоидная). После, из найденных значений, выбирается min{F(α)} и соответствующий ему α.
После этого окончательно вычисляется .
В качестве критерия окончания поиска применяется требование малости изменения значений оптимизируемой функции при перемещении к каждому новому приближению:
||f(xm+1) - f(xm)||<Δ.
Значение Δ выбирается меньшим (Δ =0,0001) для сферической функции и бОльщим – для эллипсоидной (Δ =0,1), поскольку значения эллипсоидной функции значительно превышают значения сферической в одной и той же области трёхмерного пространства.
2.1. Метод Ньютона
Метод Ньютона запишем в виде:
Новое приближение к точке оптимума x(i,j,k) находится из соотношения:
(4)
где частные производные находятся по координатам точки m-ного приближения к искомой точке оптимума;
α=[0;1] – корректирующий множитель.
Принцип его выбора практически без изменений применяется согласно алгоритму, описанному в п.п. 2.1.
Задача сводится к определению обратной матрицы Гессе, которая, в данном случае, вычислена точно, поскольку можно точно определить значения производных исследуемых функций в каждой точке области поиска.
В качестве критерия окончания поиска применяется требование малости изменения значений оптимизируемой функции при перемещении к каждому новому приближению:
||f(xm+1) - f(xm)||<Δ.
Значение Δ выбирается меньшим (Δ =0,0001) для сферической функции и бОльщим – для эллипсоидной (Δ =0,1), поскольку значения эллипсоидной функции значительно превышают значения сферической в одной и той же области трёхмерного пространства.
3. Некоторые особенности алгоритма
В работе программы используется сетка значений координат с постоянными шагами {h1, h2, h3}, которые, в обшем случае, могут быть не равными, что особенно актуально при анализе подобных заданной в работе эллипсоидной (1) функций.
После применения соотношений (3) метода спуска и (4) метода Ньютона, полученные координаты, конечно же, не совпадают с узлами сетки. Поэтому, для дальнейшего использования, производится коррекция найденного приближения координат:
- либо подискивается точка ближайшая по координатамЮ,
- либо из ближайшего окружения данной точки определяется та, в которой функция f(x) даёт меньшее значение, чем только, что вычисленное приближение к оптимуму f(xm+1).
Фактически, используются оба подхода, что, в случае применения метода наискорейшего спуска к эллипсоидной функции, позволило серьёзно увеличить значения шагового параметра и повысить скорость сходимости с одновременным обеспечением точности даже при относительно крупном шаге h=0.1 и числе узловых точек по каждой координате ≤201.
Каждая найденная точка отмечается на трёхмерном графике перекрещивающимися пунктирными линиями, показывающими привязку к кординатной сетке,а путь от точки предыдущейитерации к только что найденной точке отмечается линией синего цвета.
4. Текст программы
% Поиск минимума f1(x) и f2(x)
% Методами НаиСкорейшего спуска и
clear all; clc; close all;
% Число шагов при повороте изображения:
%Vstep=20; dn(1:3)=0;
% Тип поиска:
s='min'; znak=eye(3); dd(2)=0; dd(1)=1e8;
% Задаём координатную сетку в окрестностях {-5;5}^3 точки M0(0;0;0):
% с равными шагами по координатам:
x0(1)=0;x0(2)=0;x0(3)=0;
dx(1)=5;dx(2)=5;dx(3)=5;
h(1)=0.1;h(2)=0.1;h(3)=0.1;
for i=1:3
x(i,:)=x0(i)-dx(i):h(i):x0(i)+dx(i);
end;
% Находим максимальные значения индексов:
imax=length(x(1,:));
jmax=length(x(2,:));
kmax=length(x(3,:));
% Находим индексы середины отрезка:
i12=(imax-1)/2+1; j12=(jmax-1)/2+1; k12=(kmax-1)/2+1;
% Формируем поле значений функции:
% Начальное значение пункта меню функций
mf=1;
% Цикл по выбору функции:
while mf<3
% Счётчик общего числа произведённых итераций
% (вычислений функции f(x1,x2,x3):
summ=0;
% Счётчик числа испытаний
% (для данного метода и данной фунцкии)
NN=0;
% Вызываем меню выбора функции:
mf=menu('Select function:','1. Sphere Function','2. Ellipsoidal Function','exit program');
switch mf
case 1
% Ограничение по числу шагов поиска
N=100;
% Точность вычисления min f1(x1,x2,x3)
d=0.0001;
% Определяем сферическую функцию
for ii=1:imax
for jj=1:jmax
for kk=1:kmax
% Определяем сферическую функцию
f(ii,jj,kk)=(x(1,ii)-0)^2+(x(2,jj)-0)^2+(x(3,kk)-0)^2+0;
end;
end;
end;
case 2
% Ограничение по числу шагов поиска
N=200;
% Точность вычисления min f2(x1,x2,x3)
d=0.1;
% Определяем эллипсоидную функцию
for ii=1:imax
for jj=1:jmax
for kk=1:kmax
% Определяем эллипсоидную функцию
f(ii,jj,kk)=10^(6*(1-1)/(3-1))*(x(1,ii)-0)^2+10^(6*(2-1)/(3-1))*(x(2,jj)-0)^2+10^(6*(3-1)/(3-1))*(x(3,kk)-0)^2+0;
end;
end;
end;
case 3
% Выходим из программы
break;
end;
% Начальное значение пункта меню методов
mm=1;
% Цикл по выбору метода:
while mm<3
% Вызываем меню выбора метода:
mm=menu('Select Method:','1. Steepest descent','2. Newton''s method','Go to seLect f1 / f2');
switch mm
case 3
% Выходим в цикл выбора функций
break;
end;
close all;
Az0=-147;
Et0=22;
% Открываем графическое окно:
figure(1);
% Режим сохранения нарисованных изображений
hold on;
% Включаем координатную сетку
grid on;
% Выбираем ракурс обзора
view(Az0,Et0);
% Обозначения осей:
xlabel('x1');
ylabel('x2');
zlabel('x3');
axis([-dx(1) dx(1) -dx(2) dx(2) -dx(3) dx(3)]);
% Устанавливаем счётчик точек траектории в начальное значение
m=1;
% Начальную точку задаём индексами элементов
% массивов x1, x2, x3 так, что бы она оказалась
% на некотором расстоянии от возможного min или max:
%i(m)=imax; j(m)=jmax; k(m)=kmax;
i(m)=fix(rand*imax)+1; j(m)=fix(rand*jmax)+1; k(m)=fix(rand*kmax)+1;
text(x(1,i(m))-h(1),x(2,j(m))-h(2),x(3,k(m)),strcat(int2str(m)));
% Начинаем цикл итераций
for nn=1:N
% При mm=1 Расчёт методом наискорейшего спуска
% При mm=2 Расчёт методом Ньютона
% Первые Частные производные
% рассчитываются одинаково для обоих методов
% Вычисляем частную производную функции f(x1,x2,x3) по x1 через конечные разности:
switch i(m)
case 1
fx(1)=(f(i(m)+1,j(m),k(m))-f(i(m),j(m),k(m)))/(x(1,i(m)+1)-x(1,i(m)));
case imax
fx(1)=(f(i(m),j(m),k(m))-f(i(m)-1,j(m),k(m)))/(x(1,i(m))-x(1,i(m)-1));
otherwise
fx(1)=(f(i(m)+1,j(m),k(m))-f(i(m)-1,j(m),k(m)))/(x(1,i(m)+1)-x(1,i(m)-1))/2;
end;
% Вычисляем частную производную функции f(x1,x2,x3) по x2 через конечные разности:
switch j(m)
case 1
fx(2)=(f(i(m),j(m)+1,k(m))-f(i(m),j(m),k(m)))/(x(2,j(m)+1)-x(2,j(m)));
case jmax
fx(2)=(f(i(m),j(m),k(m))-f(i(m),j(m)-1,k(m)))/(x(2,j(m))-x(2,j(m)-1));
otherwise
fx(2)=(f(i(m),j(m)+1,k(m))-f(i(m),j(m)-1,k(m)))/(x(2,j(m)+1)-x(2,j(m)-1))/2;
end;
% Вычисляем частную производную функции f(x1,x2,x3) по x2 через конечные разности:
switch k(m)
case 1
fx(3)=(f(i(m),j(m),k(m)+1)-f(i(m),j(m),k(m)))/(x(3,k(m)+1)-x(3,k(m)));
case jmax
fx(3)=(f(i(m),j(m),k(m))-f(i(m),j(m)-1,k(m)-1))/(x(3,k(m))-x(3,k(m)-1));
otherwise
fx(3)=(f(i(m),j(m),k(m)+1)-f(i(m),j(m),k(m)-1))/(x(3,k(m)+1)-x(3,k(m)-1))/2;
end;
% Определяем шаговый параметр, корректируем его знак, если нужно:
switch mf
case 1
% Для сферической функции f1:
% Задаём ДИАПАЗОН параметра шага:
switch mm
case 1
alfa=0.1:0.05:1;
znak=eye(3);
case 2
alfa=0.1:0.1:1;
% Рассчитываем обратную матрицу Гессе
% для сферической функции для метода Ньютона
znak=inv(eye(3)*2);
end;
F=(x(1,i(m))-(znak(1,:)*fx')*alfa).^2+(x(2,j(m))-(znak(2,:)*fx')*alfa).^2+(x(3,k(m))-(znak(3,:)*fx')*alfa).^2;
case 2
% Для эллипсоидной функции f2:
% Задаём ДИАПАЗОН параметра шага:
switch mm
case 1
alfa=5e-8:5e-8:5e-7;
if dd(2)>dd(1)
znak=-eye(3);
else
znak=eye(3);
end;
case 2
alfa=0.1:0.1:1;
% Рассчитываем обратную матрицу Гессе
% для эллипсоидной функции для метода Ньютона
znak=inv(2*[10^(6*(1-1)/(3-1)) 0 0;
0 10^(6*(2-1)/(3-1)) 0;
0 0 10^(6*(3-1)/(3-1))]);
end;
F=10^(6*(1-1)/(3-1))*(x(1,i(m))-(znak(1,:)*fx')*alfa).^2+10^(6*(2-1)/(3-1))*(x(2,j(m))-(znak(2,:)*fx')*alfa).^2+10^(6*(3-1)/(3-1))*(x(3,k(m))-(znak(3,:)*fx')*alfa).^2;
end;
[Fmin,n]=min(F);
% Находим координаты следующей точки:
x11=x(1,i(m))-(znak(1,:)*fx')*alfa(n);
x21=x(2,j(m))-(znak(2,:)*fx')*alfa(n);
x31=x(3,k(m))-(znak(3,:)*fx')*alfa(n);
% Увеличиваем счётчик точек траектории
m=m+1;
% Находим индексы ближайшей по координатам точки в координатной сетке:
[x11,i(m)]=min(abs(x(1,:)-x11));
[x21,j(m)]=min(abs(x(2,:)-x21));
[x31,k(m)]=min(abs(x(3,:)-x31));
% Корерктируемся в соседнюю точку,
% где значение функции меньше:
Fmin=1e16;
for ii=i(m)-1:i(m)+1
for jj=j(m)-1:j(m)+1
for kk=k(m)-1:k(m)+1
if and(and(and(ii>0,ii<imax+1),and(jj>0,jj<jmax+1)),and(kk>0,kk<kmax+1))
if f(ii,jj,kk)<Fmin
iii=ii; jjj=jj; kkk=kk;
Fmin=f(ii,jj,kk);
end;
end;
end;
end;
end;
i(m)=iii; j(m)=jjj; k(m)=kkk;
%end;
% Рисуем точку минимума на m-ном шаге приближения
%plot3([0 x(1,i(m))],[x(2,j(m)) x(2,j(m))],[x(3,k(m)) x(3,k(m))],'k--');
%plot3([x(1,i(m)) x(1,i(m))],[0 x(2,j(m))],[x(3,k(m)) x(3,k(m))],'k--');
%plot3([x(1,i(m)) x(1,i(m))],[x(2,j(m)) x(2,j(m))],[0 x(3,k(m))],'k--');
plot3([-5 x(1,i(m))],[x(2,j(m)) x(2,j(m))],[x(3,k(m)) x(3,k(m))],'k--');
plot3([x(1,i(m)) x(1,i(m))],[-5 x(2,j(m))],[x(3,k(m)) x(3,k(m))],'k--');
plot3([x(1,i(m)) x(1,i(m))],[x(2,j(m)) x(2,j(m))],[-5 x(3,k(m))],'k--');
plot3([x(1,i(m)) 5],[x(2,j(m)) x(2,j(m))],[x(3,k(m)) x(3,k(m))],'k--');
plot3([x(1,i(m)) x(1,i(m))],[x(2,j(m)) 5],[x(3,k(m)) x(3,k(m))],'k--');
plot3([x(1,i(m)) x(1,i(m))],[x(2,j(m)) x(2,j(m))],[x(3,k(m)) 5],'k--');
% Рисуем путь от предыдущей m-ой к (m+1)-ой точке
plot3([x(1,i(m-1)) x(1,i(m))],[x(2,j(m-1)) x(2,j(m))],[x(3,k(m-1)) x(3,k(m))],'ob-','LineWidth',2);
% Подписываем номер точки
%text(x(1,i(m))+h(1),x(2,j(m))+h(2),x(3,k(m))+h(3),strcat(int2str(m)));
%text(x(1,i(m))-10*h(1),x(2,j(m))-10*h(2),x(3,k(m))-10*h(3),strcat(int2str(m)));
text(x(1,i(m)),x(2,j(m)),x(3,k(m))-10*h(3),strcat(int2str(m)));
% Вычисляем полученное приращение функции и запоминаем предыдующее:
dd(1)=dd(2);
dd(2)=abs(f(i(m),j(m),k(m))-f(i(m-1),j(m-1),k(m-1)));
% Вывод номера найденной точки приближения в заголовок графика:
title(strcat('Point number m=',int2str(m),' f',s,'=f(',num2str(x(1,i(m))),';',num2str(x(2,j(m))),';',num2str(x(3,k(m))),')=',num2str(f(i(m),j(m),k(m))),' d=',num2str(dd(2)),' alfa=',num2str(alfa(n))));
% Делаем паузу для отображения
pause(0.5);
% Проверка условия окончания счёта
%if sqrt((x(1,i(m))-x(1,i(m-1)))^2+(x(2,j(m))-x(2,j(m-1)))^2+(x(3,k(m))-x(3,k(m-1)))^2)<d
%if abs(f(i(m),j(m),k(m))-f(i(m-1),j(m-1),k(m-1)))<d
if dd(2)<d
% Прерывание цикла итераций, если условие выполняется
break;
end;
% конец цикла итераций
end;
disp('Число испытаний (для данного метода и данной фунцкии) NN=');
NN=NN+1
disp('Общее число произведённых итераций (вычислений функции f(x1,x2,x3): summ=');
summ=summ+m
% Выключаем запоминание выводимых изображений
hold off;
% конец цикла выбора метода по mm
end;
% конец цикла выбора функций по mf
end;
5. Результаты работы программы
Для каждого случая выбранной для анализа функции и используемого метода, рассчитаем ряд показателей:
N – общее число запусков метода;
Ns – число реализаций метода, приведших к попаданию в точку истинного минимума;
- среднее число вычислений целевой функции, необходимых для достижения минимума (считается по успешно завершившимся запускам);
– относительное ожидаемое время исполнения (как отношение суммы всех вычислений функции по всем и успешным, и неуспешным запускам к числу успешных запусков);
Показатели испытаний:
5.1. Метод наискорейшего спуска для сферической функции
5.2. Испытания метода Ньютона для сферической функции
5.3. Испытания метода наискорейшего спуска для эллипсоидной функции
5.4. Испытания метода Ньютона для эллипсоидной функции
Выводы
В процессе выполнения работы была разработана программа реализующая методы скорейшего спуска и Ньютона для тестирования этих алгоритмов на сферической и эллипсоидной функции в трёхмерном пространстве.
Оба метода показали себя быстро и точно сходящимися для сферической функции.
На эллипсоидной функции метод Ньютона показал ещё более лучшие результаты.
Метод наискорейшего спуска для эллипсоидной функции в некоторых испытаниях «промахивался» мимо точки минимума, а, в большинстве экспериментов, показал заметно большее число итераций (вычислений оптимизируемой функции) до сходимости.
В случае сходимости, ряд экспериментов с методом наискорейшего спуска для эллипсоидной функции закончился определением минимума функции с погрешностью, но в пределах заданной. Сама величина заданной абсолютной погрешности минимального значения, используемая для реализации КОП (критерия окончания поиска) была установлена несколько выше, чем для сферической функции, поскольку значения эллипсоидной функции на несколько порядков превышают значению сферической функции в заданной области поиска.
Что бы возвращать алгоритм на траекторию спуска была введена коррекция значений шагового параметра «альфа» по знаку и коррекция координат при их последовательном отклонении в сторону больших значений минимизируемой функции.
Данный механизм так же положительно показал себя в остальных трёх случаях, позволяя добиться абсолютной точности в определении как координат точки минимума, так и минимального значения функции.
Список использованных источников
1. Зазарыкин В.М. Численные методы : учеб. пособие / В.М. Зазарыкин, В. Г.Житомирский, М.П. Лапчик. —М. : Просвещение, 1990. — 176 с.
2. Демидович Б.П. Основы вычислительной математики / Б.П.Демидович, И.А.Марон. —М. : Наука, 1966. — 664 с.
3. МеркуловаН.Н. Методы приближенных вычислений : в 2 ч. : учеб. пособие / Н.Н.Меркулова, М.Д.Михайлов. — Томск : Томск. гос. ун-т, 2005. — Ч. 1. — 257 с.
4. Формалеев В.Ф. Численные методы / В.Ф.Формалеев, Д. Л. Ревизников. — М. : Физмалит, 2004. — 398 с.
5. Hansen, N., Finck, S., Ros, R., & Auger, A. (2009). Real-parameter black-box optimization benchmarking 2009: Noiseless functions definitions.
