% !TeX spellcheck = ru_RU

\documentclass{cmi}
\usepackage{algorithm}
\floatname{algorithm}{Алгоритм}
\usepackage{algpseudocode}
\usepackage{bm}
\usepackage{multicol}
\usepackage{upgreek}
\usepackage[justification=centering]{caption}

% Данные поля заполняются редакцией журнала
%\makeatletter
%\newcommand{\@volume}{}
%\newcommand{\@issue}{}
%\newcommand{\@num}{}
%\newcommand{\@pages}{}

% Пользовательские команды
\newcommand*\ceq{\mathop{\vcenter{\hbox{:}}{=}}}
\newcommand\doubleplus{\ensuremath{\mathbin{+\mkern-10mu+}}}
\newcommand{\algrule}[1][.2pt]{\par\vskip.01\baselineskip\hrule height #1\par\vskip.01\baselineskip}

% Подключение bib файлов
\addbblresources{Bibliography.bib}{Bibliography_sec.bib}

\begin{document}

% ЗАПОЛНЯЕТСЯ АВТОРАМИ: Релевантные коды УДК статьи
% (см. http://teacode.com/online/udc/)
\classify{УДК 519.852, 004.021, 004.032.24}
% 519.852		Линейное программирование
% 004.021		Алгоритмы
% 004.032.24	Режим обработки данных: параллельный, одновременный 

\renewcommand*{\thefootnote}{\fnsymbol{footnote}}

\title{ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ МЕТОДА ПОВЕРХНОСТНОГО ДВИЖЕНИЯ ДЛЯ РЕШЕНИЯ ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ}

% AlFaMove: Scalable Implementation of Surface Movement Method for Cluster Computing Systems

\authors{Л.Б.~Соколинский, Н.А.~Ольховский, И.М.~Соколинская}

\organizations{Южно-Уральский государственный университет\\
	(454080 Челябинск, пр.~им.~В.И.~Ленина, д.~76)}

\newcommand{\emailsinsert}{
\href{mailto:leonid.sokolinsky@susu.ru}{leonid.sokolinsky@susu.ru},
\href{mailto:olkhovskiina@susu.ru}{olkhovskiina@susu.ru}, \href{mailto:irina.sokolinskaya@susu.ru}{irina.sokolinskaya@susu.ru}
}

% ЗАПОЛНЯЕТСЯ АВТОРАМИ: Дата поступления статьи в редакцию 
\newcommand{\rdate}{10.06.2024}
\emails{\emailsinsert}
\redactiondate{\rdate}

\renewcommand*{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}

% Оформление аннотации
\begin{abstract}%
Работа посвящена численной реализации нового метода линейного программирования, получившего название <<метод поверхностного движения>>. В основе реализации лежит оригинальный алгоритм AlFaMove, который строит на поверхности допустимого многогранника оптимальный целевой путь от произвольной граничной точки до точки, являющейся решением задачи линейного программирования. Оптимальность пути заключается в том, что направление движения по грани многогранника соответствует максимальному увеличению значения целевой функции. Для вычисления оптимального направления движения используется метод, базирующийся на операции построения псевдопроекции на линейное многообразие. Операция псевдопроекции обобщает понятие ортогональной проекции и реализуется с помощью итерационного алгоритма проекционного типа. Доказано, что в случае линейного многообразия, образуемого путем пересечения гиперплоскостей, псевдопроекция совпадает с ортогональной проекцией. Также доказано, что в случае линейного многообразия метод на основе псевдопроектирования вычисляет вектор движения в направлении максимального увеличения целевой функции. Выполнена параллельная реализация алгоритма AlFaMove. Приведены результаты вычислительных экспериментов на кластерной вычислительной системе, демонстрирующие высокую масштабируемость предложенной численной реализации.

\keywords{линейное программирование, метод поверхностного движения, численная реализация, алгоритм AlFaMove, параллельная реализация, кластерная вычислительная система, исследование масштабируемости.}
\end{abstract}

%Образец цитирования статьи
\forcitation{Соколинский~Л.Б., Ольховский~Н.А., Соколинская~И.М. Численная реализация метода поверхностного движения для решения задач линейного программирования}

%Верхние колонтитулы
\markboth{Л.Б.~Соколинский, Н.А.~Ольховский, И.М.~Соколинская}{Численная реализация метода поверхностного движения}

\section*{Введение}

Эпоха больших данных и индустрия 4.0 породили задачи линейного программирования (ЛП) сверхбольших размерностей, включающих в себя миллионы переменных и миллионы ограничений~\bblcite{Fathi2019:Eng, Kopanos2019:Eng, Schlenkrich2023:Eng, Sokolinskaya2017a}. Во многих случаях объектом линейного программирования являются задачи, связанные с оптимизацией нестационарных процессов~\bblcite{Branke2002:Eng}. В нестационарных задачах ЛП целевая функция и/или ограничения изменяются в течение вычислительного процесса. Также среди этого класса задач встречаются приложения, в которых необходимо выполнять оптимизацию в режиме реального времени. Для решения таких задач необходимы масштабируемые методы и параллельные алгоритмы линейного программирования. 

Один из стандартных подходов к решению нестационарных задач оптимизации состоит в том, чтобы рассматривать каждое изменение как появление новой задачи оптимизации, которую необходимо решать с нуля~\bblcite{Branke2002:Eng}. Однако такой подход часто непрактичен, поскольку решение проблемы с нуля без повторного использования информации из прошлого может занять слишком много времени. Таким образом, желательно иметь алгоритм оптимизации, способный непрерывно адаптировать решение к изменяющейся среде, повторно используя информацию, полученную в прошлом. Этот подход применим для процессов реального времени, если алгоритм достаточно быстро отслеживает траекторию движения оптимальной точки. В случае больших задач ЛП последнее требует разработки масштабируемых методов и параллельных алгоритмов ЛП.

До настоящего времени наиболее популярными методами решения задач ЛП являются симплекс-метод~\bblcite{Dantzig1998:Eng} и методы внутренних точек~\bblcite{Zorkaltsev2018}. Эти методы способны решать задачи с десятками тысяч переменных и ограничений. Однако масштабируемость параллельных алгоритмов, основанных на симплекс-методе, в общем случае ограничивается~\mbox{16--32} процессорными узлами~\bblcite{Mamalis2015:Eng}. Что касается алгоритмов внутренних точек, они не поддаются эффективному распараллеливанию в общем случае. Это ограничивает применение указанных методов для решения сверхбольших нестационарных задач ЛП в режиме реального времени. В соответствии с этим задача разработки масштабируемых методов и эффективных параллельных алгоритмов ЛП для кластерных вычислительных систем остается актуальной.

В недавней работе~\bblcite{Olkhovsky2023a} было дано теоретическое описание нового метода ЛП --- метода поверхностного движения, строящего на поверхности допустимого многогранника\footnote{Допустимый многогранник --- область допустимых решений задачи линейного программирования.} оптимальный целевой путь к решению задачи ЛП. Под оптимальным целевым путем понимается путь по поверхности допустимого многогранника в направлении наибольшего увеличения значений целевой функции. Однако предложенный в этой статье алгоритм~1 на шаге 15 требует нахождения на границе гипердиска точки с максимальным значением целевой функции. При этом не приводится численный алгоритм, позволяющий выполнить этот шаг. В этой статье мы приводим и исследуем алгоритм AlFaMove, устраняющий допущенный пробел.

Статья организована следующим образом. Раздел~\ref{sec:theoretical_basis} содержит теоретический базис, необходимый для описания метода поверхностного движения и его численной реализации. Раздел~\ref{sec:pseudoprojection} посвящен описанию операции псевдопроекции, позволяющий найти вектор движения по оптимальному целевому пути для линейного многообразия, получаемого в результате пересечением гиперплоскостей. В разделе~\ref{sec:face-movement-algorithm} дается формализованное описание алгоритма AlFaMove, представляющего собой численную реализацию метода поверхностного движения. Раздел~\ref{sec:parallel-version} посвящен описанию параллельной версии алгоритма AlFaMove. В разделе~\ref{sec:implementation_end_experiments} представлены информация о программной реализации алгоритма AlFaMove и результаты экспериментов на кластерной вычислительной системе по исследованию его масштабируемости. В заключении суммируются полученные результаты и намечаются направления дальнейших исследований.

\section{Теоретический базис}\label{sec:theoretical_basis}

Данный раздел содержит необходимый теоретический базис, используемый для описания алгоритма движения по граням AlFaMove. Рассмотрим задачу ЛП в следующем виде:

\begin{equation}\label{eq:LP_problem}
	\bm{\bar{x}} = \arg \mathop {{\text{max}}}\limits_{\bm{x}\in\mathbb{R}^n} \left\{ {\left\langle \bm{c},\bm{x} \right\rangle |A\bm{x} \leqslant \bm{b}} \right\},	
\end{equation}
где $\bm{c}\in \mathbb{R}^n$, $\bm{b} \in \mathbb{R}^m$, $A \in \mathbb{R}^{m\times n}$, $m>1$, $\bm{c}\ne \bm{0}$. Здесь $\left\langle \cdot,\cdot \right\rangle$ обозначает скалярное произведение двух векторов. Мы предполагаем, что ограничение $\bm{x} \geqslant \bm{0}$ также включено в матричное неравенство $A\bm{x} \leqslant \bm{b}$ в форме
\begin{equation*}
	-\bm{x} \leqslant \bm{0}.
\end{equation*}
Линейная целевая функция задачи~\eqref{eq:LP_problem} имеет вид
\begin{equation*}
	f(\bm{x})=\left\langle \bm{c},\bm{x}\right\rangle.
\end{equation*}
Вектор $\bm{c}$ в данном случае является градиентом целевой функции $f(\bm{x})$.

Пусть $\bm{a}_i\in\mathbb{R}^n$ обозначает вектор, представляющий $i$-тую строку матрицы~$A$. Мы предполагаем, что $\bm{a}_i\ne \bm{0}$ для всех $i\in\lbrace 1, \dots, m \rbrace$. Обозначим через $\hat H_i$ замкнутое полупространство, определяемое неравенством $\left\langle \bm{a}_i,\bm{x}\right\rangle \leqslant b_i$, а через $H_i$ --- ограничивающую его гиперплоскость:
\begin{equation}\label{eq:half-space}
	\hat H_i=\lbrace \bm{x}\in\mathbb{R}^n | \left\langle  \bm{a}_i,\bm{x} \right\rangle \leqslant b_i \rbrace;
\end{equation}
\begin{equation}\label{eq:H_i}
	H_i=\lbrace \bm{x}\in\mathbb{R}^n | \left\langle  \bm{a}_i,\bm{x} \right\rangle = b_i \rbrace.
\end{equation}

Определим допустимый многогранник
\begin{equation}\label{eq:polytope_M}
	M = \bigcap\limits_{i = 1}^m {{\hat H}_i},
\end{equation}
представляющий множество допустимых точек задачи ЛП~\eqref{eq:LP_problem}. Заметим, что $M$ в этом случае будет замкнутым выпуклым множеством. Мы будем предполагать, что множество $M$ является ограниченным и $M\ne\emptyset$, то есть задача ЛП~\eqref{eq:LP_problem} имеет решение.

Дадим определение рецессивного полупространства~\bblcite{Sokolinsky2023}.

\begin{definition}\label{def:recessive_halfspace} 
	Полупространство $\hat H_i$ называется рецессивным, если
	\begin{equation}\label{eq:recessive_half-space}
		\forall \bm{x}\in {H_i},\forall \lambda>0:\bm{x}+\lambda \bm{c}\notin \hat H_i.
	\end{equation}
\end{definition}
\noindent Геометрический смысл этого определения состоит в том, что луч, исходящий в направлении вектора~$\bm{c}$ из любой точки гиперплоскости, ограничивающей рецессивное полупространство, не имеет общих точек с этим полупространством, за исключением начальной. Известно~\bblcite{Sokolinsky2023}, что следующее условие является необходимым и достаточным для того, чтобы полупространство $\hat H_i$ было рецессивным:
\begin{equation*}%\label{eq:recessive_half-space_condition}
	\langle \bm{a}_i,\bm{c}\rangle>0.
\end{equation*}
Определим
\begin{equation}\label{eq:I}
	\mathcal{I}=\left\lbrace i\in\lbrace 1, \dots, m \rbrace \left| \left\langle \bm{a}_i,\bm{c}\right\rangle >0 \right. \right\rbrace,
\end{equation}
то есть $\mathcal{I}$ представляет множество индексов, для которых полупространство $\hat H_i$ является рецессивным. Поскольку допустимый многогранник $M$ представляет собой ограниченное множество, имеем
\begin{equation*}%\label{Ic_ne_eptyset}
	\mathcal{I}\ne\emptyset.
\end{equation*}
Положим
\begin{equation}\label{eq:polytope_hat_M}
	\hat{M} = \bigcap\limits_{i \in \mathcal{I}} \hat H_i.
\end{equation}
Очевидно, что $\hat{M}$ является выпуклым, замкнутым, неограниченным многогранником. Будем называть его рецессивным. Из~\eqref{eq:polytope_M} и~\eqref{eq:I} следует
\begin{equation*}%\label{eq:M_subset_hat_M}
	M \subset \hat{M}.
\end{equation*}
Обозначим через $\Gamma(M)$ множество граничных точек допустимого многогранника $M$, а через $\Gamma(\hat M)$ --- множество граничных точек рецессивного многогранника $\hat{M}$\footnote{Под граничной точкой множества $\hat{M}\subset\mathbb{R}^n$ понимается точка в $\mathbb{R}^n$, для которой любая открытая ее окрестность в $\mathbb{R}^n$ имеет непустое пересечение как с множеством  $\hat{M}$, так и с его дополнением.}. Согласно утверждению~3 в~\bblcite{Sokolinsky2023} имеем
\begin{equation*}%\label{eq:}
	\bar{\bm{x}} \in \Gamma(\hat M),
\end{equation*}
то есть решение задачи ЛП~\eqref{eq:LP_problem} лежит на границе рецессивного многогранника $\hat{M}$.

Следуя~\bblcite{Maltsev1970:Rus}, дадим определение ортогональной проекции на гиперплоскость.
\begin{definition}\label{def:orthogonal_projection_on_hyperplane} 
	Пусть в пространстве $\mathbb{R}^n$ имеется гиперплоскость
	\begin{equation*}%\label{eq:}
		H = \left\lbrace \bm{x}\in\mathbb{R}^n \middle|\left\langle  \bm{a},\bm{x} \right\rangle = b \right\rbrace.
	\end{equation*}
	Ортогональная проекция $\bm{\pi}_H(\bm{v})$ точки $\bm{v} \in \mathbb{R}^n$ на гиперплоскость $H$ определяется формулой
	\begin{equation}\label{eq:pi_H}
		\bm{\pi}_H(\bm{v}) = {\bm{v}} - \frac{\left\langle \bm{a},\bm{v} \right\rangle - b} {{\| \bm{a} \|}^2}\bm{a}.
	\end{equation}
\end{definition}

Следующее утверждение дает способ вычисления оптимального пути на гиперплоскости.

\begin{proposition}\label{prp:optimal_path_on_hyperplane}
	Пусть в пространстве $\mathbb{R}^n$ задана гиперплоскость $H$ с нормалью $\bm{a}\in\mathbb{R}^n$, проходящая через точку $\bm{u}\in\mathbb{R}^n$:
	\begin{equation}\label{eq:H}
		H = \left\{ {{\bm{x}} \in {\mathbb{R}^n}\left| {\left\langle {{\bm{a}},{\bm{x}}} \right\rangle  = \left\langle {{\bm{a}},{\bm{u}}} \right\rangle } \right.} \right\}.
	\end{equation}
	Пусть задана линейная функция $f(\bm{x}):\mathbb{R}^n \to \mathbb{R}$ с градиентом $\bm{c}\in\mathbb{R}^n$:
	\begin{equation}\label{eq:f(x)}
		f(\bm{x})=\left\langle \bm{c},\bm{x}\right\rangle.
	\end{equation}
	Пусть векторы $\bm{a}$ и $\bm{c}$ линейно независимы (не коллинеарны, и среди них нет нулевого вектора). Положим
	\begin{equation}\label{eq:v}
		\bm{v} = \bm{u} + \bm{c}.
	\end{equation}
	Построим ортогональную проекцию $\bm{\pi}_H(\bm{v})$ точки $\bm{v}$ на гиперплоскость $H$:
	\begin{equation}\label{eq:w=pi_H(v}
		\bm{w} = \bm{\pi}_H(\bm{v}).
	\end{equation}
	Тогда вектор $\bm{d} = \bm{w}-\bm{u}$ однозначно задает направление максимального увеличения линейной функции $f(\bm{x})$, определяемой формулой~\eqref{eq:f(x)}.
\end{proposition}
\begin{proof} 
	Предположим противное, то есть существует точка $\bm{\tilde w} \in H$ такая, что  
	\begin{equation}\label{eq:<c,dot_w>>=<c,w>}
		\left\langle \bm{c}, \bm{\tilde w} \right\rangle \geqslant \left\langle \bm{c}, \bm{w} \right\rangle,
	\end{equation}
	$\| \bm{\tilde w} - \bm{u}\| = \| \bm{w} - \bm{u}\|$ и $\bm{\tilde w} \ne \bm{w}$ (см.~рис.~\ref{fig:Optimal_Path_on_Hyperplane}). Здесь и далее $\|  {}\cdot{}  \|$ обозначает евклидову норму.	
	%
	\begin{figure}[t]
				\includegraphics[height=4.5 cm]{Fig_Optimal_Path_on_Hyperplane}
				\caption{Иллюстрация к доказательству предложения~\ref{prp:optimal_path_on_hyperplane}. \\
					Пунктиром обозначена гиперокружность с радиусом $\| \bm{w} - \bm{u} \| $.}
				\label{fig:Optimal_Path_on_Hyperplane}
	\end{figure}
	%
	Вычислим $\left\langle \bm{c}, \bm{w} \right\rangle$. В соответствии с определением~\ref{def:orthogonal_projection_on_hyperplane} ортогональная проекция $\bm{\pi}_H(\bm{v})$ точки $\bm{v}$ на гиперплоскость $H$, задаваемую формулой~\eqref{eq:H}, вычисляется следующим образом:
	\begin{equation*}%\label{eq:}
		{\bm{w}} = {\bm{v}} - \frac{{\left\langle {{\bm{a}},{\bm{v}} - {\bm{u}}} \right\rangle }}{{{{\| {\bm{a}} \|}^2}}}{\bm{a}}.
	\end{equation*}
	Подставив вместо $\bm{v}$ правую часть формулы~\eqref{eq:v}, получим
	\begin{equation}\label{eq:w}
		{\bm{w}} = {\bm{u}} + {\bm{c}} - \frac{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }}{{{{\| {\bm{a}} \|}^2}}}{\bm{a}}.
	\end{equation}
	Используя~\eqref{eq:w}, находим
	\begin{equation}\label{eq:<c,w>}
		\left\langle {{\bm{c}},{\bm{w}}} \right\rangle  = \left\langle {{\bm{c}},{\bm{u}}} \right\rangle  + {\| {\bm{c}} \|^2} - \frac{{{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }^2}}}{{{{\| {\bm{a}} \|}^2}}}.
	\end{equation}
	Поскольку $\bm{a}$ и $\bm{c}$ линейно независимы, в соответствии с неравенством Коши--Буняковского имеем
	\begin{equation*}%\label{eq:}
		{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle ^2} < {\| {\bm{a}} \|^2} \cdot {\| {\bm{c}} \|^2}.
	\end{equation*}
	Отсюда следует
	\begin{equation}\label{eq:positive_term}
		{\| {\bm{c}} \|^2} - \frac{{{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }^2}}}{{{{\| {\bm{a}} \|}^2}}} > 0.
	\end{equation}
	
	Теперь вычислим $\left\langle \bm{c}, \bm{\tilde w} \right\rangle$. Определим $\bm{\tilde u} = \bm{\bm{\pi}_L}(\bm{\tilde w})$ --- ортогональная проекция точки $\bm{\tilde w}$ на прямую~$L$, проходящую через точки $\bm{u}$ и $\bm{w}$. По построению существует число $\delta \in \mathbb{R}$, удовлетворяющее условию 	
	\begin{equation}\label{eq:delta}
		-1 \leqslant \delta < 1,  
	\end{equation}
	такое, что 
	\begin{equation*}%\label{eq:}
		{\bm{\tilde u}} = \bm{u} + \delta ({\bm{w}} - {\bm{u}}).
	\end{equation*}
	Положим
	\begin{equation}\label{eq:dot_v}
		\bm{\tilde v} = \bm{u} + \delta (\bm{v} - \bm{u}).
	\end{equation}
	Тогда точка $\bm{\tilde u}$ является ортогональной проекцией точки $\bm{\tilde v}$ на гиперплоскость $H$, определяемую формулой~\eqref{eq:H}, и может быть вычислена следующим образом:
	\begin{equation*}%\label{eq:}
		{\bm{\tilde u}} = {\bm{\tilde v}} - \frac{{\left\langle {{\bm{a}},{\bm{\tilde v}} - {\bm{u}}} \right\rangle }}{{{{\| {\bm{a}} \|}^2}}}{\bm{a}}.
	\end{equation*}
	Подставив вместо $\bm{\tilde v}$ правую часть формулы~\eqref{eq:dot_v}, отсюда получим
	\begin{equation*}%\label{eq:dot_w}
		{\bm{\tilde u}} = {\bm{u}} + \delta \left( {{\bm{v}} - {\bm{u}} - \frac{{\left\langle {{\bm{a}},{\bm{v}} - {\bm{u}}} \right\rangle }}{{{{\| {\bm{a}} \|}^2}}}{\bm{a}}} \right).
	\end{equation*}
	Используя~\eqref{eq:v}, произведем замену $\bm{v}-\bm{u} = \bm{c}$:
	\begin{equation}\label{eq:dot_u}
		{\bm{\tilde u}} = {\bm{u}} + \delta \left( {{\bm{c}} - \frac{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }}{{{{\| {\bm{a}} \|}^2}}}{\bm{a}}} \right).
	\end{equation}
	Очевидно, 
	\begin{equation}\label{eq:dot_w}
		\bm{\tilde w} = (\bm{\tilde w} - \bm{\tilde u}) + \bm{\tilde u}.
	\end{equation}
	Заменим второе слагаемое в~\eqref{eq:dot_w} на правую часть формулы~\eqref{eq:dot_u}:
	\begin{equation*}%\label{eq:}
		\bm{\tilde w} = (\bm{\tilde w} - \bm{\tilde u}) + {\bm{u}} + \delta \left( {{\bm{c}} - \frac{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }}{{{{\| {\bm{a}} \|}^2}}}{\bm{a}}} \right).
	\end{equation*}
	Отсюда получаем
	\begin{equation*}%\label{eq:}
		\left\langle {{\bm{c}},{\bm{\tilde w}}} \right\rangle  = \left\langle {{\bm{c}},{\bm{\tilde w}} - {\bm{\tilde u}}} \right\rangle  + \left\langle {{\bm{c}},{\bm{u}}} \right\rangle  + \delta \left( {{{\| {\bm{c}} \|}^2} - \frac{{{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }^2}}}{{{{\| {\bm{a}} \|}^2}}}} \right).
	\end{equation*}
	По построению вектор $\bm{\tilde w} - \bm{\tilde u}$ ортогонален вектору $\bm{v}-\bm{u} = \bm{c}$. Поэтому $\left\langle {{\bm{c}},{\bm{\tilde w}}- \bm{\tilde u}} \right\rangle = 0$. Таким образом, последняя формула преобразуется к виду
	\begin{equation}\label{eq:<c,dot_w>}
		\left\langle {{\bm{c}},{\bm{\tilde w}}} \right\rangle  = \left\langle {{\bm{c}},{\bm{u}}} \right\rangle  + \delta \left( {{{\| {\bm{c}} \|}^2} - \frac{{{{\left\langle {{\bm{a}},{\bm{c}}} \right\rangle }^2}}}{{{{\| {\bm{a}} \|}^2}}}} \right).
	\end{equation}
	Сопоставляя~\eqref{eq:<c,w>} и~\eqref{eq:<c,dot_w>}, с учетом~\eqref{eq:positive_term} и~\eqref{eq:delta} получаем
	\begin{equation*}%\label{eq:}
		\left\langle \bm{c}, \bm{\tilde w} \right\rangle < \left\langle \bm{c}, \bm{w} \right\rangle,
	\end{equation*}
	что противоречит~\eqref{eq:<c,dot_w>>=<c,w>}.	 
\end{proof}

Возвращаясь к задаче ЛП~\eqref{eq:LP_problem}, можно сказать следующее. Пусть $\bm{u} \in M \cap \Gamma( \hat{M})$ и существует единственная рецессивная гиперплоскость $H_{i'}$ ($i' \in \mathcal{I}$) такая, что $\bm{u} \in H_{i'}$. В этом случае вектор~$\bm{d}$, определяющий направление оптимального целевого пути из точки $\bm{u}$, в соответствии с утверждением~\ref{prp:optimal_path_on_hyperplane} вычисляется по формуле
\begin{equation}\label{eq:d}
	{\bm{d}} = {\bm{c}} - \frac{{\left\langle {{{\bm{a}}_{i'}},{\bm{u}} + {\bm{c}}} \right\rangle  - {b_{i'}}}}{{{{\| {{{\bm{a}}_{i'}}} \|}^2}}}{{\bm{a}}_{i'}}.
\end{equation}
В следующем разделе мы рассмотрим общий случай, когда через точку $\bm{u}$ проходит несколько гиперплоскостей.

\section{Операция псевдопроектирования на линейное многообразие}\label{sec:pseudoprojection}

Пусть $\mathcal{J} \subseteq \lbrace 1, \dots, m \rbrace$, $\mathcal{J} \ne \emptyset$, и $\bigcap\limits_{i \in \mathcal{J}}{H_i} \ne \emptyset$. В этом случае множество индексов $\mathcal{J}$ определяет в пространстве $\mathbb{R}^n$ линейное многообразие 
\begin{equation}\label{eq:L_cap_hyperplanes}
	L =  \bigcap\limits_{i\in \mathcal{J}}{H_i}.
\end{equation}
Обозначим $k_L$ --- размерность линейного многообразия $L$. При $k_L < n - 1$ многообразие~$L$ не является гиперплоскостью, и для определения вектора движения $\bm{d}$ по этому многообразию в направлении максимального увеличения целевой функции нельзя применить формулу~\eqref{eq:d}, так как такое линейное многообразие невозможно задать одним линейным уравнением в пространстве $\mathbb{R}^n$. Однако мы можем найти указанный вектор $\bm{d}$  с помощью операции псевдопроектирования~\bblcite{Sokolinsky2023}. Определим проекционное отображение $\bm{\varphi}(\cdot)$:
\begin{equation}\label{eq:psi}
	\bm{\varphi}(\bm{x}) = \frac{1}{|\mathcal{J}|}\sum\limits_{i\in \mathcal{J}} {\bm{\pi}_{H_i}(\bm{x})}.
\end{equation}
Известно~\bblcite{Vasin2005}, что отображение $\bm{\varphi}(\bm{x})$ является непрерывным $L$-фейеровским отображением, и последовательность
точек
\begin{equation}\label{eq:x_k}
	\mathop{\left\lbrace \bm{x}_k=\bm{\varphi}^k\left(\bm{x}_0\right)\right\rbrace}\nolimits_{k=1}^{\infty},
\end{equation}
порождаемая этим отображением, начиная с произвольной точки $\bm{x}_0 \in \mathbb{R}^n$, сходится к точке, принадлежащей $L$:
\begin{equation*}
	\bm{x}_k \to \bm{\tilde x} \in L.
\end{equation*}
Теперь мы готовы дать определение псевдопроекции на линейное многообразие, образуемое пересечением гиперплоскостей.

\begin{definition}\label{def:pseudoprojection}
	Пусть $\mathcal{J} \subseteq \lbrace 1, \dots, m \rbrace$, $\mathcal{J} \ne \emptyset$, $\bigcap\nolimits_{i \in \mathcal{J}}{H_i} \ne \emptyset$, $\bm{\varphi}(\cdot)$ --- проекционное отображение, определяемое формулой~\eqref{eq:psi}. Псевдопроекцией $\bm{\rho}_{\mathcal{J}}(\bm{x})$ точки $\bm{x}\in \mathbb{R}^n$ на линейное многообразие $L=\bigcap\nolimits_{i \in \mathcal{J}}{H_i}$ называется предельная точка последовательности~\eqref{eq:x_k}:
	\begin{equation*}
		\mathop {\lim }\limits_{k \to \infty } \left\| {\bm{\rho}_{\mathcal{J}}(\bm{x}) - {{\bm{\varphi}}^k}(\bm{x})} \right\| = 0.
	\end{equation*}	
\end{definition}

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

\begin{lemma}\label{lem:pseudoprojection=orthogonal_projection}
	Пусть в пространстве $\mathbb{R}^n$ заданы гиперплоскость $H_{i'}$ и линейное многообразие $L$, принадлежащее этой гиперплоскости:
	\begin{eqnarray*}
		&& H_{i'} = \{ \bm{x} \in {\mathbb{R}^n} |\, \langle {\bm{a}}_{i'},{\bm{x}} \rangle  = b_{i'} \};%\label{eq:}
		\\&& L=\bigcap\nolimits_{i \in \mathcal{J}}{H_i};
		\\&& i' \in \mathcal{J}.%\label{eq:} 
	\end{eqnarray*}
	Обозначим через $P$ линейное многообразие, являющееся ортогональным дополнением к $L$:
	\begin{equation}\label{eq:P}
		P = L^\bot.
	\end{equation}
	Тогда для любой точки $\bm{v}$, принадлежащей линейному многообразию $P$, ее ортогональная проекция $\bm{\pi}_{H_{i'}}(\bm{v})$ на гиперплоскость $H_{i'}$ также будет принадлежать линейному многообразию $P$: 
	\begin{equation*}%\label{eq:}
		\forall \bm{v} \in P: \bm{\pi}_{H_{i'}}(\bm{v}) \in P.
	\end{equation*}
\end{lemma}
\begin{proof} 
	Обозначим $\bm{p}$ --- ортогональная проекция точки $\bm{v}$ на гиперплоскость $H_{i'}$:
	\begin{equation} \label{eq:p=pi(v)}
	\bm{p} = \bm{\pi}_{H_{i'}}(\bm{v}).
	\end{equation}
	Без ограничения общности мы можем считать, что $\bm{p} \in L$ (см. рис~\ref{fig:Projecting_plane}).
	%
	\begin{figure}[t]
				\includegraphics[height=7 cm]{Fig_Projecting_plane}
				\caption{Иллюстрация к лемме~\ref{lem:pseudoprojection=orthogonal_projection}.}
				\label{fig:Projecting_plane}
	\end{figure}
	%
	Предположим, что точка $\bm{p}$ не принадлежит линейному многообразию $P$. Обозначим $\bm{q}$ --- точка пересечения линейного многообразия $L$ с его ортогональным дополнением $P$:
	\begin{equation*}
		\bm{q} = L \cap P.
	\end{equation*}
	Поскольку $P$ является ортогональным дополнением к линейному многообразию $L$, точка $\bm{q}$ является ортогональной проекцией точки~$\bm{p}$ на линейное многообразие $P$:
	\begin{equation}\label{eq:q}
		\bm{q} = \bm{\pi}_P(\bm{p}).
	\end{equation}
	Рассмотрим треугольник $\triangle(\bm{v},\bm{p},\bm{q})$. В силу~\eqref{eq:p=pi(v)} угол $\angle\bm{p}$ с вершиной в точке $\bm{p}$ является прямым. В в соответствии с~\eqref{eq:q} угол $\angle\bm{q}$ с вершиной в точке $\bm{q}$ также является прямым. Но это возможно только в случае, когда $\bm{p} = \bm{q}$, то есть
	\begin{equation*}%\label{eq:}
		\bm{p} \in P.
	\end{equation*}
	Лемма доказана.	
\end{proof}

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

\begin{proposition}\label{prp:pseudoprojection=orthogonal_projection}
	Пусть выполняются следующие условия:
	\begin{eqnarray}
		&& \mathcal{J} \subseteq \lbrace 1, \dots, m \rbrace,\;\label{eq:J=1,...,m;|J|=n-1}
		\\&& L = \bigcap\limits_{i \in \mathcal{J}}{H_i},\, L \ne \emptyset \label{eq:L};
	\end{eqnarray}
	где $H_i=\lbrace \bm{x}\in\mathbb{R}^n | \left\langle  \bm{a}_i,\bm{x} \right\rangle = b_i \rbrace$. Обозначим $\bm{\pi}_L(\bm{x})$ --- ортогональная проекция точки $\bm{x} \in \mathbb{R}^n$ на линейное многообразие~$L$.
	Тогда 
	\begin{equation*}
		{\bm{\rho}}_L(\bm{x}) = {\bm{\pi}}_L(\bm{x}),
	\end{equation*}	
	то есть, псевдопроекция на линейное многообразие~$L$ совпадает с ортогональной проекцией.
\end{proposition}
\begin{proof} 
	Зафиксируем произвольную точку $\bm{v}_0 \in \mathbb{R}^n$. Рассмотрим линейное многообразие $P$, содержащее точку $\bm{v}_0$, и являющееся ортогональным дополнением к линейному многообразию $L$:
	\begin{equation} \label{eq:v_0_in_P;P=D_bot}
		\bm{v}_0 \in P = L^\bot.
	\end{equation}
	Обозначим через $\bar{\bm{v}}$ точку пересечения линейного многообразия $L$ с его ортогональным дополнением $P$:
	\begin{equation*}%\label{eq:D_cap_P}
		L \cap P = \{ \bar{\bm{v}} \}
	\end{equation*}
	(см. рис.~\ref{fig:Fig_Pseudoprojection_on_edge}).
	%
	\begin{figure}[t]
				\includegraphics[height=6.5 cm]{Fig_Pseudoprojection_on_edge}
				\caption{Иллюстрация к доказательству утверждения~\ref{prp:pseudoprojection=orthogonal_projection} при  $n=3$.}
				\label{fig:Fig_Pseudoprojection_on_edge}
	\end{figure}
	%
	Построим ортогональную проекцию точки $\bm{v}_0$ на гиперплоскость $H_j$ для произвольного $j \in \mathcal{J}$:
	\begin{equation*}
		\bm{p}_0 = \bm{\pi}_{H_j}(\bm{v}_0).
	\end{equation*}
	В соответствии с леммой~\ref{lem:pseudoprojection=orthogonal_projection} имеем
	\begin{equation*}
		\bm{\pi}_{H_j}(\bm{v}_0) \in P.
	\end{equation*}
	Отсюда и из~\eqref{eq:psi} следует, что 
	\begin{equation*}
		\bm{v}_1 = \varphi(\bm{v}_0)  \in P.
	\end{equation*}
	Это означает, что  последовательность
	\begin{equation*}%\label{eq:v_k}
		\mathop{\left\lbrace \bm{v}_k=\bm{\varphi}^k\left(\bm{v}_0\right)\right\rbrace}\nolimits_{k=1}^{\infty}
	\end{equation*}
	сходится к точке $\bar{\bm{v}}$ пересечения линейного многообразия $L$ с линейным многообразием $P$, то есть, $\bm{\rho}_L(\bm{v}_0) = \bar{\bm{v}}$.
	С~другой стороны, в силу~\eqref{eq:v_0_in_P;P=D_bot} имеем $\bm{\pi}_L(\bm{v}_0) = \bar{\bm{v}}$.
	Следовательно
	\begin{equation*}
		\forall \bm{x} \in \mathbb{R}^n:\, \bm{\rho}_L(\bm{x}) = \bm{\pi}_L(\bm{x}).
	\end{equation*}
	Утверждение доказано.
\end{proof}

Процедура приближенного вычисления псевдопроекции на линейное многообразие представлена в виде алгоритма~\ref{alg:pseudoprojection}.
%
\begin{algorithm}[t]\caption{Вычисление псевдопроекции $\bm{\rho}_{\mathcal{J}}(\bm{x})$}\label{alg:pseudoprojection}				\begin{algorithmic}[1]
					\Require $H_i=\left\lbrace \bm{x}\in\mathbb{R}^n \middle|\left\langle  \bm{a}_i,\bm{x} \right\rangle = b_i \right\rbrace;\; \mathcal{J} \subseteq \lbrace 1, \dots, m \rbrace;\; \mathcal{J} \ne \emptyset;\; \bigcap\limits_{i \in \mathcal{J}}{H_i} \ne \emptyset$
					\Function {$\bm{\rho}_{\mathcal{J}}$}{$\bm{x}$}
					\State $k\ceq 0$
					\State $\bm{x}_0 \ceq \bm{x}$
					\Repeat
					\State $\Sigma\ceq 0$
					\For{$i\in\mathcal{J}$}
					\State $\Sigma \ceq \Sigma + (\langle \bm{a}_i,\bm{x}_k \rangle - b_i) \bm{a}_i / {\| \bm{a}_i \|}^2$
					\EndFor
					\State $\bm{x}_{(k+1)}\ceq \bm{x}_k - \Sigma / \left| \mathcal{J}\right|$ 
					\State $\xi_{max} \ceq 0$ \Comment{Максимальная невязка}
					\For{$i\in\mathcal{J}$}
					\State $\xi_i \ceq \| \left\langle \bm{a}_i,\bm{x}_{k+1} \right\rangle - b_i \|$
					\If{$\xi_i > \xi_{max}$}
					\State $\xi_{max} \ceq \xi_i$
					\EndIf
					\EndFor	
					\State $k\ceq k+1$
					\Until{$\xi_{max} < \epsilon_\xi$} \Comment{Малый параметр $\epsilon_\xi > 0$}
					\State \Return $\bm{x}_k$
					\EndFunction					
\end{algorithmic}\end{algorithm}
%
Дадим краткий комментарий по шагам этого алгоритма. На шаге~2 счетчик итераций $k$ устанавливается в значение ноль. Шаг~3 задает начальное приближение $\bm{x}_0$. Шаг~4 открывает итерационный цикл вычисления псевдопроекции. Шаги~5\nobreakdash--\hspace{0pt}8 для текущего приближения~$\bm{x}_k$ вычисляют сумму из правой части формулы~\eqref{eq:psi}. Шаг~9 находит следующее приближение $\bm{x}_{k+1}$. Шаги~10\nobreakdash--\hspace{0pt}16 определяют максимальную невязку $\xi$ для нового приближения относительно всех гиперплоскостей $H_i$, участвующих в вычислении. На шаге~17 счетчик итераций увеличивается на единицу. Шаг~18 проверяет условие выхода из итерационного цикла. Шаг~19 возвращает в качестве результата полученное приближение~$\bm{x}_k$.

\section{Алгоритм движения по граням AlFaMove}\label{sec:face-movement-algorithm}

В этом разделе мы опишем новый алгоритм AlFaMove (Along Faces Movement), представляющий собой численную реализацию метода поверхностного движения~\bblcite{Olkhovsky2023a}. Алгоритм AlFaMove строит на поверхности допустимого многогранника путь из произвольной граничной точки ${\bm{u}}_0\in M \cap\Gamma(\hat{M})$ до точки $\bar{\bm{x}}$, являющейся решением задачи ЛП~\eqref{eq:LP_problem}. Перемещение по граням допустимого многогранника происходит в направлении наибольшего увеличения значения целевой функции. Путь, построенный в результате такого движения, называется оптимальным целевым путем.

В основе алгоритма AlFaMove лежит процедура $\bm{D}(\bm{u})$ вычисления вектора движения $\bm{\bar{d}}$ по грани допустимого многогранника $M$ из точки  $\bm{u}$ в направлении максимального увеличения значения целевой функции. Процедура $\bm{D}(\bm{u})$ представлена в виде алгоритма~\ref{alg:direction-vector}. 
%
\begin{algorithm}[t]\caption{Вычисление вектора движения $\bm{\bar d} = \bm{D}(\bm{u})$}\label{alg:direction-vector}\begin{algorithmic}[1]
					\Require $H_i=\left\lbrace \bm{x}\in\mathbb{R}^n \middle|\left\langle  \bm{a}_i,\bm{x} \right\rangle = b_i \right\rbrace$; $\bm{u} \in \Gamma(M)$
					\Function {$\bm{D}$}{$\bm{u}$}
					\State $\mathcal{U} \ceq \emptyset$ \Comment{$\mathcal{U}$ --- множество индексов гиперплоскостей $H_i$, проходящих через точку $\bm{u}$}
					\For{$i = 1 \dots m$}
					\If{$\langle  \bm{a}_i,\bm{u} \rangle = b_i$}
					\State $\mathcal{U} \ceq \mathcal{U} \cup \lbrace i \rbrace $
					\EndIf
					\EndFor
					\State $\bm{\bar d} \ceq \bm{0}$
					\State $f \ceq -\infty$ \Comment{$f$ --- значение целевой функции $f(\bm{x}) = \langle \bm{c}, \bm{x}\rangle$}
					\State $\bm{e_c} \ceq \bm{c}/\| \bm{c} \|$
					\State $\bm{v} \ceq \bm{u} + \delta\bm{e_c}$ \Comment{Большой параметр $\delta > 0$}
					\For{$\mathcal{J}\in\mathcal{P}(\mathcal{U})\backslash\emptyset$} \Comment{$\mathcal{P}(\mathcal{U})$ --- множество всех подмножеств множества $\mathcal{U}$}
					\State $\bm{w} \ceq \bm{\rho}_{\mathcal{J}}(\bm{v})$
					\State $\bm{d} \ceq \bm{w} - \bm{u}$
					\State $\bm{e_d} \ceq \bm{d}/\| \bm{d} \|$
					\If{$(\bm{u} + \tau\bm{e_d}) \in M$} \Comment{Малый параметр $\tau > 0$}
					\If{$\left\langle \bm{c}, \bm{u} + \tau\bm{e_d}\right\rangle > f$}
					\State $f \ceq \langle \bm{c}, \bm{u} + \tau\bm{e_d}\rangle$
					\State $\bm{\bar d} \ceq \bm{d}$
					\EndIf
					\EndIf
					\EndFor
					\State \Return $\bm{\bar d}$
					\EndFunction
\end{algorithmic}\end{algorithm}
%
Схематично работа этого алгоритма изображена на рис.~\ref{fig:Direction_vector}. 
%
\begin{figure}[!ht]
	\includegraphics[height=6 cm]{Fig_Direction_vector}
	\caption{Вычисление вектора направления движения~$\bm{\bar d}$ по грани линейного многообразия~$L$.}
	\label{fig:Direction_vector}
\end{figure}
%
Дадим краткий комментарий по шагам алгоритма~\ref{alg:direction-vector}.
Шаги~2--7 строят множество $\mathcal{U}$, включающее в себя индексы всех гиперплоскостей $H_i$, проходящих через точку $\bm{u}$. На шаге~8 вектору направления движения $\bm{\bar d}$ в качестве начального значения присваивается нулевой вектор. На шаге~9 значению целевой функции $f$ присваивается бесконечно малое число\footnote{В случае 64-разрядной арифметики с плавающей точкой, бесконечно малое число равно $-1\cdot 10^{308}$}. На шаге~10 строится единичный вектор $\bm{e_c}$, сонаправленный с вектором $\bm{c}$. На шаге~11 строится точка $\bm{v}$ путем прибавления вектора $\delta\bm{e_c}$ к вектору $\bm{u}$ (см. рис.~\ref{fig:Direction_vector}). Здесь $\delta$ --- <<большой>> положительный параметр: чем больше~$\delta$, тем точнее будет вычислен вектор направления движения $\bm{\bar d}$. Однако при увеличении параметра~$\delta$ будет увеличиваться и время вычисления псевдопроекции $\bm{\rho}_{\mathcal{J}}(\bm{v})$ на шаге~13. Цикл \textbf{for} на шаге~12 перебирает все возможные комбинации индексов гиперплоскостей, проходящих через точку $\bm{u}$. Каждой такой комбинации $\mathcal{J}$ соответствует линейное многообразие $L=\bigcap\nolimits_{i \in \mathcal{J}}{H_i}$, которое также проходит через точку $\bm{u}$. Шаг~13 вычисляет точку $\bm{w}$, выполняя псевдопроекцию точки $\bm{v}$ на линейное многообразие, соответствующее комбинации $\mathcal{J}$. На шаге~14 вычисляется возможный вектор движения $\bm{d}$ по текущему линейному многообразию в направлении максимального увеличения значений целевой функции. Шаг~15 вычисляет единичный вектор $\bm{e_d}$, сонаправленный с вектором $\bm{d}$. Шаг~16 проверяет, что малое движение из точки  $\bm{u}$ в направлении $\bm{d}$ не выходит за границы допустимого многогранника. Шаг~17, в свою очередь, проверяет, что значение целевой функции в точке $(\bm{u} + \bm{e_d})$ превышает максимальное значение, полученное на предыдущих итерациях цикла \textbf{for} (шаг~12). В этом случае это значение запоминается как максимальное (шаг~18), а найденное направление присваивается вектору $\bm{\bar{d}}$ (шаг~19). После того, как проверены все возможные комбинации, вектор $\bm{\bar{d}}$ возвращается в качестве результата (шаг~23). Если ни одна комбинация не прошла проверку на шагах~16--17, то в качестве результата будет возвращен нулевой вектор. Это означает, что любое движение из точки $\bm{u}$ по поверхности допустимого многогранника не приводит к увеличению значения целевой функции.

Теперь мы готовы описать алгоритм движения по граням AlFaMove, решающий задачу ЛП~\eqref{eq:LP_problem}. За основу мы берем алгоритм~1 из работы~\bblcite{Olkhovsky2023a}. Реализация алгоритма AlFaMove на псевдокоде представлена в виде алгоритма~\ref{alg:AlFaMove}.
%
\begin{algorithm}[b]\caption{Алгоритм движения по граням AlFaMove}\label{alg:AlFaMove}\begin{algorithmic}[1]
					\Require $\hat H_i=\left\lbrace \bm{x}\in\mathbb{R}^n \middle|\left\langle  \bm{a}_i,\bm{x} \right\rangle \leqslant b_i \right\rbrace;\; M = \bigcap\limits_{i = 1}^m {{\hat H}_i};\; \hat M=\bigcap\limits_{i \in \mathcal{I}} {\hat H_i};\; i \in \mathcal{I} \Leftrightarrow 	\langle \bm{a}_i,\bm{c}\rangle>0$ 
					\State \textbf{input} ${\bm{u}}_0$
					\State \textbf{assert} ${\bm{u}}_0\in M \cap \Gamma(\hat{M})$
					\State $\bm{d}_0 \ceq \bm{D}(\bm{u}_0)$
					\State \textbf{assert} $\bm{d}_0 \ne \bm{0}$
					\State $k\ceq 0$
					\Repeat
					\State ${\bm{u}}_{k+1} \ceq {\bm{\mu}}(\bm{u}_k,\bm{d}_k)$
					\State $\bm{d}_{k+1} \ceq \bm{D}(\bm{u}_{k+1})$
					\State $k\ceq k+1$
					\Until{$\bm{d}_k = \bm{0}$}
					\State \textbf{output} $\bm{u}_k$ \Comment Решение задачи ЛП~\eqref{eq:LP_problem}
					\State \textbf{stop}
\end{algorithmic}\end{algorithm}
%
Прокомментируем шаги этого алгоритма. На шаге~1 вводится начальное приближение ${\bm{u}}^{(0)}$. Это может быть произвольная граничная точка рецессивного многогранника $\hat M$, удовлетворяющая условию
\begin{equation*}%\label{eq:u^(0)_in_hat_M_cap_Gamma(M)} 
	\bm{u}_0\in M \cap\Gamma(\hat{M}).
\end{equation*}
Это условие проверяется на шаге~2. Для получения подходящего начального приближения может применяться алгоритм, реализующий стадию Quest апекс-метода~\bblcite{Sokolinsky2023}. Шаг~3 вычисляет начальный вектор движения $\bm{d}_0$. Для этого используется функция $\bm{D}(\cdot)$, реализованная в алгоритме~\ref{alg:direction-vector}. Предполагается, что $\bm{d}_0 \ne \bm{0}$\footnote{Равенство вектора $\bm{d}_0$ нулевому вектору означает, что точка $\bm{u}_0$ является решением задачи ЛП~\eqref{eq:LP_problem}.}. Это условие контролируется на шаге~4. На шаге~5 счетчик итераций~$k$ устанавливается в значение ноль. Шаг~6 открывает цикл \textbf{repeat}, выполняющий движение по граням допустимого многогранника до тех пор, пока очередной вектор движения $\bm{d}_k$ не окажется нулевым вектором. Это означает, что последнее найденное приближение $\bm{u}_k$ является решением задачи ЛП~\eqref{eq:LP_problem}. Шаг~7 в теле цикла вычисляет следующее приближение $\bm{u}_{k+1}$ с помощью вектор-функции ${\bm{\mu}}$ (определение см. ниже). Шаг~8 вычисляет вектор движения $\bm{d}_{k+1}$ для вновь найденного приближения $\bm{u}_{k+1}$. Шаг~9 увеличивает счетчик итераций $k$ на единицу. Если последний вектор движения равен нулевому вектору, то на шаге~10 происходит выход из цикла \textbf{repeat/until}, после чего на шаге~11 выводятся координаты точки $\bm{u}_k$ в качестве решения задачи ЛП~\eqref{eq:LP_problem}. Шаг~12 завершает работу алгоритма AlFaMove.

Вектор-функция ${\bm{\mu}}(\cdot)$, используемая на шаге 7, определяется следующим образом. Обозначим
\begin{equation}\label{eq:Q}
	\mathcal{Q} =\left\lbrace i\in\lbrace 1, \dots, m \rbrace \left| \; \left\langle \bm{a}_i,\bm{u} \right\rangle < b_i \wedge \left\langle \bm{a}_i,\bm{d} \right\rangle > 0 \right. \right\rbrace.
\end{equation}
Тогда
\begin{equation}\label{eq:mu(u,d)}
	{\bm{\mu}}(\bm{u},\bm{d}) = \arg \mathop {{\text{min}}}\limits_{i\in \mathcal{Q}} \left\{ \left\| \bm{u} - \bm{x} \right\| | \; \bm{x} = \bm{\gamma}_i(\bm{u},\bm{d}) \right\}.	
\end{equation}
Здесь $\bm{\gamma}_i(\bm{u},\bm{d})$ обозначает вектор-функцию, вычисляющую косоугольную проекцию точки~$\bm{u}$ на гиперплоскость $H_i$ по направлению $\bm{d}$:
\begin{equation*}%\label{eq:}
	\bm{\gamma}_i(\bm{u},\bm{d}) = \bm{u} - \frac{\left\langle \bm{a}_i, \bm{u}\right\rangle - b_i}{\left\langle \bm{a}_i, \bm{d}\right\rangle}\bm{d}.
\end{equation*}
Действие функции ${\bm{\mu}}$ проиллюстрировано на рис.~\ref{fig:Jump}.
%
\begin{figure}[t]
	\includegraphics[height=6 cm]{Fig_Jump}
	\caption{Действие функции ${\bm{\mu}}$: \\ ${\bm{\mu}}(\bm{u},\bm{d}) = \bm{\gamma}_1(\bm{u},\bm{d})$.}
	\label{fig:Jump}
\end{figure}
%
В соответствии с рисунком гиперплоскости $H_4$ и $H_5$ не удовлетворяют неравенству $\left\langle \bm{a}_i,\bm{u} \right\rangle < b_i$ в~\eqref{eq:Q}. Гиперплоскости $H_3$ и~$H_6$ не удовлетворяют неравенству $\left\langle \bm{a}_i,\bm{d} \right\rangle > 0$ в~\eqref{eq:Q}. Таким образом, $\mathcal{Q} = \{1, 2\}$. Так как $\left\| \bm{u} - \bm{\gamma}_1(\bm{u},\bm{d}) \right\| < \left\| \bm{u} - \bm{\gamma}_2(\bm{u},\bm{d}) \right\|$, то ${\bm{\mu}}(\bm{u},\bm{d}) = \bm{\gamma}_1(\bm{u},\bm{d})$.

Сходимость алгоритма~\ref{alg:AlFaMove} к решению задачи ЛП~\eqref{eq:LP_problem} за конечное число итераций обеспечивается следующей теоремой.

\begin{theorem}{(Сходимость алгоритма AlFaMove)}\label{trm:alfamove_convergence} Пусть допустимый многогранник~$M$ задачи ЛП~\eqref{eq:LP_problem} является ограниченным непустым множеством. Пусть~$\bar{\bm{x}}$ --- решение задачи ЛП~\eqref{eq:LP_problem}. Тогда последовательность приближений~$\left\lbrace \bm{u}_k\right\rbrace_{k = 1}^K$, генерируемых алгоритмом~\ref{alg:AlFaMove}, является конечной ($K<+\infty$), и $\left\langle \bm{c}, {\bm{u}}_K \right\rangle = \left\langle \bm{c}, \bar{\bm{x}} \right\rangle$, то есть $\bm{u}_K$ является решением задачи ЛП~\eqref{eq:LP_problem}.
\end{theorem}
\begin{proof} 
	Обозначим через $\bm{d}_{AlFaMove}$ вектор $\bm{d}_{k+1}$, вычисляемый на шаге~8 алгоритма~\ref{alg:AlFaMove}. В соответствии с шагами~13,~14 алгоритма~\ref{alg:direction-vector} имеем
	\begin{equation*}%\label{eq:d_k}
		\bm{d}_{AlFaMove}=\bm{\rho}_{\mathcal{J}}(\bm{v}) - \bm{u}.
	\end{equation*}
	Согласно утверждению~\ref{prp:pseudoprojection=orthogonal_projection} отсюда следует
		\begin{equation}\label{eq:d_k}
		\bm{d}_{AlFaMove}=\bm{\pi}_{\mathcal{J}}(\bm{v}) - \bm{u},
	\end{equation}
	где $\bm{\pi}_{\mathcal{J}}(\bm{v})$ обозначает ортогональную проекцию точки $\bm{v}$ на линейное многообразие \mbox{$L=\bigcap\nolimits_{i \in \mathcal{J}}{H_i}$}. В соответствии с утверждением~\ref{prp:optimal_path_on_hyperplane} это означает, что вектор $\bm{d}_{AlFaMove}$ 	однозначно задает направление максимального увеличения целевой функции $f(\bm{x})$ задачи~\eqref{eq:LP_problem}. А это, в свою очередь, означает, что алгоритм~\ref{alg:AlFaMove} является численной реализацией метода поверхностного движения~\bblcite{Olkhovsky2023a}, то есть последовательности приближений $\left\lbrace \bm{u}_{AlFaMove}^k\right\rbrace$ и $\left\lbrace \bm{u}_{SMM}^k\right\rbrace$, генерируемые алгоритмом~\ref{alg:AlFaMove} из настоящей статьи и алгоритмом~1 из~\bblcite{Olkhovsky2023a} соответственно, совпадают. Таким образом сходимость алгоритма AlFaMove непосредственно следует из сходимости алгоритма поверхностного движения, гарантируемой теоремой~1 из статьи~\bblcite{Olkhovsky2023a}.
\end{proof}

\section{Параллельная версия алгоритма AlFaMove}\label{sec:parallel-version}

Наиболее ресурсоемкой операцией, используемой в алгоритме AlFaMove (алгоритм~\ref{alg:AlFaMove}), является операция вычисления вектора движения, выполняемая в цикле на шаге~8. При решении больших задач ЛП она занимает более 90\% процессорного времени. Это объясняется тем, что функция вычисления вектора движения $\bm{D}(\cdot)$, реализованная в виде алгоритма~\ref{alg:direction-vector}, использует в цикле на шаге~13 операцию псевдопроектирования, заключающуюся в последовательном применении отображения $\varphi(\cdot)$, задаваемого формулой~\eqref{eq:psi}, к исходной точке (см. алгоритм~\ref{alg:pseudoprojection}, реализующий вычисление псевдопроекции). Известно, что в случае больших задач ЛП проекционный метод может потребовать значительных временных затрат~\bblcite{Gould2008:Eng}. Кроме того, следует отметить, что алгоритм~\ref{alg:direction-vector} в цикле на шаге~12 перебирает все непустые подмножества множества $\mathcal{U}$, включающего в себя индексы гиперплоскостей, проходящих через точку $\bm{u}$. Если, например, через точку проходит 30 гиперплоскостей, то у нас получится $2^{30} -1= 1\,073\,741\,823$ непустых подмножества. Перебор такого количества подмножеств потребует суперкомпьютерных мощностей. Потому мы разработали параллельную версию алгоритма AlFaMove, представленную в виде алгоритма~\ref{alg:direction-vector_parallel}.

\begin{algorithm}[t]\caption{Параллельная версия алгоритма AlFaMove}\label{alg:direction-vector_parallel}\begin{multicols}{2}\begin{center}
						\textbf{мастер} \\
						\textbf{$l$-тый рабочий $(l=0,\ldots,L-1)$}
				\end{center}\end{multicols}\algrule\begin{multicols}{2}\begin{algorithmic}[1]
						\State \textbf{input} ${n,m,A,\bm{b},\bm{u}}_0$
						\State $k\ceq 0$
						\Repeat
						\State  \textbf{Broadcast} $\bm{u}_k$
						\State 
						\State 
						\State
						\State 
						\State 
						\State
						\State
						\State 
						\State 
						\State
						\State
						\State \textbf{Gather} $\mathcal{L}_{reduce}$
						\State $(\bm{d}_k,f_k) \ceq Reduce (\oplus,\mathcal{L}_{reduce}) $
						\If{$\bm{d}_k = \bm{0}$}
						\State $exit \ceq \textbf{true}$
						\Else
						\State ${\bm{u}}_{k+1} \ceq {\bm{\mu}}(\bm{u}_k,\bm{d}_k)$
						\State $k\ceq k+1$
						\State $exit \ceq \textbf{false}$
						\EndIf
						\State \textbf{Broadcast} $exit$
						\Until {$exit$}
						\State \textbf{output} $\bm{u}_k,f_k$
						\State \textbf{stop}
					\end{algorithmic}\begin{algorithmic}[1]
						%---------------------------------------------------------------
						\State \textbf{input} $n,m,A,\bm{b},\bm{c}$
						\State
						\Repeat
						\State \textbf{RecvFromMaster} $\bm{u}_k$
						\State $\mathcal{U} \ceq [\;]$
						\For{$i = 1 \dots m$}
						\If{$\langle  \bm{a}_i,\bm{u}_k \rangle = b_i$}
						\State $\mathcal{U} \ceq \mathcal{U} \doubleplus [i]$
						\EndIf
						\EndFor
						\State $K \ceq 2^{|\mathcal{U}|}-1$
						\State $L \ceq \mathrm{NumberOfWorkers}$
						\State $\mathcal{L}_{map(l)} \ceq [lK/L,\ldots,(l+1)K/L-1]$
						\State $\mathcal{L}_{reduce(l)} \ceq Map(\operatorname{F}_{\bm{u}_k}, \mathcal{L}_{map(l)})  $
						\State $(\bm{d}_l,f_l) \ceq Reduce (\oplus, \mathcal{L}_{reduce(l)}) $
						\State \textbf{SendToMaster} $(\bm{d}_l,f_l)$
						\State
						\State
						\State
						\State
						\State
						\State
						\State
						\State
						\State \textbf{RecvFromMaster} $exit$
						\Until{$exit$}
						\State
						\State \textbf{stop}
\end{algorithmic}\end{multicols}\end{algorithm}
%
Параллельный алгоритм~\ref{alg:direction-vector_parallel} построен на основе модели параллельных вычислений \sloppy BSF~\bblcite{Sokolinsky2021:Eng}, ориентированной на кластерные вычислительные системы. Модель BSF использует схему распараллеливания <<мастер--рабочие>> и требует представление алгоритма в виде операций над списками с использованием функций высшего порядка \emph{Map} и \emph{Reduce}.

В качестве второго параметра функции высшего порядка \emph{Map} в алгоритме~\ref{alg:direction-vector_parallel} используется список $\mathcal{L}_{map}=\left[1,\ldots, K\right]$, содержащий порядковые номера всех подмножеств множества~$\mathcal{U}$, за исключением пустого. Здесь $K=2^{|\mathcal{U}|}-1$. В качестве первого параметра фигурирует параметризованная функция
\[\operatorname{F}_{\bm{u}}:\left\lbrace 1,\ldots, K \right\rbrace  \to \mathbb{R}^n\times\mathbb{R}, \]
определенная следующим образом:
\begin{equation}\label{eq:F_x}
	\begin{gathered}
		\operatorname{F}_{\bm{u}}(j)=(\bm{d}_j,f_j);\hfill \\
		\bm{d}_j = 
		\begin{cases}
			\bm{e}_{\bm{d}}, & \text{if} \; (\bm{u} + \tau\bm{e}_{\bm{d}}) \in M \wedge  \langle  \bm{c}, \bm{w}\rangle > \langle \bm{c}, \bm{u} \rangle; \\
			\bm{0}, & \text{if} \; (\bm{u} + \tau\bm{e}_{\bm{d}}) \notin M \vee \langle  \bm{c}, \bm{w}\rangle \leqslant \langle \bm{c}, \bm{u} \rangle;
		\end{cases} \hfill \\
		f_j\, = 
		\begin{cases}
			\langle  \bm{c}, \bm{u} + \bm{e}_{\bm{d}}\rangle, & \text{if} \; (\bm{u} + \tau\bm{e}_{\bm{d}}) \in M \wedge  \langle  \bm{c}, \bm{w}\rangle > \langle \bm{c}, \bm{u} \rangle;\\
			-\infty, & \text{if} \; (\bm{u} + \tau\bm{e}_{\bm{d}}) \notin M \vee \langle  \bm{c}, \bm{w}\rangle \leqslant \langle \bm{c}, \bm{u} \rangle,
		\end{cases}
	\end{gathered}
\end{equation}
где
\begin{equation}\label{eq:w=rho...}
	\bm{w} = \bm{\rho}_{\sigma(j)}(\bm{u} + \delta\bm{c}/\| \bm{c} \|)
\end{equation}
и 
\[\bm{e}_{\bm{d}} = \frac{\bm{w} - \bm{u}}{\|\bm{w} - \bm{u}\|}.\]
Семантика функции $\operatorname{F}_{\bm{u}}(\cdot)$ однозначно определяется алгоритмом~\ref{alg:direction-vector}. Функция~$\sigma(\cdot)$, используемая в формуле~\eqref{eq:w=rho...}, отображает натуральное число $j \in \left\lbrace 1,\ldots, K \right\rbrace$ в $j$-тое подмножество множества элементов списка $\mathcal{U}$ следующим образом. Число~$j$ преобразуется в двоичное представление, состоящее из $|\mathcal{U}|$ разрядов. Каждому разряду соответствует индекс гиперплоскости из списка $\mathcal{U}$ в естественном порядке. Если разряд содержит единицу, то соответствующий индекс входит в подмножество~$\sigma(j)$, если --- ноль, то не входит. Например, пусть через точку $\bm{u}$ проходят гиперплоскости $H_2,H_4,H_7,H_9$. В этом случае $\mathcal{U} =\left[2,4,7,9 \right]$ и $K=2^4-1=15$, то есть из множества элементов списка~$\mathcal{U}$ может быть сформировано 15 различных непустых подмножеств.  Найдем, например, пятое подмножество. Функция $\sigma(\cdot)$ преобразует число~$5$ в двоичную запись из 4-x разрядов $0101$ и возвращает в качестве результата подмножество $\left\lbrace 4,9\right\rbrace$. Таким образом  функция высшего порядка $Map\left(\operatorname{F}_{\bm{u}}, \mathcal{L}_{map}\right)$ преобразует список $\mathcal{L}_{map}$ номеров подмножеств в список пар $(\bm{d}_j,f_j)$:
\begin{equation*}%\label{eq:key}
	Map\left(\operatorname{F}_{\bm{u}}, \mathcal{L}_{map}\right)=\left[\operatorname{F}_{\bm{u}}(1),\ldots,\operatorname{F}_{\bm{u}}\left(K\right) \right]=\left[(\bm{d}_1,f_1),\ldots,(\bm{d}_K,f_K) \right].
\end{equation*}
Здесь $\bm{d}_j$ ($j=1,\ldots, K$) является единичным вектором движения, а $f_j$ --- значение целевой функции, которое достигается в точке $\bm{u}+\bm{d}_j$. 

Обозначим через $\mathcal{L}_{reduce}$ список пар, генерируемый функцией высшего порядка $Map$:
\[\mathcal{L}_{reduce}=Map\left(\operatorname{F}_{\bm{u}}, \mathcal{L}_{map}\right) =\left[(\bm{d}_1,f_1),\ldots,(\bm{d}_K,f_K) \right].\]
Определим бинарную ассоциативную операцию
\[\oplus:\mathbb{R}^n\times\mathbb{R}\to \mathbb{R}^n\times\mathbb{R},\]
являющуюся первым параметром функции высшего порядка \emph{Reduce}:
\begin{equation}\label{eq:oplus}
	\left(\bm{d}',f'\right) \oplus \left(\bm{d}'',f''\right) =
	\begin{cases}
		\left(\bm{d}',f'\right),  & \mbox{если} f' \geqslant f'';\\
		\left(\bm{d}'',f''\right),& \mbox{если} f' < f''.
	\end{cases}
\end{equation}
Функция высшего порядка $Reduce$ редуцирует список пар $\mathcal{L}_{reduce}$ к одной паре путем последовательного применения операции $\oplus$ ко всем элементам списка: 
\begin{equation*}%\label{key}
	Reduce\left(\oplus, \mathcal{L}_{reduce}\right)=(\bm{d}_1,f_1)\oplus\ldots\oplus(\bm{d}_K,f_K) =(\bm{d}_{j'},f_{j'}),
\end{equation*}
где согласно~\eqref{eq:oplus}
\begin{equation*}
	j'= \arg \mathop{\max}\limits_{1\leqslant j\leqslant K} f_j.
\end{equation*}

Параллельная работа алгоритма~\ref{alg:direction-vector_parallel} организована по схеме <<мастер--рабочие>> и включает в себя $L+1$ процесс: один процесс--мастер и $L$ процессов--рабочих. Процесс--мастер осуществляет общее управление вычислениями, распределяет работу между процессами--рабочими, получает от них результаты и формирует итоговый результат. Для простоты будем предполагать, что количество подмножеств $K$ кратно количеству рабочих $L$. На шаге~1 мастер вводит исходные данные задачи ЛП и начальную точку~$\bm{u}_0$. На шаге~2 мастер устанавливает счетчик итераций $k$ в значение ноль. Шаги 3--26 реализуют основной цикл \textbf{repeat/until}, вычисляющий решение задачи ЛП~\eqref{eq:LP_problem}. На шаге~4 мастер рассылает текущее приближение $\bm{u}_k$ всем рабочим. На шаге~16 он получает от рабочих частичные результаты, которые на шаге~17 редуцируются в пару $(\bm{d}_k,f_k)$. Если на шаге~18 выполняется условие $\bm{d}_k= \bm{0}$, то решение найдено (мы предполагаем $\bm{d}_0 \ne \bm{0}$). В этом случае на шаге~19 мастер присваивает логической переменной $exit$ значение \textbf{true}. Если $\bm{d}_k \ne \bm{0}$, то мастер на шаге~21 вычисляет следующее приближение  ${\bm{u}}_{k+1}$,  увеличивает на шаге~22 счетчик итераций $k$ на единицу и на шаге~23 присваивает логической переменной $exit$ значение \textbf{false}. На шаге~25 мастер рассылает всем рабочим значение логической переменной $exit$. Если логическая переменная  $exit$ принимает значение <<истина>>, цикл \textbf{repeat/until} завершается на шаге~26. На шаге~27 мастер выводит в качестве результата последнее приближение ${\bm{u}}_k$ и оптимальное значение целевой функции $f_k$. Шаг~28 завершает работу процесса--мастера.

Все рабочие выполняют один и тот же код, но над различными данными. На шаге~1 $l$\nobreakdash-\hspace{0pt}тый рабочий \mbox{($l=1,\ldots,L$)} вводит исходные данные задачи ЛП. Цикл \textbf{repeat/until}  рабочего (шаги~3\nobreakdash-\hspace{0pt}26) соответствует циклу \textbf{repeat/until} мастера. На шаге~4 $l$\nobreakdash-\hspace{0pt}тый рабочий получает от мастера текущее приближение ${\bm{u}}_k$. Затем он формирует подсписок $\mathcal{L}_{map(l)}$ своих порядковых номеров подмножеств для последующей обработки (шаги 5--13). Подсписки различных рабочих не пересекаются:
\begin{equation}%\label{key}
	l'\ne l'' \Leftrightarrow \mathcal{L}_{map(l')} \ne \mathcal{L}_{map(l'')},
\end{equation}
а их конкатенация дает полный список:
\begin{equation}%\label{key}
	\mathcal{L}_{map}=\mathcal{L}_{map(0)}\doubleplus\ldots\doubleplus\mathcal{L}_{map(L-1)}.
\end{equation}
На шаге~14 рабочий вызывает функцию высшего порядка $Map$, которая, в свою очередь, применяет параметризованную функцию  $\operatorname{F}_{\bm{u}_k}$, определенную формулами~\eqref{eq:F_x}, ко всем элементам подсписка $\mathcal{L}_{map(l)}$, формируя на выходе подсписок пар $\mathcal{L}_{reduce(l)}$. Этот подсписок на шаге~15 редуцируется рабочим в единственную пару  $(\bm{d}_l,f_l)$ с помощью функции высшего порядка $Reduce$, которая последовательно применяет бинарную операцию $\oplus$, определенную по формуле~\eqref{eq:oplus}, ко всем элементам подсписка  $\mathcal{L}_{reduce(l)}$. Результат пересылается мастеру на шаге~16. На шаге~25 рабочий получает от мастера значение логической переменной $exit$. Если эта переменная принимает значение \textbf{true}, то рабочий процесс завершается. В противном случае продолжает выполняться цикл \textbf{repeat/until}. Операторы обмена \textbf{Broadcast}, \textbf{Gather}, \textbf{RecvFromMaster} и \textbf{SendToMaster} обеспечивают неявную синхронизацию работы процесса--мастера и процессов--рабочих.

\section{Вычислительные эксперименты}\label{sec:implementation_end_experiments}

Мы реализовали параллельную версию алгоритма AlFaMove  на языке C++ с использованием программного BSF-каркаса~\bblcite{Sokolinsky2021a:Eng}, базирующегося на модели параллельных вычислений BSF~\bblcite{Sokolinsky2021:Eng}. BSF-каркас инкапсулирует все аспекты, связанные с распараллеливанием программы на основе библиотеки MPI. Исходные коды параллельной реализации алгоритма AlFaMove решения задач ЛП свободно доступны в репозитории GitHub по адресу \url{https://github.com/leonid-sokolinsky/AlFaMove}. Разработанная программа была протестирована на большом количестве задач ЛП, взятых из различных источников. Все эти задачи в формате MTX~\bblcite{Boisvert1996:Eng} доступны по адресу \url{https://github.com/leonid-sokolinsky/Set-of-LP-Problems}. В качестве тестов мы также использовали искусственные задачи, полученные с помощью генератора случайных задач ЛП FRaGenLP~\bblcite{Sokolinsky2021b}. Эти задачи доступны по адресу \url{https://github.com/leonid-sokolinsky/Set-of-LP-Problems/tree/main/Rnd-LP}. Мы не смогли протестировать реализацию AlFaMove на задачах из репозитория \mbox{Netlib-LP}~\bblcite{Gay1985:Eng}, так как во всех этих задачах количество гиперплоскостей, проходящих через начальную точку $\bm{u}_0$, превышало число 30, что соответствует количеству возможных комбинаций, равному $1\,073\,741\,824$. Массивы таких размеров не допускаются доступными нам компиляторами C++.

С~помощью разработанной программы мы исследовали масштабируемость алгоритма AlFaMove. В экспериментах мы использовали параметризованную задачу ЛП <<гиперкуб с отсеченной вершиной>>, для которой размерность пространства $n$ является параметром. Ограничения этой задачи содержат $2n+1$ неравенств следующего вида:
\begin{equation}\label{sokol_Formula1}
	\left\{ {\begin{array}{*{20}{l}}
			{{x_1}}&{}&{}&{}& \leqslant & 200  \\
			{}&{{x_2}}&{}&{}& \leqslant & 200  \\
			{}&{}& &{}& ~\vdots \\
			{}&{}&{}&{{x_n}}& \leqslant & 200  \\
			{{x_1~~~+}}&{ {x_2}}& \cdots &{ +~~~~{x_n}}& \leqslant &{200(n - 1) + 100}\\
			{{x_1} \geqslant 0,}&{{x_2} \geqslant 0,}& \cdots &, {{x_n} \geqslant 0}.
	\end{array}} \right.
\end{equation}
Градиент целевой функции задается вектором
\begin{equation}\label{sokol_Formula3}
	\bm{c} = \left( {1,2, \ldots ,n} \right).
\end{equation}
Задача предполагает нахождение максимума целевой функции и имеет единственное решение в точке $(100,200,\ldots,200)$ со значением целевой функции, равным  $100(n^2+n-1)$. Для произвольного $n$ эта задача может быть получена в формате MTX с помощью генератора FRaGenLP, если в качестве количества случайных неравенств задать $0$. Эти задачи ЛП для различных $n$ доступны по адресу \url{https://github.com/leonid-sokolinsky/Set-of-LP-Problems/tree/main/Rnd-LP} под именами \texttt{lp\_rnd<n>-0}, где в качестве \texttt{<n>} указана размерность пространства. Вычислительные эксперименты проводились на суперкомпьютере <<Торнадо \mbox{ЮУрГУ}>>~\bblcite{Dolganina2022:Eng}, характеристики которого представлены в табл.~\ref{tbl:Tornado_SUSU}. 
%	
\begin{table}[t]
	\caption{Характеристики кластера <<Торнадо ЮУрГУ>>}
	{\begin{tabular}{|l|l|}
			\hline
			Параметр & Значение \\
			\hline
			Количество процессорных узлов & 480 \\
			Процессоры & Intel Xeon X5680 (6 cores, 3.33 GHz) \\
			Число процессоров на узел & 2\\
			Память на узел & 24 GB DDR3\\
			Соединительная сеть & InfiniBand QDR (40 Gbit/s) \\
			Операционная система & Linux CentOS\\
			\hline
		\end{tabular}\label{tbl:Tornado_SUSU}}
\end{table}
%
Все вычисления выполнялись с двойной точностью, при которой число с плавающей точкой занимает в оперативной памяти 64 бита.

В первой серии экспериментов исследовалась зависимость ускорения и эффективности распараллеливания алгоритма AlFaMove от количества используемых процессорных узлов кластера при решении задач <<гиперкуб с отсеченной вершиной>>. Результаты этих экспериментов представлены на рис.~\ref{fig:speedup_diagramm}.
%
\begin{figure}[b]
	\includegraphics[height=6 cm]{Fig_Diagrams}
	\caption{Ускорение и эффективность распараллеливания алгоритма AlFaMove.}
	\label{fig:speedup_diagramm}
\end{figure}
%
В данном контексте ускорение $\alpha(L)$ определялось как отношение времени~$T(1)$ решения задачи на конфигурации с узлом-мастером и единственным узлом-рабочим ко времени~$T(L)$ решения той же задачи на конфигурации с узлом-мастером и $L$ узлами-рабочими:
\begin{equation*}%\label{eq:}
	\alpha(L)=\frac{T(1)}{T(L)}.
\end{equation*}
Эффективность распараллеливания $\beta(L)$ вычислялась по формуле
\begin{equation*}%\label{eq:}
	\beta(L)=\frac{T(1)}{L \cdot T(L)}.
\end{equation*}
Вычисления проводились для размерностей $16$, $18$ и $20$.  Число ограничений соответственно составило $33$, $37$ и $41$. В качестве начальной точки всегда выбиралась вершина допустимого многогранника со следующими координатами:
\begin{equation}\label{eq:start_vertex}
	x_1=0,\quad\ldots\quad x_{n/2}=0,\quad x_{n/2+1}=200,\quad\ldots\quad x_n=200.
\end{equation}
%
Эксперименты продемонстрировали хорошую масштабируемость алгоритма AlFaMove на задачах <<гиперкуб с отсеченной вершиной>>, начиная с размерности $n=18$. В этом случае алгоритм демонстрировал ускорение, близкое к линейному. На меньших размерностях затраты на обмены и латентность начинают превалировать над вычислительными затратами, что приводит к резкому уменьшению границы масштабируемости алгоритма\footnote{Под границей масштабируемости понимается максимальное количество процессорных узлов, на котором наблюдается рост ускорения.}. Для $n=16$ эта граница оказалась равной 180 узлам. Эксперименты также показали, что с увеличением размерности задачи эффективность распараллеливания на малом количестве процессорных узлов (меньше 120) падает. Однако, при большем количестве процессорных узлов наблюдается обратная тенденция. Так, для размерности $n=16$ эффективность распараллеливания на 20 процессорных узлах составила 61\%, после чего снизилась до 23\% на 220 узлах. При этом, для $n=20$ эффективность распараллеливания оказалась равной 50\% и 40\% соответственно.

В следующей серии экспериментов исследовалась зависимость времени вычислений от размерности задачи <<гиперкуб с отсеченной вершиной>> на различных многопроцессорных конфигурациях с количеством рабочих узлов $L=60$, $L=120$ и $L=180$. Результаты этих экспериментов представлены на рис.~\ref{fig:time}.
%
\begin{figure}[t]
	\includegraphics[height=6 cm]{Fig_Time}
	\caption{Зависимость времени работы AlFaMove от размерности задачи на различных многопроцессорных конфигурациях (\textit{L} --- количество рабочих узлов).}
	\label{fig:time}
\end{figure}
%
Размерность варьировалась от $n=16$  до $n=24$ с шагом 2. Для размерности $n=24$ каждый из списков $\mathcal{L}_{map}$ и $\mathcal{L}_{reduce}$ включал в себя~$16\,777\,215$ элементов. Это --- максимальный размер, допускаемый используемым компилятором. В качестве начальной точки всегда выбиралась вершина допустимого многогранника с координатами~\eqref{eq:start_vertex}. На всех исследуемых конфигурациях эксперименты показали экспоненциальный рост времени вычислений с увеличением размерности задачи. Однако, конфигурации с б\'{о}льшим количеством рабочих узлов демонстрировали существенно меньшее время работы алгоритма AlFaMove.

В третьей серии экспериментов мы исследовали поведение алгоритма AlFaMove на кубе Кле-Минти (Klee-Minty), представляющем собой параметризованную задачу ЛП с размерностью пространства в качестве параметра. Область допустимых решений этой задачи представляет собой гиперкуб с перекошенными углами и может быть задана следующей системой неравенств:
\begin{equation*}%\label{eq:klee-minty}
	\left\{ {\begin{array}{*{20}{l}}
			{{x_1}}&{}&{}&{}& & \leqslant &  5  \\
			{4x_1~~+}&{{x_2}}&{}&{}& &\leqslant & 25  \\
			{8x_1~~+}&{4x_2~~+}& x_3 &{}& & \leqslant & 125 \\
			& & & & & ~\vdots \\
			{{2^n x_1~+}}&{ {2^{n-1} x_2}}& \cdots &{ +~~~4{x_{n-1}}} &{ +~~~~{x_n}}& \leqslant & 5^n\\
			{{x_1} \geqslant 0,}&{{x_2} \geqslant 0,}& \cdots &, {{x_n} \geqslant 0}.
	\end{array}} \right.
\end{equation*}
Градиент целевой функции задается вектором
\begin{equation*}%\label{eq:}
	\bm{c} = \left( {2^{n-1},\,2^{n-2},\, \ldots \,,2,\,1} \right).
\end{equation*}
Задача предполагает нахождение максимума целевой функции и имеет единственное решение в точке $(0,\,\ldots\,,0,\,5^n)$ со значением целевой функции, равным  $5^n$. Кле и Минти в статье~\bblcite{Klee1972:Eng} показали, что в случае, когда начальной вершиной является центр координат, классический симплекс-метод при решении этой задачи посещает все $2^n$ вершин, делая $2^n-1$ итераций. Известно, что большинство оптимизационных алгоритмов демонстрируют плохую производительность применительно к кубу Кле-Минти. Мы применили алгоритм AlFaMove для решения кубов Кле-Минти размерности от 5 до 9. Результаты экспериментов, приведенные в табл.~\ref{tbl:Klee-Minty}, показывают, что алгоритм AlFaMove во всех случаях находил решение за $2n-1$ итераций, в то время, как симплекс-метод совершал $2^n-1$ итераций.\linebreak
%
\begin{table}[H]
	\caption{Эксперименты с кубами Кле-Минти.}
	{\begin{tabular}{|>{\centering}p{0.1\textwidth}|>{\centering}p{0.11\textwidth}|>{\centering}p{0.07\textwidth}|>{\centering}p{0.12\textwidth}|>{\centering}p{0.11\textwidth}|p{1.7 cm}|}
			\hline
			Размер-  & \multicolumn{4}{c|}{AlFaMove} & \, Simplex \\
			\cline{2-6}
			ность & Граница масшт-ти & Время (сек.) & Относит. погреш-ть & Кол-во итераций & \phantom{k}Кол-во итераций \\
			\hline
			5 & 10	& 0.2	& $0.9 \cdot 10^{-12}$	& 9		& ~~~~~31 \\
			6 & 15 	& 2 	& $0.2 \cdot 10^{-12}$	& 11	& ~~~~~63 \\
			7 & 20 	& 13 	& $0.8 \cdot 10^{-11}$ 	& 13	& \,~~~~127 \\
			8 & 25 	& 126 	& $0.8 \cdot 10^{-11}$	& 15	& \,~~~~255 \\
			9 & 30 	& 1445 	& $0.2 \cdot 10^{-10}$	& 17 	& \,~~~~511 \\
			\hline		\end{tabular}\label{tbl:Klee-Minty}}
\end{table}
%
\noindent Относительная погрешность вычислялась по формуле
\begin{equation*}
	\delta = \left| \frac{f_{exact}-f_{approx}}{f_{exact}}\right| ,
\end{equation*}
где $f_{exact}$ --- точное максимальное значение целевой функции, $f_{approx}$ --- значение, вычисленное алгоритмом AlFaMove. Подсчет итераций симплекс-метода производился с помощью онлайн калькулятора, доступного по адресу \url{https://www.pmcalculators.com/simplex-method-calculator}. Эксперименты также показали, что порог масштабируемости алгоритм AlFaMove на кубах Кле-Минти рос линейно с ростом размерности. При этом одновременно наблюдался экспоненциальный рост времени вычислений. 

\section*{Заключение}

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

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

Алгоритмы проекционного типа характеризуются низкой скоростью сходимости, зависящей от углов между гиперплоскостями, образующими линейное многообразие. Также отмечено, что при вычислении вектора движения возникает переборная задача комбинаторного типа, имеющая высокую пространственную и временную сложность. Представлена параллельная версия алгоритма AlFaMove, ориентированная на кластерные вычислительные системы. Параллельная версия реализована на языке C++ с использованием программного BSF-каркаса, основанного на модели параллельных вычислений BSF. Проведены эксперименты по исследованию масштабируемости алгоритма AlFaMove на кластерной вычислительной системе. Вычислительные эксперименты показали, что задача линейного программирования с $24$ переменными и $49$ ограничениями демонстрирует ускорение близкое к линейному на $320$ процессорных узлах кластера. Задачи большей размерности приводили к ошибке компилятора, связанной с превышением максимального допустимого размера массивов.

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

\vspace{1em}
{\it Исследование выполнено при финансовой поддержке РНФ (проект №\,23"~21"~00356).}

\insertbiblio

% Сведения об авторах (Ф.И.О.(полностью), степень, звание, кафедра (отдел), вуз, e-mail)
Ольховский Николай Александрович, аспирант, кафедра системного программирования, Южно-Уральский государственный университет (национальный исследовательский университет) (Челябинск, Российская Федерация)

Соколинский Леонид Борисович, д.ф.-м.н., профессор, заведующий кафедрой системного программирования, Южно-Уральский государственный университет (национальный исследовательский университет) (Челябинск, Российская Федерация)

{\phantom{eee} \hrule \vskip 4 pt}
\classify{} % Код MSC не нужен
%\classify{MSC 00000} % Коды MSC см. http://www.ams.org/msc/

\title{NUMERICAL IMPLEMENTATION OF SURFACE MOVEMENT METHOD IN LINEAR PROGRAMMING}

\authorsen{L.B.~Sokolinsky, N.A.~Olkhovsky, I.M.~Sokolinskaya}

\organizations{South Ural State University (pr. Lenina 76, Chelyabinsk, 454080 Russia)}

\emails{\emailsinsert}
\redactiondateen{\rdate}

\begin{abstract}%
The article is devoted to the numerical implementation of a new linear programming method called the surface movement method. The implementation is based on the AlFaMove algorithm, which builds, on the surface of the feasible polytope, the optimal objective path from an arbitrary boundary point to a point that is a solution to the linear programming problem. The optimality of the path lies in the fact that the direction of movement along the face of the polytope corresponds to the maximum increase in the value of the objective function. To calculate the optimal direction of movement, a method based on the operation  of pseudoprojection onto a linear manifold is used. The pseudoprojection operation is a generalization of the notion of orthogonal projection. Pseudo projection is implemented using an iterative projection-type algorithm. The proposition is proved that the pseudoprojection coincides with the orthogonal projection in the case of a linear manifold that is the intersection of hyperplanes. It is also proved that, in the case of a linear manifold, the pseudoprojection method calculates a movement vector that coincides with the direction of maximum increase of the objective function. A parallel implementation of the AlFaMove algorithm is performed. The results of computational experiments on a cluster computing system are presented, which demonstrate the high scalability of the proposed numerical implementation.

\keywordsen{linear programming, surface movement method, numerical implementation, AlFaMove algorithm, parallel implementation, cluster computing system, scalability evaluation.}
\end{abstract}

\forcitationen{Sokolinsky~L.B., Olkhovsky~N.A., Sokolinskaya~I.M. Numerical Implementation of Surface Movement Method in Linear Programming}

\openaccess

\insertbiblioen
% В Bibliography_sec.bib
% 1) заменить все <<:Eng,>> на <<:Eng:TR,>>;
% 2) заменить все <<:Rus,>> на <<:Rus:TR,note = {(in Russian)},>>.

%Англоязычная публикация~\bblcite{Keil2006:Eng} должна иметь citation key вида Keil2006:Eng (<фамилия 1-го автора><год публикации>:Eng). Помещается в оба списка.

%Русскоязычная публикация~\bblcite{Zorkaltsev2018}, имеющая англоязычную версию должна иметь citation key вида Zorkaltsev2018 (<фамилия 1-го автора><год публикации>). Помещается в первый список. Англоязычная версия должна иметь citation key вида Zorkaltsev2018:TR (<фамилия 1-го автора><год публикации>:TR). Помещается во второй список.

%Русскоязычная публикация~\bblcite{Eremin1979:Rus}, не имеющая англоязычной версии должна иметь citation key вида Eremin1979:Rus (<фамилия 1-го автора><год публикации>:Rus). Помещается в первый список. Англоязычная транскрипция ссылки на такую публикацию должна иметь такой же citation Eremin1979:Rus (<фамилия 1-го автора><год публикации>:Rus). Помещается во второй список.


\end{document}
