% !TEX encoding = UTF-8 Unicode
\documentclass[a4paper, 12pt]{article}
\usepackage[left=20mm, top=15mm, right=10mm, bottom=15mm]{geometry}
\usepackage[parfill]{parskip}
\usepackage{graphicx}
\usepackage[14pt]{extsizes}
\usepackage{setspace,amsmath}
\usepackage{mathtools}
\usepackage{ dsfont }
\usepackage{amsmath,amssymb}
\usepackage[unicode]{hyperref}
\usepackage{xcolor}
\usepackage{color}
\usepackage{minted}
\usepackage{caption}
\usepackage{array}
\newcolumntype{P}[1]{>{\centering\arraybackslash}p{#1}}
\usepackage{cmap}
\usepackage[T2A]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[english, russian]{babel}
\usepackage{amssymb}
\begin{document}
\begin{titlepage}
\thispagestyle{empty}
\begin{center}
Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования Московский государственный технический университет имени Н.Э. Баумана
\end{center}
\vfill
\centerline{\large{Лабораторная работа №5. Вариант 8.}}
\centerline{\large{«Методы поиска условного экстремума»}}
\centerline{\large{по курсу}}
\centerline{\large{«Методы оптимизации»}}
\vfill
Студент группы ИУ9-81 \hfill Ковинько А.В.
Преподаватель \hfill Каганов Ю.T.
\vfill
\centerline{Москва, 2019}
\clearpage
\end{titlepage}
\newpage
\setcounter{page}{2}
\tableofcontents
\newpage
\section{Цель работы}
\begin{enumerate}
\item Изучение алгоритмов условной оптимизации.
\item Разработка программ реализации алгоритмов условной оптимизации.
\item Нахождение оптимальных условий решений для задач с учетом ограничений.
\end{enumerate}
\newpage
\section{Постановка задачи}
\textbf {Дано}: Функция Розенброка на множестве $R^2$:
\begin{equation*}
f(x) = \sum^{n-1}_{i=1}{[a(x_i^2 - x_{i+1})^2 + b(x_i - 1)^2] + f_0},
\end{equation*}
где
\begin{equation*}
a = 158, \quad b = 2, \quad f_0 = 40, \quad n = 2,
\end{equation*}
тогда функция $f(x)$ будет выглядеть следующим образом:
\begin{equation*}
f(x) = 158 * (x_0^2 - x_1)^2 + 2 * (x_0 - 1)^2 + 40
\end{equation*}
Функции ограничений:
\begin{equation*}
\begin{cases}
g_1(x_1, x_2) = x_1^2+x_2^2 - 10 \leq 0 \\
g_2(x_1, x_2) = -x_1 \leq 0 \\
g_3(x_1, x_2) = -x_2 \leq 0
\end{cases}
\end{equation*}
\begin{enumerate}
\item Найти условный экстремум методами:
\begin{enumerate}
\item Штрафных функций;
\item Барьерных функций;
\item Комбинированным методом штрафных и барьерных функций;
\item Модифицированных функций Лагранжа;
\item Проекции градиента.
\end{enumerate}
\item Найти все стационарные точни и значения функций, соотвестсвующие этим точкам.
\item Оценить скорость сходимости указанных алгоритмов.
\item Реализовать алгоритмы с помощью языка программирования высокого уровня.
\end{enumerate}
\newpage
\section{Исследование}
Найдем глобальные экстремумы функции
\begin{equation*}
f(x) = 158 * (x_0^2 - x_1)^2 + 2 * (x_0 - 1)^2 + 40
\end{equation*}
с помощью сервиса WolframAlpha:
\begin{equation*}
min(f(x)) = 40,\quad (x_0, x_1) = (1, 1)
\end{equation*}
\begin{center}
\begin{minipage}{1\linewidth}
\includegraphics[width=\linewidth]{plt.png}
\captionof{figure}{График функции $f(x)$}
\end{minipage}
\end{center}
\subsection{Метод штрафных функций.}
Идея метода заключается в сведении задачи на условный минимум к решению последовательности задач поиска безусловного минимума вспомогательной функции:
\begin{equation}
F(x, r^k) = f(x) + P(x, r^k) \rightarrow \min_{x \in R^n},
\end{equation}
где $P(x, r^k)$ - штрафная функция, $r^k$ - параметр штрафа, задаваемый на каждой k- й итерации. Это связано с возможностью применения эффективных и надежных методов поиска безусловного экстремума,
\subsection{Метод барьерных функций.}
Идея метода заключается в сведении задачи на условный минимум к решению последовательности задач поиска безусловного минимума вспомогательной функции:
\begin{equation}
F(x, r^k) = f(x) + P(x, r^k) \rightarrow \min_{x \in R^n},
\end{equation}
где $P(x, r^k)$ - штрафная функция, $r^k$ - параметр штрафа. Используется обратная штрафная функция $P(x, r^k) = -r^k \sum_{j=1}^m \frac{1}{g_i(x)}$.
\subsection{Метод модифицированных функций Лагранжа.}
Стратегия аналогична используемой в методе внешних штрафов, только штрафная функция добавляется не к целевой функции, а к классической функции Лагранжа. В результате задача на условный минимум сводится к решению последовательности задач поиска безусловного минимума модифицированной функции Лагранжа:
\begin{multline}
L(x, \lambda ^k, \mu ^k, r^k) = f(x) + \sum_{j=1}^l \lambda_j^k g_j(x) + \frac{r^k}{2} \sum_{j=1}^{l} g_j^2(x) + \\ \frac{1}{2r^k} \sum_{j=l+1}^{m} (max^2(0, \mu_j^k + r^k g_j(x)) - \mu^2)
\end{multline}
где $\lambda ^k$ - векторы множителей Лагранжа; $r^l$ - параметр штрафа; k - номер итерации.
\subsection{Метод проекции градиента.}
Стратегия поиска решения задачи учитывает тот факт, что решение $x^*$ может лежать как внутри, так и на границе множества допустимых решений.
Для определения приближенного решения $x^*$ строится последовательность точек
\begin{equation}
\{x^*\}: x^{k+1} = x^k + \delta x^k, \quad k=1,..,
\end{equation}
где приращение $\delta x^k$ определяется в каждой точке $x^k$ в зависимости от того, где ведется поиск – внутри или на границе множества допустимых решений.
\newpage
\section{Практическая реализация}
Все методы были реализованы на языке программирования \textbf{Python}.
\textbf{Листинг 1.} Метод штрафных функций.
\begin{minted}[frame=single, framesep=10pt, fontsize = \footnotesize, linenos=false, breaklines]{python}
def p(x, r):
return r/2 * (max(0, (x[0] ** 2 + x[1] ** 2 - 10)) ** 2 + max(0, (-x[0]))** 2 + max(0, -x[1]) ** 2)
def penalty_out_methods(f, x0, r0, c=10, e=10**-8, max_iter=10):
r = r0
k = 0
while k < max_iter:
k += 1
help_f = lambda x: f(x) + p(x, r)
search_res = optimize.minimize(help_f, x0, method='CG')
x0 = search_res.x
res_p = p(x0, r)
if res_p <= e:
break
r *= c
return [k, x0, f(x0)]
\end{minted}
\textbf{Листинг 2.} Метод барьерных функий.
\begin{minted}[frame=single, framesep=10pt, fontsize = \footnotesize, linenos=false, breaklines]{python}
def p_in(x, r):
return -r*(1/(x[0] ** 2 + x[1] ** 2 - 10) - 1/x[0] - 1/x[1])
def penalty_in_method(f, x0, r0, c=10, e=10**-8, max_iter=100):
r = r0
x = x0
k = 0
while k < max_iter:
k += 1
help_f = lambda x: f(x) + p_in(x, r)
search_res = optimize.minimize(help_f, x, method='CG')
x = search_res.x
res_p = p_in(x, r)
if res_p < e:
break
r /= c
return [k, x, f(x)]
\end{minted}
\textbf{Листинг 3.} Метод модифицированных функций Лагранжа.
\begin{minted}[frame=single, framesep=10pt, fontsize = \footnotesize, linenos=false, breaklines]{python}
def lagrange(f, x0, r0=1, c=3, e=10 ** -8, max_iter=1000):
g = g_array()
r = r0
u = np.random.random(len(g))
x = x0
k = 0
while k < max_iter:
k += 1
l = lambda x: f(x) + 1 / (2 * r) * sum([max([0, u[i] + r * gf(x)]) ** 2 - u[i] ** 2 for i, gf in enumerate(g)])
search_res = optimize.minimize(l, x, method='CG')
next_x = search_res.x
x = next_x
p = 0.5 * sum([max([0, u[i] + r * gf(x)]) ** 2 - u[i] ** 2 for i, gf
in enumerate(g)]) / r
if p <= e:
return [k, x, f(x)]
r *= c
u = [max(0, u[i] + r * gf(x)) for i, gf in enumerate(g)]
return [k, x, f(x)]
\end{minted}
\textbf{Листинг 4.} Комбинированный метод штрафных и барьерных функций.
\begin{minted}[frame=single, framesep=10pt, fontsize = \footnotesize, linenos=false, breaklines]{python}
def penalty_barrier(f, x, e=10 ** -8, max_iter=1000):
r_in = 0.5
r_out = 1
c = 10
k = 0
while k < max_iter:
k += 1
combined = lambda x: f(x) + p_out(x, r_out) + p_in(x, r_in)
x = optimize.minimize(combined, x, method='CG').x
if (p_out(x, r_out) + p_in(x, r_in)) <= e:
break
r_in /= c
r_out *= c
return [k, x, f(x)]
\end{minted}
\textbf{Листинг 5.} Метод проекции градиента.
\begin{minted}[frame=single, framesep=10pt, fontsize = \footnotesize, linenos=false, breaklines]{python}
def gradient(f, x0, e=10 ** -8, max_iter=100):
x = x0
k = 0
case = 0
deltax = 0.0001
g = g_array()
dg = g_der()
gradf = grad(x)
id = np.eye(3)
np.delete(id, 2)
while k < max_iter:
for gf in g:
if -e <= gf(x) and gf(x) <= 0:
case = 1
break
if case:
gf = gradf(x)
A = dg(x)
if gf == 0:
la = np.matmul(np.matmul(inv(np.matmul(dg(x), A.transpose())), A),
np.array([gradf(x).tolist()])
.transpose())
.transpose()
.reshape(-1)
if la <= 0:
return [k, x, f(x)]
max_la = 0
max_i = 0
for i, elem in enumerate(la):
if elem < 0:
max_la = min([max_la, elem])
max_i = i
np.delete(A, max_i)
else:
deltax = - inv(id - np.matmul(A.transpose(), np.matmul(A, A.transpose())))
if deltax < e:
la = np.matmul(np.matmul(inv(np.matmul(dg(x), A.transpose())), A), np.array([gradf(x).tolist()])
.transpose())
.transpose()
.reshape(-1)
if la <= 0:
return [k, x, f(x)]
max_la = 0
max_i = 0
for i, elem in enumerate(la):
if elem < 0:
max_la = min([max_la, elem])
max_i = i
np.delete(A, max_i)
helper = lambda alpha: x + alpha * deltax
alpha = optimize.minimize(helper, x, 'CG')
x = x + alpha * deltax
return [k, x, f(x)]
\end{minted}
\newpage
\section{Результаты.}
При последовательном запуске всех алгоритмов со следующими параметрами $\epsilon = 10^{-8}$ были получены следующие результаты:
\textbf{Листинг 6.} Результаты выполнения программ.
\begin{minted}[frame=single, framesep=10pt, fontsize = \footnotesize, linenos=false, breaklines]{text}
Метод штрафных функций: [1, array([1.00000329, 1.00000666]), 40.000000000022816]
Метод внутренних штрафов[10, array([0.99999609, 0.99999217]), 40.00000000003059]
Метод Лагранжа: [1, array([0.99999985, 0.9999997 ]), 40.00000000000004]
Проекция градиента: [8, array([0.99999276, 0.99998549]), 40.000000000105025]
Комбинированный метод штрафных и барьерных функций: [10, array([1.00000117, 1.00000236]), 40.00000000000279]
\end{minted}
Все результаты с небольшой погрешностью совпадают с результатами полученными с помощью сервиса WolframAlpha.com в пункте 3.
\end{document}