ПРОГРАММНАЯ РЕАЛИЗАЦИЯ ВОССТАНОВЛЕНИЯ ВЫПУКЛОГО МНОЖЕСТВА С ПОСТОЯННЫМ УГЛОМ ПОВОРОТА

Рассмотрим и реализуем алгоритм восстановления выпуклого множества по его опорной функции с помощью программных средств Anaconda (дистрибутив Python) и GNU Octave. Данный алгоритм должен реализовываться с заданной точностью, так как вычислительные погрешности решения вспомогательных подзадач алгоритма могут оказаться весьма значительными.

Исследуем основную задачу, по которой строится алгоритм восстановления выпуклого множества. Имеется опорная функция c(F,ψ) непустого выпуклого компактного множества , которая находится в пространстве  и задается произвольное положительное число ε. Алгоритм реализует построение множество  так, чтобы выполнялось включение  и хаусдорфово расстояние между множествами  и  не превосходило числа ε.

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

Теорема 1. Пусть множество  и для углов  и  выполняется неравенства:

тогда справедливы оценка h (,  и включение

 – пространство, которое состоит из всех выпуклых компактных подмножеств пространства ;

 – базовый угол поворота;

 – базовый угол поворота, для размерности ;

Лемма 1. Пусть множество  и для углов  выполняется равенство:

Тогда справедливы оценка h (,  и включение .  Данная Лемма 1. – это частный случай Теоремы 1 для размерности .

Для построения множества  в пространстве  вводим сферическую систему координат:

, при

Сферический радиус при этом будет , углы сферы – ,

 , при . [1]

Через данную систему координат выразим вектор единичный длины :

Теперь через единичный вектор определим , имеющий угол , где .

Для ознакомления потребуются общий вид координат точки, через опорную функцию и вектор . [2]

где  = .

Для  точка будет иметь координаты:.

Определяем выпуклое множество, рассматриваемое в данном алгоритме,

 = conv { | . [3]

Основные этапы восстановления выпуклого множества:

• Задать начальное приближение для угла/углов  / и ;
• Итеративно уменьшать углы с некоторым мелким шагом, чтобы выполнялись условия Леммы 1;
• Основной цикл: вычислять координаты очередной точки выпуклой оболочки с помощью численного нахождения максимума опорной функции и описанных формул.

Основная программа написана на Python версии 3.6 в среде Anaconda 3, использует стандартные для научных расчётов библиотеки numpy, scipy, matplotlib. Благодаря numpy появляется возможность векторизовать вычисления, это существенно ускоряет работу алгоритма, если размерность n пространства достаточно большая. [4]

Библиотека scipy содержит множество реализаций методов оптимизации первого и второго порядков и в данном случае применяется для максимизации опорной функции на выпуклом множестве (эллипсе/эллипсоиде), что необходимо при вычислении координат точек выпуклой оболочки. Библиотеку matplotlib можно использовать для отрисовки результатов работы программы в случае n=2, однако для визуализации была выбрана более удобная среда GNU Octave 4.4. Octave является свободным программным обеспечением и призван стать альтернативой среде MATLAB. [5]

Язык среды Octave также интерпретируемый, а большинство команд совместимы с её коммерческим аналогом. Octave нативно поддерживает, в отличие от сторонних библиотек Python, целый ряд методов для работы с 2D и 3D графикой, а кроме того, как и MATLAB, интерактивные инструменты визуализации.

Для удобства было создано 3 файла с расширением «.ipynb» (IPython Notebook) : основной, где описаны функции, универсальные в использовании – подходящие  к восстановлению множеств размерности  и множеств размерности , и отдельные файлы с функциями, использующиеся только в каждой из размерностей.

Опишем основные функции, которые имеются в программе.

Def make_ellipsoid(axis_coefficients), где axis_coefficients – массив длины 2 или 3 с коэффициентами эллипса/эллипсоида.  Функция возвращает функцию, которая в заданной точке формирует   – return lambda x: (x**2).dot(1.0 / (axis_coefficients**2)) - 1.

Функция с параметрами - хаусдорфово расстояние до 0, точность ε и размерностью, которая будет 2 или 3. Данная функция возвращает функции, которые проверяют, удовлетворяет ли заданные углы  и  условиям теоремы 1. Если размерность будет равна 2, то  и  не используется, его проверять не нужно.

Функция, которая проверяет угол по условию Леммы 1.: check_lemma_alpha = lambda a: a <= asin(0.25*tol/(area+0.5*tol).

sqrt_ndim = np.sqrt(ndim) # ndim - размерность,  или   в нашем случае, берем квадратный корень.

check_lemma_gamma = lambda g: g<=2*asin(tol/(6*sqrt_ndim*(area+0.5*tol))) if ndim >= 3 else None – аналогично, как и в случае check lemma alpha. Если , в эту переменную ничего не записывается. Последнее – возвращаем пару значений, а именно, массив из 2 элементов (функций) return check_lemma_alpha, check_lemma_gamma.

Используем функции def psi2d(angle) и def psi3d(angles) для высчитывания координат для точки, где Angle – скаляр, а Angles – массив из 2 элементов.

Функция, которая находит максимум опорной функции (sup_func), с помощью стандартного метода minimize из библиотеки scipy.

Minimize численно решает задачу условной минимизации с ограничением–равенством (constraints, type='eq'). Другими словами, просто находим максимум опорной функции на эллипсе/эллипсоиде.

Def maximize_sup_func(psi, ellipsoid_func, x0), где

psi – массив точек, которые считаются с помощью psi2d или psi3d;

ellipsoid_func – функция, возвращаемая make_ellipsoid;

x0 – некоторая точка, с который начнётся оптимизационный процесс.

Функция, которая считает скалярное произведение (psi, x) в произвольной x. Она же опорная функция – sup_func = lambda x: -psi.dot(x), максимизируем опорную функцию на эллипсе/эллипсоиды, возвращаем найденный максимум

return -minimize(sup_func, x0, constraints={'type': 'eq', 'fun': ellipsoid_func})['fun'].

Функция, которая будет считать по формулам координаты точки def cchp2d(cfk, cfk1, angles, a), здесь использовалась cchp – аббревиатура Calculate Convex Hull Point.

Функция которая будет возвращать массив длины 2

arr = zeros_(2)

arr[0] = (cfk1 * cos_(angles[0]) - cfk * cos_(angles[1])) / sin_(a)

arr[1] = (cfk * sin_(angles[1]) - cfk1 * sin_(angles[0])) / sin_(a)

return arr

Для массива длины 3 функция будет записываться аналогично.

Далее кратко опишем файл для размерности , для этого инициализируем переменные: argv[1] – точность, argv[2],[3] - растяжения вдоль осей эллипсоид, argv[4] - файл, в которой сохранится файл с выпуклой оболочкой eps = float(sys.argv[1]), ax_coefs = np.array([float(sys.argv[2]), float(sys.argv[3])]), output_file = sys.argv[4], x0 = np.copy(ax_coefs), ellip_func = make_ellipsoid(ax_coefs) и добавим некоторое начальное значение угла step = 1e-4, соответственно, alpha = np.pi / 2 – step.

Затем уменьшаем угол, пока не выполнится условия Теоремы 1:

check_alp, _ = make_alpha_gamma_checker(ax_coefs.max(), eps, ndim=2)  while not check_alp(alpha): alpha - = step print('alpha_valid={}'.format(alpha)).

Построим сетку углов, её размер = количеству точек оболочки: grid_alpha = np.arange(int(2 * np.pi / alpha) + 1); print(len(grid_alpha)). Массив размера (количество точек оболочки)х(2) points = np.zeros((len(grid_alpha), 2), dtype='float32'); angles = np.array([0.0, alpha]).

Основные вычисления будут следующими:

for k in grid_alpha:

cfk = maximize_sup_func(psi2d(k * alpha), ellip_func, x0) + 0.5*eps

cfk1 = maximize_sup_func(psi2d((k+1) * alpha), ellip_func, x0) + 0.5*eps

 points[k] = cchp2d(cfk, cfk1, [k * alpha, (k+1) * alpha], alpha)

Файл с оболочкой сохраняется np.savetxt(output_file, points, fmt='%.10f').

Для множества с размерностью  будет аналогично, только добавится gamma = np.pi / 2 – steр:

 points = np.zeros((len(grid_gamma) * len(grid_alpha), 3), dtype='float32')

for i, a1 in enumerate(grid_gamma):

angles = np.array([a1, 0.0, alpha])

for j, a2 in enumerate(grid_alpha):

cfk = maximize_sup_func(psi3d([angles[0],angles[1]]), ellip_func, x0) + 0.5*eps

cfk1 = maximize_sup_func(psi3d([angles[0],angles[2]]), ellip_func, x0) + 0.5*eps

 points[i*len(grid_alpha)+j] = cchp3d(cfk, cfk1, angles, alpha)

angles[-2:] += alpha

np.savetxt(output_file, points, fmt='%.10f')

Чтобы полученные значения выпуклых оболочек визуализировать, воспользуемся средой Octave.

Нарисуем основные для обоих случаев параметры:

function [x, y] = ellipse(a, b, n); angle = linspace(0, 2*pi, n).'; x = a * cos(angle);y = b * sin(angle).

Далее на примере размерности  визуализируем выпуклую оболочку для множества. Рисунок 2 эллипсоид с a=3,b=2,c=3, epsilon=1, в сетке 5162 точки (58 по альфа х 89 по гамма), угол альфа=0.071, гамма=0.055.

ch = importdata("tmp_new/ch");

fig = get(0,"currentfigure");

clf;

ax = gca;

hold(ax);

[x,y,z] = ellipsoid(0,0,0,3,2,2,200);

surf(ax,x,y,z,"edgecolor","none");

j = convhulln(ch);

srf = trisurf(j,ch(:,1),ch(:,2),ch(:,3));

set(srf,"facecolor","y","edgecolor","k")

p1 = max(ch(:));

p2 = min(ch(:));

p = max(p1,abs(p2));

shift = 0.03;

axis([-(1+shift)*p (1+shift)*p ...

  -(1+shift)*p (1+shift)*p -(1+shift)*p (1+shift)*p],"equal")

set(ax,"fontsize",12)

view(ax,3)

zoom(1.5)

grid on

Для размерности  будет записано аналогично и проще. Рисунок 1 эллипс с a=2,b=10, epsilon=1, точек 265, альфа=0.024.

Рисунок 1. Эллипс с a=2, b=10, epsilon=1, точек 265, альфа=0.024

Рисунок 2. Эллипсоид с a=3, b=2, c=3, epsilon=1, в сетке 5162 точки (58 по альфа х 89 по гамма), угол альфа=0.071, гамма=0.055

Список литературы:

1. Васильев Ф.П. Численные методы решения экстремальных задач: учеб. пособие для вузов. 2-е изд., перераб. и доп. М.: Наука. Гл. ред. физ.-мат. лит., 1988, С.148-162.
2. Воеводин В.В. Линейная алгебра. 2-е изд., перераб. и доп. М.: Наука. Гл. ред. физ.-мат. лит., 1980, С.33-40,79-99.
3. Понтрягин Л.С., Болтянский В.Г., Гамкрелидзе Р.В., Мищенко Е.Ф. Математическая теория оптимальных процессов. – М.: Наука. Гл. ред. физ.-мат. лит., 1961. 392 с.
4. GNU Octave [Электронный ресурс]: сайт. - URL: https://www.gnu.org/software/octave/ (дата обращения: 23.03.19)
5. Anaconda [Электронный ресурс]: сайт. - URL: https://www.anaconda.com/distribution/ (дата обращения: 23.03.19)