Методические указания к лабораторной работе по дисциплине "Численные методы оптимизации сложных систем". Тема: Безусловная минимизация сложной функции многих переменных методами второго порядка в среде MATLAB

В настоящее время разработано огромное количество алгоритмов оптимизации сложных функций многих переменных, базирующихся на различных классических и эвристических подходах и реализованных в разнообразных системах компьютерной математики (СКМ).

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

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

Кроме того, немаловажным преимуществом СКМ Matlab при online-обучении прикладным математическим дисциплинам является удобство создания программ с возможностью перевода на другие языки программирования и наглядность представления результатов вычислений.

Целью работы является исследование возможностей СКМ Matlab для решения практической задачи нахождения численного локального минимума целевой функции двух переменных методами второго порядка - на основе алгоритмов Ньютона - Рафсона и Марквардта (Левенберга). 

                                                         Теоретическая часть

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

Примером ЦФ, для которой практически "не работают" градиентные методы, является функция Розенброка (см. Рисунок 1), выражение для которой имеет вид

оптимальная точка .

                             

                                           Рисунок 1 - График функции Розенброка

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

К численным методам второго порядка относятся методы Ньютона, Ньютона-Рафсона и квазиньютоновские алгоритмы, например, метод Марквардта (Левенберга).

Оптимизационные методы второго порядка используют квадратичную аппроксимацию ЦФ.

В дальнейшем под оптимизацией будем понимать минимизацию функций.

 

Основное выражение для выполнения итераций методом Ньютона имеет следующий вид:

                                 

где обратная матрица Гессе (гессиан или матрица вторых частных производных), вычисленная в точке  .

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

Для решения проблем, связанных с вырожденностью или отрицательной определенностью гессиана ЦФ, а также с влиянием начального приближения на скорость сходимости метода Ньютона, была создана усовершенствованная модификация рассмотренного алгоритма – метод Ньютона-Рафсона (см. (2), Рисунок 2). Главным преимуществом данного метода является использование градиентного алгоритма (в частности, метода наискорейшего спуска) в случаях возникновения указанных выше проблем на каком-либо шаге алгоритма Ньютона.

Базовая расчетная формула для метода Ньютона-Рафсона приведена в выражении (2), где в вышеуказанных вырожденных случаях применяется "метод наискорейшего спуска" Коши.

 

Блок-схема алгоритма Ньютона-Рафсона показана на Рисунке 2.

                         

                                    Рисунок 2 - Блок-схема метода Ньютона-Рафсона

Основным достоинством метода Ньютона-Рафсона является быстрая сходимость в окрестности точки экстремума, а главными недостатками – необходимость, помимо градиента, на каждой итерации вычислять обратную матрицу Гессе и зависимость сходимости метода от выбора начальной точки (эта зависимость проявляется, прежде всего, в том, что улучшение значения целевой функции с каждым шагом гарантируется только достаточным условием положительной определенности матрицы Гессе  ).

Квазиньютоновский алгоритм Марквардта (Левенберга) был разработан для решения указанных выше проблем, возникающих при использовании классического метода Ньютона.

Итерации данного алгоритма строятся согласно схеме

       

где  Е – единичная матрица;  – последовательность неотрицательных чисел, гарантирующих положительную определенность матрицы ; начальная точка  подлежит выбору.

Обычно начальное значение задается как минимум на порядок больше максимального по абсолютной величине элемента h_ij матрицы Гессе, т.е. , так что вступает в силу «закон больших чисел», который гарантирует положительную определенность специальной матрицы:

                                           

Таким образом, большим значениям  соответствует направление поиска, совпадающее с направлением антиградиента ЦФ  в точке , а при уменьшении  до нуля – направление поиска становится ньютоновским.

Если после текущей итерации получена точка с меньшим (при минимизации ЦФ) значением целевой функции (шаг результативен), то параметр  сокращают, а градиент и гессиан в новой точке пересчитывают, иначе – увеличивают значение и выполняют шаг алгоритма в направлении антиградиента, соответствующего прошлому результативному направлению поиска минимума ЦФ (градиент и гессиан не переопределяют).

Блок-схема алгоритма Марквардта (Левенберга) показана на Рисунке 3.

                   

                               Рисунок 3 - Блок-схема Марквардта (Левенберга) 

Достоинства метода Марквардта (Левенберга):

1. простота (по сравнению с методом Ньютона) реализации алгоритма;
2. высокая скорость сходимости в окрестности точки минимума ЦФ ;
3. отсутствие затруднений с выбором шага поиска минимума (данный недостаток присутствует у метода наискорейшего спуска);
4. данный метод позволяет найти локальный минимум функции и при невыполнении достаточного условия оптимальности – положительной определенности матрицы Гессе.

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

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

 

                                                              Практическая часть

Постановка задачи

Для заданной в варианте задания функции двух переменных из начальной точки выполнить поиск численного значения безусловного локального минимума ЦФ с заданной точностью вычислений с использованием алгоритмов наискорейшего спуска Коши, Ньютона-Рафсона и методом Марквардта (Левенберга) и сравнить результаты, полученные с применением стандартных средств СКМ Matlab.

Задать предельное количество итераций M = 100.

Начальные значения коэффициентов алгоритма Марквардта: С=2, . 

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

 

Порядок выполнения работы

1. Изучить теоретический материал.
2. Составить программы, реализующие оптимизационные алгоритмы Марквардта, метода наискорейшего спуска [https://hub.exponenta.ru/post/laboratornaya-rabota-po-distsipline-chislennye-metody-optimizatsii-slozhnykh-sistem-tema-bezuslovnaya-minimizatsiya-slozhnoy-funktsii-mnogikh-peremennykh-gradientnymi-metodami-v-srede-matlab3] и Ньютона-Рафсона.
3. Произвести запуск и отладку программ.
4. Выполнить численные расчеты всеми разработанными алгоритмами и сравнить результаты безусловной минимизации заданной целевой функции с использованием стандартных средств СКМ Matlab. Расхождения оценить при помощи вычислений абсолютных погрешностей по аргументам и значениям целевых функций, полученных при запуске каждого из алгоритмов и стандартных оптимизационных методов СКМ Matlab (функций fminsearch, fminunc).
5. Построить графики зависимостей количества реализаций алгоритма Марквардта и абсолютной погрешности решения задачи от принятых значений параметра С указанного алгоритма (см. Рисунок 3).
6. Оформить отчет по лабораторной работе и защитить ее.

Содержание отчета по работе

1. Титульный лист лабораторной работы.
2. Блок-схемы оптимизационных алгоритмов.
3. Тексты программ оптимизационных алгоритмов, реализованных в СКМ Matlab.
4. Результаты проведенных расчетов в виде Workspace - скрипта и графиков.
5. Выводы по работе.
6. Список использованной литературы и источников.

Пример выполнения лабораторной работы

Выполним решение задачи поиска локального минимума функции Розенброка (см. Рисунок 1).

Составим для заданной ЦФ программу в редакторе Editor СКМ Matlab с переменным набором выходных значений, в составе которых могут быть значение ЦФ, ее градиент (вектор первых частных производных) и матрица вторых частных производных (гессиан) функции:

function [f,g,H] = Rozenbrock(x)
  f = 100*(x(2) - x(1)^2)^2 +(1 - x(1))^2;
  if nargout > 1
      g = [-400*(x(2)-x(1)^2)*x(1)-2+2*x(1); 200*x(2)-200*x(1)^2];
  end
  if nargout > 2
      H = [1200*x(1)^2-400*x(2)+2 -400*x(1); -400*x(1) 200];
  end    

Тексты программ численных алгоритмов Марквардта и Ньютона-Рафсона прикреплены в приложении к Методическим указаниям.

Напишем скрипт вызова всех разработанных методов:

clear
X0 =[-2.1 2];
eps0 = 1e-5;
[XX, FF, iter_count, modgr] = SteepestDescent(@Rozenbrock, X0, eps0,100,0.01);
[XX2, FF2, iter_count2, modgr2] = Newton(@Rozenbrock, X0, eps0,50);
[~,~,H0] = Rozenbrock(X0);
mu0 = 10*max(max(abs(H0)));
[XX3, FF3, iter_count3, modgr3] = Marquardt(@Rozenbrock, X0, eps0, 5, mu0, 100);

Реализуем вызов стандартных оптимизационных методов СКМ Matlab и расчет абсолютных разностей полученных значений с запрограммированными алгоритмами:

[x,y] = fminsearch(@Rozenbrock,X0);
options = optimset('MaxFunEvals',3000, 'TolX', eps0);
[x1,y1] = fminunc(@Rozenbrock, X0, options);

delY11 = abs(FF-y); delY12 = abs(FF-y1);
delY21 = abs(FF2-y); delY22 = abs(FF2-y1);
delY31 = abs(FF3-y); delY32 = abs(FF3-y1);

После запуска скрипта на выполнение были получены следующие результаты в Workspace:

FF = 2.0880036638666111;

FF2 = 1.0608050361346755E-20;

FF3 = 5.0568688896959729E-14;

H0 = ...
  [4494 840;
   840 200];

X0 = [-2.1 2];

XX = [-0.44420212815939619 0.19253653820844233];

XX2 = [0.99999999989700461 0.99999999979400933];

XX3 = [0.99999981729272613 0.99999962147566457];

delY11 = 2.0880036637277173;

delY12 = 2.0880036638419273;

delY21 = 1.3889366658754818E-10;

delY22 = 2.4683776302227065E-11;

delY31 = 1.3884309790925927E-10;

delY32 = 2.4633207623938158E-11;

eps0 = 1E-5;

iter_count = 100;

iter_count2 = 5;

iter_count3 = 39;

modgr = 3.8578172312943875;

modgr2 = 2.0603519087026568E-10;

modgr3 = 5.5384641175017473E-6;

mu0 = 44940;

x = [0.999988584770166 0.99997687664855084];

x1 = [0.99999503641892407 0.999990051266371];

y = 1.3889366659815623E-10;

y1 = 2.4683776312835117E-11;

Примечание

В ряде случаев представленные численные алгоритмы и стандартные процедуры Matlab могут сойтись к разным локальным минимумам   ЦФ с близкими величинами  (в случае мультимодальной "овражной" целевой функции), в зависимости от заданных начальных условий поиска.

Вывод:

в ходе лабораторной работы проведена безусловная минимизация функции Розенброка 2-х переменных методами первого, второго порядка и стандартными процедурами СКМ Matlab. 

Из показанного в Workspace-скрипта видно, что градиентный метод наискорейшего спуска задачу минимизации ЦФ Розенброка не решает (modgr = 3.8578172312943875 > eps0).

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

Быстрее всего к оптимальному решению сошелся метод Ньютона-Рафсона (за 5 итераций, modgr2 = 2.06E-10); алгоритм Марквардта гарантированно сошелся к минимуму при значении параметра С=5 за 39 шагов. Слишком высокие значения коэффициента С (по сравнению с абсолютными значениями элементов гессиана  ) могут привести к получению отрицательно определенной матрицы , что отрицательно скажется на сходимости данного метода.

 

Литература

1. Медынский М. М., Дьячук А. К. Численные методы оптимизации с использованием системы Maple 11. – М.: Изд-во МАИ-ПРИНТ, 2009. 287 c., ил.
2. Дьяконов В. П. MATLAB. Полный самоучитель.– М.: ДМК Пресс, 2012. – 768 с., ил.
3. Конспект лекций по курсу «Численные методы оптимизации сложных систем». Автор: доцент, к.т.н. Дьячук А.К.