Максиминное расположение точек на сфере

[B@R5uk]

Пусть имеется N точек. Требуется расположить их на единичной сфере так, чтобы минимальное расстояние от каждой точки до любой другой было максимально. Какова будет симметрия получающегося многогранника? Какого будет количество рёбер минимальной длины (жёстких рёбер)? Какова будет длина такого ребра? Кто-нибудь интересовался такой постановкой задачи?

Сразу скажу, это НЕ задача Томсона. Если уж и искать физический аналог, который будет расталкивать точки по сфере, то им будет набор упругих стержней. И то, не между любой парой точек, а только вдоль "жёстких" ребер. Тем не менее, в некоторых простых случаях, ответы на эти обе задачи совпадают. Например, тетраэдр, октаэдр и икосаэдр являются решениями. Но наибольший интерес, на мой взгляд, представляют промежуточные случаи.

[Утундрий]

Расстояния меняем по сфере или в объемлющем пространстве?

[B@R5uk]

Да без разницы особо.

$$d=2\sin\frac{\varphi}{2}$$

Жёсткие рёбра, которые задают форму, будут все равны. Тут есть одно исключение для случая 5 точек. Оптимальным решением является фигура с ребром корень из двух (или пи пополам, если в угловой мере), где две точки расположены на противоположных сторонах сферы (назовём их полюсами), а три другие болтаются на экваторе. Формулировки задачи недостаточно, чтобы их зафиксировать в однозначную фигуру, например, две склеенные треугольные пирамиды. Но это единственное такое исключение, на сколько мне известно.

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

[Утундрий]

Разница возникнет, если длина ребра-дуги больше половины периметра большого круга.

[мат-ламер]

B@R5uk
Посмотрите задачу 120 из книги Шклярский и др. "Геометрические оценки ...". Но эта задача непростая. В книге есть картинки для $n<10$ . Наверное сейчас с помощью компьютера решили (хотя бы приблизительно) эту задачу и для больших $n$ . Возможно тут на форуме что-то похожее и было. Попробуйте поискать.

[Alex Krylov]

Во книжоночка интересная попалась Thomas Ericson, Victor Zinoviev, Codes on Euclidean Spheres и там 3 глава Codes in dimension n=3

[Sender]

Гугл говорит, что это проблема Таммеса. Любопытно, что он был не математиком, а ботаником, исследовавшим пыльцевые зёрна.
https://pdodds.w3.uvm.edu/teaching/courses/2009-08UVM-300/docs/others/1994/mooers1994a.pdf

[Dmitriy40]

Разве это не Задача Томсона?
Примерно по ней же была тема «Экстремальное расположение зарядов на сфере», в которой полезные ссылки:

мат-ламер в сообщении #1584244 писал(а):

На всякий случай приведу ссылку на журнал "Математическое просвещение" https://www.mccme.ru/free-books/matpros2.html. Смотрите первую статью в разделе "Наш семинар".

мат-ламер в сообщении #1584301 писал(а):

Пока читаю следующие статьи (вышел через статьи из Википедии про задачу Томсона): статья 1 и статья 2 . Вторая статья вскоре будет доступна и на русском в ЖВМиМФ.

manul91 в сообщении #1584326 писал(а):

Min-Energy Configurations of Electrons On A Sphere
Уже для 16-ти зарядов две разные равновесные конфигурации.

[B@R5uk]

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

Dmitriy40, я же в первом посте написал, что это не она. Вон выше тоже мне указали, что задача носит название проблемы Таммеса. Ради интереса я сравнил решения, и уже для семи точек они отличаются принципиально: проблема Таммеса даёт фигуру на рисунке ниже слева, а задача Томсона в качестве решения имеет пятиугольную бипирамиду (рисунок в центре). Для восьми точек проблема Таммеса и задача Томсона приводят к похожим решениям: квадратной антипризме. Но в случае первой — это честный полуправильный многогранник (рисунок справа), а случае второй — он вытянут по вертикали (треугольные грани являются не равносторонними, а равнобедренными треугольниками с углом при основании 71,694°).

мат-ламер в сообщении #1680918 писал(а):

Но эта задача непростая. В книге есть картинки для $n<10$ . Наверное сейчас с помощью компьютера решили (хотя бы приблизительно) эту задачу и для больших $n$ .

Я так понимаю, в задачнике делался упор на доказательство оптимальности? Потому что там только до шести точек включительно, а самое интересное начинается от семи и более. Я тоже эту задачу пытаюсь численно решать. Сначала отжигом со случайными блужданиями с симуляцией расталкивания каждой точки от ближайшей соседней, чтобы приблизительно прикинуть топологию минимальных рёбер, а потом уже численно точно Нелдером-Мидом подгоняю значения сферических координат, чтобы все жёсткие рёбра были одинаковы. Но отжиг надёжно сходится только для небольшого количества точек, меньше двух десятков.

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

[B@R5uk]

B@R5uk в сообщении #1680963 писал(а):

жёсткие рёбра (минимальной длины) принадлежат выпуклой оболочке вершин.

Если обозначить длину ребра l, а расстояние от центра единичной сферы до этого ребра d, то будет выполняться

$$\frac{l^2}{4}+d^2=1$$

Это означает, если длина ребра минимально возможная, то расстояние максимально, что с необходимостью приводит к тому, что ребро принадлежит выпуклой оболочке. Потому что в противном случае найдётся пара точек на сфере, расстояние между которыми меньше величины l, что будет противоречить требованию её минимальности.

Тут у меня возник другой вопрос: есть ли какой-нибудь алгоритмизируемый способ убедиться, что заданная конструкция жёстких рёбер сама по себе жесткая? То есть, что многогранник задаёт локальный минимум задаче. А то я довольно долго (ошибочно) полагал, что многогранник для 10 точек на рисунке слева ниже является оптимальным решением, и только недавно понял, что он не просто не является оптимальным, он даже локально оптимальным не является. Если разорвать жёсткое ребро 1-2 и, рассматривая угол 1-0-2 как параметр, начать его увеличивать, то весь многогранник поплывёт. При этом длина жёстких рёбер будет расти до тех пор, пока вершины 1 и 2 не упрутся жёсткими рёбрами в соответсвующие вершины под ними (на рисунке слева они соединены парой тонких чёрных линий). В результате получится фигура на рисунке справа. У неё уже 19 жёстких рёбер (на одно больше, чем у предыдущей) и она доставляет задаче глобальный максимум.

Можно ли как-то это (локальную оптимальность) проверить с помощью систематического расчёта? Может быть через матрицу производных как-нибудь?


-------------------------
Код для расчёта картинок причесал, теперь не стыдно показать. Если кому вдруг тоже будет интересно поиграться.

Файл расчёта координат многогранника для 10 точек (оптимальное решение) calc_10.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact

points_func = @ (w) [
    w(1)        0
    w(1)        pi
    w(2)        pi / 2 - w(3)
    w(2)        pi / 2 + w(3)
    pi - w(4)   0
    pi - w(5)   pi / 2
];
dist_select_func = @ (d) sum (([d(1, 3), d(3, 4), d(3, 5), d(3, 6), d(5, 6)] - d(1, 2)) .^ 2);
target_func = @ (w) dist_select_func (calc_dist_func (calc_sphere_to_cartesian (points_func (w))));

w0 = [0.6 1.4 0.6 1.1 0.4];

fmso = optimset ('fminsearch');
fmso = optimset (fmso, 'Display', 'final', 'TolFun', 1e-15, 'TolX', 1e-15, ...
    'MaxIter', 20000, 'MaxFunEvals', 4000);
w = fminsearch (@(w) target_func (w), w0, fmso);

tp = [
    w(5)        pi / 2
    w(5)       -pi / 2
    w(4)        0
    w(4)        pi
    pi - w(2)   pi / 2 - w(3)
    pi - w(2)   pi / 2 + w(3)
    pi - w(2)  -pi / 2 - w(3)
    pi - w(2)  -pi / 2 + w(3)
    pi - w(1)   0
    pi - w(1)   pi
];
num = size (tp, 1);
coords = calc_sphere_to_cartesian (tp);

calc_plot_func (coords)
view ([-24 15])
calc_save_image (num)
 

Файл расчёта координат многогранника для 10 точек (неверное решение) calc_10_wrong.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact

points_func = @ (w) [
    w(1)        0
    w(1)        2 * pi / 3
    w(2)        pi / 3
    pi - w(3)   0
    pi          0
];
dist_select_func = @ (d) sum (([d(1, 3), d(3, 4), d(4, 5)] - d(1, 2)) .^ 2);
target_func = @ (w) dist_select_func (calc_dist_func (calc_sphere_to_cartesian (points_func (w))));

w0 = [0.6 1.4 0.8];

fmso = optimset ('fminsearch');
fmso = optimset (fmso, 'Display', 'final', 'TolFun', 1e-15, 'TolX', 1e-15, ...
    'MaxIter', 20000, 'MaxFunEvals', 4000);
w = fminsearch (@(w) target_func (w), w0, fmso);

tp = [
    w(1)        0
    w(1)        2 * pi / 3
    w(1)       -2 * pi / 3
    w(2)        pi / 3
    w(2)       -pi / 3
    w(2)        pi
    pi - w(3)   0
    pi - w(3)   2 * pi / 3
    pi - w(3)  -2 * pi / 3
    pi          0
];
num = size (tp, 1);
coords = calc_sphere_to_cartesian (tp);

calc_plot_func (coords)
view ([-24 15])
calc_save_image ('points_10_wrong')
 

Файл отрисовки многогранника (наиболее трудоёмкая в плане программирования часть) calc_plot_func.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160089.html
%   B@R5uk, 2025, April

function calc_plot_func (coords, pcolor, hcolor, ecolor)

subplot (111)

%   Прозрачность граней
alpha_value = 0.6;
%   Цвет граней
if 1 == nargin
    pcolor = [0.8, 0.9, 1.0];
end
%   Цвет "жёстких" рёбер
if 2 >= nargin
    hcolor = 'b';
end
%   Цвет "мягких" рёбер
if 3 >= nargin
    ecolor = 'k';
end

%   Допустимая погрешность для компланарности вершин
tol_plane = 1e-10;
%   Допустимая погрешность для длины жёстких рёбер
tol_lngth = 1e-10;
%   Количество вершин многогранника
num = size (coords, 1);
%   Центр масс вершин для определения внутренней и внешней стороны каждой грани
center = mean (coords);
%   Матрица расстояний между каждой парой вершин (бесконечность на главной диагонали матрицы)
dists = calc_dist_func (coords);
%   Минимальное расстояние между точками для определения жестких рёбер
min_dist = min (dists (:));

plot3 (0, 0, 0)
hold on
%   Массив для хранения индексов вершин текущей грани
points_list = zeros (num, 1);
%   Список кодов троек вершин, проверку которых необходимо пропустить,
%   потому что они уже входят в грань с числом вершин четыре и более
skip_list = [];
skip_code_vector = [num ^ 2 num 1];
%   Массив флагов рёбер для отрисовки
edge_flags = zeros (num);
%   Построение многогранника как выпуклой оболочки его вершин
for k = 1 : num - 2
    %   Направление на центр
    cntr_vect = center - coords (k, :);
    for l = k + 1 : num - 1
        for m = l + 1 : num
            %   Пропустить тройку, если она в списке
            if 0 ~= sum (num * (num * k + l) + m == skip_list)
                continue
            end
            %   Нормаль к текущей грани
            norm_vect = cross_product (coords (l, :) - coords (k, :), coords (m, :) - coords (k, :))';
            tmp_value = cntr_vect * norm_vect;
            %   Количество вершин текущей грани
            points_count = 0;
            %   Флаг, помечающий, что текущая тройка вершин не принадлежит внешней грани
            flag = 0;
            %   Цикл ищет вершины в плоскости текущей тройки точек и
            %   проверяет, что все остальные лежат по ту же сторону от неё,
            %   что и центр масс вершин многогранника
            for n = 1 : num
                tmp_vect = coords (n, :) - coords (k, :);
                %   Для вершин в плоскости текущей тройки точек скалярное
                %   произведение нормали на направление на вершину должно
                %   равняться почти нулю
                if abs (tmp_vect * norm_vect) < tol_plane
                    points_count = points_count + 1;
                    points_list (points_count) = n;
                    continue
                end
                %   Скалярные произведения нормали на вектора направлений
                %   на вершину и на центр масс должны быть одного знака
                if 0 >= tmp_value * tmp_vect * norm_vect;
                    flag = 1;
                    break
                end
            end
            if flag
                continue
            end
            %   Координаты вершин грани
            face_coords = coords (points_list (1 : points_count), :);
            %   Для числа вершин грани четыре и более требуется
            %   нахождение выпуклой оболочки вершин в плоскости грани
            if 3 == points_count
                %   Грань является треугольником
                new_list = [1 2 3];
            else
                %   Затравка для построения -- угол наибольшей величины
                new_list = max_angle (face_coords);
                new_list (points_count) = 0;
                %   Обход всех вершин
                for n = 3 : points_count
                    new_list (n) = max_angle_next (face_coords, new_list (n - 2), new_list (n - 1));
                end
                %   Пополняем список троек вершин к пропуску
                skip_list = [skip_list, list_of_codes(points_list (1 : points_count), skip_code_vector, 3)]; %#ok<AGROW>
            end
            %   Отрисовка многоугольника с вершинами в правильном порядке
            fill3 (face_coords (new_list, 1), ...
                   face_coords (new_list, 2), ...
                   face_coords (new_list, 3), pcolor, 'EdgeAlpha', 0, 'FaceAlpha', alpha_value)
            %   Помечаем рёбра для отрисовки
            for n = 2 : points_count
                edge_flags (points_list (new_list (n - 1)), points_list (new_list (n))) = 1;
                edge_flags (points_list (new_list (n)), points_list (new_list (n - 1))) = 1;
            end
            edge_flags (points_list (new_list (1)), points_list (new_list (points_count))) = 1;
            edge_flags (points_list (new_list (points_count)), points_list (new_list (1))) = 1;
        end
    end
end
%   Количество жёстких рёбер
num_hard = 0;
%   Отрисовка рёбер
for k = 1 : num - 1
    for l = k + 1 : num
        %   Жесткие рёбра
        if dists (k, l) - min_dist < tol_lngth
            num_hard = num_hard + 1;
            plot3 (coords ([k l], 1), coords ([k l], 2), coords ([k l], 3), hcolor, 'LineWidth', 3)
            continue
        end
        %   Мягкие рёбра
        if edge_flags (k, l)
            plot3 (coords ([k l], 1), coords ([k l], 2), coords ([k l], 3), ecolor)
        end
    end
end
hold off
xlim ([-1, 1])
ylim ([-1, 1])
zlim ([-1, 1])
axis vis3d
axis off
set (gcf, 'Color', 'w')

disp (['Number of points:     ', num2str(num)])
disp (['Number of hard edges: ', num2str(num_hard)])
disp (['Minimal   distance:   ', num2str(min_dist, 15)])
disp (['Spherical distance:   ', num2str(2 * asind (min_dist / 2), 15), ' deg'])
disp (char (13))

end

function result = cross_product (a, b)

result = [a(2) * b(3) - a(3) * b(2), a(3) * b(1) - a(1) * b(3), a(1) * b(2) - a(2) * b(1)];

end

function result = max_angle (coords)

num = size (coords, 1);
min_value = 1;
result = [];
for m = 2 : num
    tmp_vect_1 = coords (m, :) - coords (1, :);
    tmp_vect_1 = tmp_vect_1' / sqrt (sum (tmp_vect_1 .^ 2));
    for n = 2 : num
        %   Максимальный угол соответсвует наименьшему скалярному
        %   произведению нормированных векторов
        tmp_vect_2 = coords (n, :) - coords (1, :);
        tmp_vect_2 = tmp_vect_2 / sqrt (sum (tmp_vect_2 .^ 2));
        value = tmp_vect_2 * tmp_vect_1;
        if min_value > value
            min_value = value;
            result = [m 1 n];
        end
    end
end

end

function result = max_angle_next (coords, k, l)

num = size (coords, 1);
min_value = 1;
result = [];
tmp_vect_1 = coords (k, :) - coords (l, :);
tmp_vect_1 = tmp_vect_1' / sqrt (sum (tmp_vect_1 .^ 2));
for m = 1 : num
    tmp_vect_2 = coords (m, :) - coords (l, :);
    tmp_vect_2 = tmp_vect_2 / sqrt (sum (tmp_vect_2 .^ 2));
    value = tmp_vect_2 * tmp_vect_1;
    if min_value > value
        min_value = value;
        result = m;
    end
end

end

function result = list_of_codes (list, vector, num)

result = vector * nchoosek (list, num)';

end
 

Файл сохраняющий картинку в файл calc_save_image.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160089.html
%   B@R5uk, 2025, April

function calc_save_image (indx)

color_map = [
    255 255 255
    0   0   0
    0   0   255
    128 144 160
    208 240 255
    128 144 255
] / 255;
pause (0.5)
frmdata = getframe (gcf ());
imgdata = rgb2ind (frmdata .cdata, color_map, 'nodither');
[height, width] = size (imgdata);
imgdata = imcrop (imgdata, [floor(width / 2) - 128, floor(height / 2) - 128, 256, 256]);
if isnumeric (indx)
    filename = ['points_', num2str(indx), '.gif'];
else
    filename = [indx, '.gif'];
end
imwrite (imgdata, color_map, filename, 'GIF', 'TransparentColor', 0)

end
 

Файл вспомогательной функции (этот и следующий) calc_sphere_to_cartesian .m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160089.html
%   B@R5uk, 2025, April

function result = calc_sphere_to_cartesian (tp)

result = [sin(tp (:, 1)) .* cos(tp (:, 2)) sin(tp (:, 1)) .* sin(tp (:, 2)) cos(tp (:, 1))];

end
 

Файл calc_dist_func .m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160089.html
%   B@R5uk, 2025, April

function result = calc_dist_func (coords)

num = size (coords, 1);
result = repmat (permute (coords, [3, 1, 2]), num, 1);
result = sqrt (sum ((permute (result, [2, 1, 3]) - result) .^ 2, 3));
result = result + diag (inf (num, 1));

end
 

[мат-ламер]

B@R5uk в сообщении #1680963 писал(а):

Я так понимаю, в задачнике делался упор на доказательство оптимальности? Потому что там только до шести точек включительно, а самое интересное начинается от семи и более.

Да, на доказательство. Задачник для школьников. Авторы решили, что для $n>6$ что-то доказать для школьников будет сложновато.

B@R5uk в сообщении #1680963 писал(а):

Я тоже эту задачу пытаюсь численно решать.

У вас в планах доказать что-то или просто пока поэкспериментировать численно? Думаю, что что-то доказать будет достаточно сложно. Хотя строго доказанные конструкции есть для $n<11$ и $n=12$ , $n=14$ и $n=24$ . Может ещё для каких $n$ , но я не в курсе.

B@R5uk в сообщении #1680963 писал(а):

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

Не понял трудности. Ребро, соединяющее две вершины, принадлежит выпуклой оболочке этих двух вершин. Но может я не так понял проблему.

B@R5uk в сообщении #1681006 писал(а):

Тут у меня возник другой вопрос: есть ли какой-нибудь алгоритмизируемый способ убедиться, что заданная конструкция жёстких рёбер сама по себе жесткая? То есть, что многогранник задаёт локальный минимум задаче.

Тут я тоже не понял. Проверить, является ли многогранник жёстким, наверное можно с помощью подпрограмм MatLaba. Но вроде вам нужно что-то другое. Мы ведь находимся в рамках исходной задачи. И у нас точки находятся на сфере. Наверное вам нужно:

B@R5uk в сообщении #1681006 писал(а):

То есть, что многогранник задаёт локальный минимум задаче

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

[B@R5uk]

мат-ламер в сообщении #1681050 писал(а):

У вас в планах доказать что-то или просто пока поэкспериментировать численно?

Численно. Хочется построить алгоритм, который бы по заданному N быстро находил какой-нибудь локальный минимум. Пока мой отжиг зависит от кучи эвристических параметров, которые надо подкручивать в зависимости от N. Но даже с этим, для определённых значений N он не сходится точно в локальный экстремум или делает это очень плохо. Вот тут по ссылке есть некоторые известные результаты, которыми я себя проверяю. Точность там, правда, не очень высокая.

мат-ламер в сообщении #1681050 писал(а):

Не понял трудности. Ребро, соединяющее две вершины, принадлежит выпуклой оболочке этих двух вершин. Но может я не так понял проблему.

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

мат-ламер в сообщении #1681050 писал(а):

Тут я тоже не понял.

Я имею в виду следующее. Если пошевелить точки на сфере и при этом найдётся такое их движение, что минимальное ребро увеличится, то ранее найденное решение не является локальным максимумом. Хотелось бы легко угадывать "направление" в котором надо шевелить точки, чтобы построить быстро сходящийся алгоритм, или убедиться, что расположение точек локально оптимально.

мат-ламер в сообщении #1681050 писал(а):

Но основная серьёзная трудность будет - как убедиться. что найденный локальный максимум является глобальным?

Мне бы пока с локальным экстремумом разобраться.

мат-ламер в сообщении #1681050 писал(а):

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

Да, я как раз имел в виду что-нибудь в духе матрицы Гёссе для метода множителей Лагранжа. Однако проблема Таммеса в точке экстремума существенно нелинейна. Задачу можно формализовать следующим способом. Пусть вектор P — это список параметров, задающих на единичной сфере координаты точек (число которых N). Тогда задача является поиском экстремума целевой функции f:

$$\mathbf{\widehat{P}}=\max\limit_{\mathbf{P}}f\bigl(\mathbf{P}\bigr)$$

где целевая функция выражается как

$$f\bigl(\mathbf{P}\bigr)=\min\limit_{1\le m<n\le N}\Bigl\lVert\mathbf{r_m}\bigl(\mathbf{P}\bigr)-\mathbf{r_n}\bigl(\mathbf{P}\bigr)\Bigr\rVert$$

или вот так (чтобы сэкономить на вычислении квадратных корней):

$$f\bigl(\mathbf{P}\bigr)=\min\limit_{1\le m<n\le N}\Bigl\lVert\mathbf{r_m}\bigl(\mathbf{P}\bigr)-\mathbf{r_n}\bigl(\mathbf{P}\bigr)\Bigr\rVert^2$$

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

$$\mathbf{P}=\left(\begin{array}{l}p_1\\p_2\\\ldots\\p_{2N-1}\\p_{2N}\end{array}\right),\qquad\mathbf{r_n}\bigl(\mathbf{P}\bigr)=\left(\begin{array}{l}\mathbf{e_x}\\\mathbf{e_y}\\\mathbf{e_z}\end{array}\right)^T\left(\begin{array}{l}\sin p_{2n-1}\cos p_{2n}\\\sin p_{2n-1}\sin p_{2n}\\\cos p_{2n-1}\end{array}\right),\qquad 1\le n\le N$$

В точке экстремума целевая функция имеет значение:

$$d=f\bigl(\mathbf{\widehat{P}}\bigr)$$

это означает, что некоторый набор рёбер имеет в качестве длины как раз это значение. Я называю их жёсткими:

$$d=\Bigl\lVert\mathbf{r_m}\bigl(\mathbf{\widehat{P}}\bigr)-\mathbf{r_n}\bigl(\mathbf{\widehat{P}}\bigr)\Bigr\rVert,\qquad\bigl\{\;m,\;n\;|\;1\le m<n\le N\;\bigr\}\in S$$

Здесь S — это множество пар индексов, указывающих эти жёсткие рёбра.

Проблема здесь в том, что целевая функция f представляет собой лоскутное одеяло, скроенное из кусков цилиндров (или параболоидов, если в определении стоит квадрат). Разумеется, сшивка лоскутков не является гладкой, там у неё изломы. А в точке экстремума (не важно, локального или глобального) у неё вообще вершина, в которой сходится куча рёбер и лоскутков. Причём не имеет значения, стоит ли в формуле для целевой функции квадрат или нет: квадрат превращает остриё каждого лоскута в чашку на достаточном удалении от искомых локальных максимумов целевой функции (удалены они по крайней мере на величину d). Угадать направление роста целевой функции тоже не представляется возможным. В многомерном пространстве вектора параметров P искомое направление может быть очень малым подмножеством всех возможных направлений.

Наглядно поведение целевой функции можно рассмотреть на такой модельной функции двух переменных:

$$f(x,y)=|x|+|y|$$

Её график в трёхмерном пространстве представляет собой сшитые куски четырёх плоскостей, что выглядит как перевёрнутая пирамида с основанием на бесконечности и рёбрами в точности над осями координат OXY. В начале координат у неё экстремум, но производной в этой точке у функции нет. Так же, находясь на осях координат на удалении от начала координат значение функции можно улучшать, двигаясь к началу. Но такое положение находится в "желобе" графика, на ребре сшивки, где производной тоже нет. Поэтому применить градиентный спуск (что я бы очень хотел сделать) здесь в лоб не получится.

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