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

[мат-ламер]

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

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

Думаю попробовать такой подход. Вводим новую переменную $\rho$ , которую будем максимизировать $\rho \to \max$ . При этом возникают ограничения в форме неравенств $d^2(R_i,R_j) \ge \rho^2$ , где $d^2$ - квадрат расстояния между точками $R_i$ и $R_j$ . Также должны быть выполнены ограничения в форме равенств $\left\lVert R_i \right\rVert^2 =1$ .

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

(Оффтоп)

Чуток подправил обозначения.

[drzewo]

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

Сразу скажу, это НЕ задача Томсона.

ok, а для этой задачи вариационный принцип есть? Функционал критической точкой которого была бы требуемая конфигурация?

[B@R5uk]

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

$$\begin{array}{ccccll} \text{Кол-во}&\text{Кол-во}&\text{Величина}&\text{Излишек}&\text{Длина жёстких}&\text{Сферическое}\\ \text{вершин }N&\text{связей}&2N-2&\text{связей}&\text{рёбер (связей)}&\text{расстояние, }^\circ\\ \hline 6&12&10&2&1.41421356237309&90\\ 7&12&12&0&1.25687046848066&77.8695421554234\\ 8&16&14&2&1.21556252413132&74.8584921856155\\ 9&18&16&2&1.15470053837925&70.5287793655093\\ 10&19&18&1&1.09142629140382&66.1468219878905\\ 11&25&20&5&1.05146222423827&63.4349488229220\\ 12&30&22&8&1.05146222423827&63.4349488229220\\ 13&24&24&0&0.956413627224499&57.1367030775783\\ 14&28&26&2&0.933862623344028&55.6705699956595\\ 15&30&28&2&0.902656188229966&53.6578501299326\\ 16&32&30&2&0.880574112176032&52.2443957527206\\ 17&34&32&2&0.862444879316907&51.0903285516277\\ \hline \end{array}$$

Визуализация этих многогранников ("толстые" рёбра синего цвета являются жёсткими) и точно рассчитанные (на сколько разрядность позволяет) координаты в конце под катом:

$\begin{array}{ccccc}6&\qquad\qquad\qquad\qquad\qquad\qquad&7&\qquad\qquad\qquad\qquad\qquad\qquad&8\\9&&10&&11\end{array}$

$\begin{array}{ccccc}12&\qquad\qquad\qquad\qquad\qquad\qquad&13&\qquad\qquad\qquad\qquad\qquad\qquad&14\\15&&16&&17\end{array}$

Весьма замечательным является 15-вершинник (первая картинка в последней строке). Этот многогранник не является симметричным (и поэтому имеет хиральность), не смотря на то, что он имеет целых 11 граней, являющихся равносторонними треугольниками. Эти треугольники группируются в два блока по 5 и 6 штук в каждом (на рисунке снизу и сверху, соответственно). При этом, если удалить три точки так, чтобы сверху и снизу осталось по 4 треугольника (в конфигурации, образующей один большой треугольник с двойной стороной), то получающаяся фигура внезапно становится симметричной. Она будет иметь одну вращательную симметрию третьего порядка (вдоль оси OZ) и три вращательных симметрий второго порядка, оси которых лежат в плоскости OXY. Центральные треугольники этих двух блоков параллельны друг другу и этой самой горизонтальной плоскости OXY (и повёрнуты друг к другу на некоторый угол). Фигуру с аналогичной симметрией так же можно получить, если добавить определённым образом (игнорируя минимальные расстояния) три точки так, чтобы сверху и снизу в блоках было по 7 равносторонних треугольников. Эта "скрытая" симметрия приводит к двум "лишним" жёстким ребрам: друг относительно друга блоки фиксируются 6-ю рёбрами, хотя для этого вполне было бы достаточно 4-х.

(Оффтоп)

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

=============================================================
Number of points:     6
Number of hard edges: 12
Minimal   distance:   1.41421356237309
Spherical distance:   90 deg
Coordinates:
   0.816496580927726                   0   0.577350269189626
  -0.408248290463863   0.707106781186548   0.577350269189626
  -0.408248290463863  -0.707106781186548   0.577350269189626
   0.408248290463863   0.707106781186547  -0.577350269189626
   0.408248290463863  -0.707106781186547  -0.577350269189626
  -0.816496580927726   0.000000000000000  -0.577350269189626

=============================================================
Number of points:     7
Number of hard edges: 12
Minimal   distance:   1.25687046848066
Spherical distance:   77.8695421554234 deg
Coordinates:
   0.725654503313800                   0   0.688059257491971
  -0.362827251656900   0.628435234240330   0.688059257491971
  -0.362827251656900  -0.628435234240330   0.688059257491971
   0.488835833773143   0.846688500655378  -0.210138312730603
   0.488835833773143  -0.846688500655378  -0.210138312730603
  -0.977671667546286   0.000000000000000  -0.210138312730603
   0.000000000000000                   0  -1.000000000000000

=============================================================
Number of points:     8
Number of hard edges: 16
Minimal   distance:   1.21556252413132
Spherical distance:   74.8584921856155 deg
Coordinates:
   0.859532503769496                   0   0.511081084529394
   0.000000000000000   0.859532503769496   0.511081084529394
   0.000000000000000  -0.859532503769496   0.511081084529394
  -0.859532503769496   0.000000000000000   0.511081084529394
   0.607781262065662   0.607781262065662  -0.511081084529394
   0.607781262065662  -0.607781262065662  -0.511081084529394
  -0.607781262065662   0.607781262065662  -0.511081084529394
  -0.607781262065662  -0.607781262065662  -0.511081084529394

=============================================================
Number of points:     9
Number of hard edges: 18
Minimal   distance:   1.15470053837925
Spherical distance:   70.5287793655093 deg
Coordinates:
   0.666666666666667                   0   0.745355992499930
  -0.333333333333333   0.577350269189626   0.745355992499930
  -0.333333333333333  -0.577350269189626   0.745355992499930
   0.500000000000000   0.866025403784439   0.000000000000000
   0.500000000000000  -0.866025403784439   0.000000000000000
  -1.000000000000000   0.000000000000000   0.000000000000000
   0.666666666666667                   0  -0.745355992499930
  -0.333333333333333   0.577350269189626  -0.745355992499930
  -0.333333333333333  -0.577350269189626  -0.745355992499930

=============================================================
Number of points:     10
Number of hard edges: 19
Minimal   distance:   1.09142629140382
Spherical distance:   66.1468219878905 deg
Coordinates:
   0.000000000000000   0.609945710115265   0.792443203461286
   0.000000000000000  -0.609945710115265   0.792443203461286
   0.859988553050598                   0   0.510313323970615
  -0.859988553050598   0.000000000000000   0.510313323970615
   0.545713145701910   0.828261362520315  -0.127201721545415
  -0.545713145701910   0.828261362520315  -0.127201721545415
  -0.545713145701910  -0.828261362520315  -0.127201721545415
   0.545713145701910  -0.828261362520315  -0.127201721545415
   0.545713145701911                   0  -0.837972053596136
  -0.545713145701911   0.000000000000000  -0.837972053596136

=============================================================
Number of points:     11
Number of hard edges: 25
Minimal   distance:   1.05146222423827
Spherical distance:   63.434948822922 deg
Coordinates:
   0.894427190999916                   0   0.447213595499958
   0.276393202250021   0.850650808352040   0.447213595499958
   0.276393202250021  -0.850650808352040   0.447213595499958
  -0.723606797749979   0.525731112119134   0.447213595499958
  -0.723606797749979  -0.525731112119134   0.447213595499958
   0.723606797749979   0.525731112119134  -0.447213595499958
   0.723606797749979  -0.525731112119134  -0.447213595499958
  -0.276393202250021   0.850650808352040  -0.447213595499958
  -0.276393202250021  -0.850650808352040  -0.447213595499958
  -0.894427190999916   0.000000000000000  -0.447213595499958
   0.000000000000000                   0  -1.000000000000000

=============================================================
Number of points:     12
Number of hard edges: 30
Minimal   distance:   1.05146222423827
Spherical distance:   63.434948822922 deg
Coordinates:
                   0                   0   1.000000000000000
   0.894427190999916                   0   0.447213595499958
   0.276393202250021   0.850650808352040   0.447213595499958
   0.276393202250021  -0.850650808352040   0.447213595499958
  -0.723606797749979   0.525731112119134   0.447213595499958
  -0.723606797749979  -0.525731112119134   0.447213595499958
   0.723606797749979   0.525731112119134  -0.447213595499958
   0.723606797749979  -0.525731112119134  -0.447213595499958
  -0.276393202250021   0.850650808352040  -0.447213595499958
  -0.276393202250021  -0.850650808352040  -0.447213595499958
  -0.894427190999916   0.000000000000000  -0.447213595499958
   0.000000000000000                   0  -1.000000000000000

=============================================================
Number of points:     13
Number of hard edges: 24
Minimal   distance:   0.956413627224499
Spherical distance:   57.1367030775783 deg
Coordinates:
   0.676286561429666                   0   0.736638640603138
   0.000000000000000   0.676286561429666   0.736638640603138
   0.000000000000000  -0.676286561429666   0.736638640603138
  -0.676286561429666   0.000000000000000   0.736638640603138
   0.704230456710091   0.704230456710091   0.090105092440959
   0.704230456710091  -0.704230456710091   0.090105092440959
  -0.704230456710091   0.704230456710091   0.090105092440959
  -0.704230456710091  -0.704230456710091   0.090105092440959
   0.839967644115645                   0  -0.542636486829639
   0.000000000000000   0.839967644115645  -0.542636486829639
   0.000000000000000  -0.839967644115645  -0.542636486829639
  -0.839967644115645   0.000000000000000  -0.542636486829639
   0.000000000000000                   0  -1.000000000000000

=============================================================
Number of points:     14
Number of hard edges: 28
Minimal   distance:   0.933862623344028
Spherical distance:   55.6705699956595 deg
Coordinates:
                   0                   0   1.000000000000000
   0.681127894674377   0.466931311672014   0.563950300360505
   0.681127894674377  -0.466931311672014   0.563950300360505
  -0.681127894674377  -0.466931311672014   0.563950300360505
  -0.681127894674376   0.466931311672014   0.563950300360505
   0.000000000000000   0.982441144573429   0.186572767169410
   0.000000000000000  -0.982441144573429   0.186572767169410
   0.982441144573429                   0  -0.186572767169410
  -0.982441144573429   0.000000000000000  -0.186572767169410
  -0.466931311672014   0.681127894674376  -0.563950300360505
   0.466931311672014   0.681127894674377  -0.563950300360505
   0.466931311672014  -0.681127894674377  -0.563950300360505
  -0.466931311672014  -0.681127894674376  -0.563950300360505
   0.000000000000000                   0  -1.000000000000000

=============================================================
Number of points:     15
Number of hard edges: 30
Minimal   distance:   0.902656188229966
Spherical distance:   53.6578501299326 deg
Coordinates:
   0.260574396630126   0.451328094114983   0.853465837209306
   0.260574396630126  -0.451328094114983   0.853465837209306
  -0.521148793260253   0.000000000000000   0.853465837209306
   0.908985923112752                   0   0.416826812456754
  -0.454492961556376   0.787204901098092   0.416826812456754
  -0.454492961556376  -0.787204901098092   0.416826812456754
   0.376835267797373   0.921711972287309   0.091881560099513
   0.609808349074382  -0.787204901098091   0.091881560099513
  -0.965392708383380  -0.244079285302712  -0.091881560099513
   0.843082907030753   0.339803796755889  -0.416826812456754
  -0.715820173808380   0.560229316607122  -0.416826812456754
  -0.127262733222372  -0.900033113363010  -0.416826812456754
   0.072963527992517   0.516015879890826  -0.853465837209306
   0.410401096745376  -0.321196208736669  -0.853465837209306
  -0.483364624737893  -0.194819671154156  -0.853465837209306

=============================================================
Number of points:     16
Number of hard edges: 32
Minimal   distance:   0.880574112176032
Spherical distance:   52.2443957527206 deg
Coordinates:
   0.622659926056996                   0   0.782492566407309
   0.000000000000000   0.622659926056996   0.782492566407309
   0.000000000000000  -0.622659926056996   0.782492566407309
  -0.622659926056996   0.000000000000000   0.782492566407309
   0.687190046418906   0.687190046418906   0.235668581286442
   0.687190046418906  -0.687190046418906   0.235668581286442
  -0.687190046418906   0.687190046418906   0.235668581286442
  -0.687190046418906  -0.687190046418906   0.235668581286442
   0.971833483573413                   0  -0.235668581286442
   0.000000000000000   0.971833483573413  -0.235668581286442
   0.000000000000000  -0.971833483573413  -0.235668581286442
  -0.971833483573413   0.000000000000000  -0.235668581286442
   0.440287056088016   0.440287056088016  -0.782492566407309
   0.440287056088016  -0.440287056088016  -0.782492566407309
  -0.440287056088016   0.440287056088016  -0.782492566407309
  -0.440287056088016  -0.440287056088016  -0.782492566407309

=============================================================
Number of points:     17
Number of hard edges: 34
Minimal   distance:   0.862444879316907
Spherical distance:   51.0903285516277 deg
Coordinates:
   0.000000000000000   0.553910528530291   0.832576198544790
   0.000000000000000  -0.553910528530291   0.832576198544790
   0.656416466468028                   0   0.754398715898713
  -0.656416466468028   0.000000000000000   0.754398715898713
   0.638757213967231   0.717900538765430   0.276781570999816
   0.638757213967231  -0.717900538765430   0.276781570999816
  -0.638757213967231  -0.717900538765430   0.276781570999816
  -0.638757213967231   0.717900538765430   0.276781570999816
   0.999317174356973                   0  -0.036948410455591
  -0.999317174356973   0.000000000000000  -0.036948410455591
   0.000000000000000   0.969459500917063  -0.245251454800248
   0.000000000000000  -0.969459500917063  -0.245251454800248
   0.605300639615620   0.488987261018904  -0.628094415070021
   0.605300639615620  -0.488987261018904  -0.628094415070021
  -0.605300639615620  -0.488987261018904  -0.628094415070021
  -0.605300639615620   0.488987261018904  -0.628094415070021
   0.000000000000000                   0  -1.000000000000000
 

[мат-ламер]

drzewo в сообщении #1681277 писал(а):

ok, а для этой задачи вариационный принцип есть? Функционал критической точкой которого была бы требуемая конфигурация?

Он тут обсуждается последние две страницы. Исходный функционал $\min d(R_i,R_j) \to \max$ , $\left\lVert R_i\right\rVert=1$ (максимизируется минимальное расстояние между точками $R_i$ на сфере) не очень приятный - он негладкий, не выпуклый и не вогнутый, многоэкстремальный. Условие экстремума (локального) пока для него в теме не выложено. Но как-нибудь надо будет этим заняться. В прошлом своём посту задачу чуток переформулировал (естественно, с сохранением решения), чтобы функционал стал гладкий. Но появилась куча новых ограничений, причём типа неравенств. Для него можно выписать своё условие экстремума.

(Оффтоп)

Чуток подправил обозначения.

[B@R5uk]

мат-ламер, лучше бы вы точки в своих обозначениях записывали как радиус-вектора $\mathbf{R}_i$, потому что буква P у меня зарезервирована под столбец всех координат/параметров системы, тем более, что вы уже трактуете своё обозначение как вектор, записывая $\bigl\lVert\mathbf{R}_i\bigr\rVert=1$. А то будем путаться и других путать, обозначая разные вещи одной буквой.

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

Для начала, ещё раз постановка (немного модифицированная). Пусть имеется N точек на единичной сфере:

$$\Bigl\lVert\mathbf{R}\bigl(P,\;n\bigr)\Bigr\rVert^2=x_n^2+y_n^2+z_n^2=1$$

где P — это столбец координат системы:

$$P=\bigl(x_1,\;y_1,\;z_1,\;x_2,\;y_2,\;\ldots,\;z_{N-1},\;x_N,\;y_N,\;z_N\bigr)^T$$

а функция R достаёт из этого столбца три координаты n-ой точки и превращает их в вектор:

$$\mathbf{R}\bigl(P,\;n\bigr)=x_n\mathbf{e}_x+y_n\mathbf{e}_y+z_n\mathbf{e}_z$$

Задача состоит в том, чтобы найти максимум целевой функции:

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

которая равна

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

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

Обозначая величину целевой функции в точке P пространства параметров греческой буквой ро, а расстояние между точками с индексами m и n прописной буквой l имеем одно равенство и серию неравенств:

$$\rho=f\bigl(P\bigr)=\min\limit_{1\le m<n\le N}l\bigl(m,\;n\bigr)$$

$$\rho\le l\bigl(m,\;n\bigr)=\Bigl\lVert\mathbf{R}\bigl(P,\;m\bigr)-\mathbf{R}\bigl(P,\;n\bigr)\Bigr\rVert^2$$

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

$$\rho=l\bigl(m_k,\;n_k\bigr)=l_k,\qquad 1\le m_k<n_k\le N,\qquad 1\le k\le K$$

Выше все переменные l с индексом равны по величине, но это разные переменные, потому что они по-разному зависят от столбца параметров P. К чему это приведёт будет понятно ниже. Составим из них столбец и обозначим его L:

$$L=\bigl(l_1,\;l_2,\;\ldots,\;l_{K-1},\;l_K\bigr)^T$$

Теперь меня интересуют такие вопросы:
1) Для заданного расположения точек P фиксируют ли жёсткие рёбра, заданные парами индексов вершин $\bigl(m_k,\;n_k\bigr)$ форму многогранника или нет?
2) Какие из заданных жёстких рёбер являются следствием наличия других жёстких рёбер и прочих условий задачи?
3) Непрерывно двигаясь из исходной точки P пространства параметров системы существует ли возможность улучшить значение целевой функции?

Первый и последний вопросы явно подразумевают некого рода движение точек или, что то же самое, движение точки P в пространстве параметров. Логично это движение обозначать дифференциалом этого столбца параметров:

$$dP=\bigl(dx_1,\;dy_1,\;dz_1,\;dx_2,\;dy_2,\;\ldots,\;dz_{N-1},\;dx_N,\;dy_N,\;dz_N\bigr)^T$$

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

$$\frac{1}{2}d\left(\Bigl\lVert\mathbf{R}\bigl(P,\;n\bigr)\Bigr\rVert^2\right)=\frac{1}{2}d\Bigl(x_n^2+y_n^2+z_n^2\Bigr)=x_ndx_n+y_ndy_n+z_ndz_n=0,\qquad 1\le n\le N$$

Получается система из N уравнений, связывающая 3N координат системы с их дифференциалами. Как и исходную систему ограничений, этот набор равенств можно записывать по-разному, в том числе через свёртку с тензором третьего ранга. Поскольку меня сейчас в первую очередь интересуют дифференциалы, а параметры на данном этапе являются константами, я запишу эту систему в матричной форме, где матрица Q является функцией точки P пространства параметров:

$$Q_S\bigl(P\bigr)dP=0_N$$

где индекс S у матрицы означает, что она получается как следствие факта принадлежности точек сфере, а индекс N у нуля — что это на самом деле столбец из N нулей. Сама матрица Q имеет N строк и 3N столбцов и представляет собой кучу нулей с координатами точек на скошенной диагонали одна тройка за другой с верхнего левого угла в правый нижний по три величины в строке по одной величине в столбце:

$$Q_S\bigl(P\bigr)=\left(\begin{matrix} x_1&y_1&z_1&0&\ldots&0&0&0&0\\ 0&0&0&x_2&\ldots&0&0&0&0\\ \ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots\\ 0&0&0&0&\ldots&z_{N-1}&0&0&0\\ 0&0&0&0&\ldots&0&x_N&y_N&z_N\\ \end{matrix}\right)$$

Во-вторых, движение точек ограничено наличием жёстких рёбер. Рассмотрим дифференциал dL столбца L одинаковых величин, введённого выше. Его компоненты:

$$dl_k=d\Bigl(l\bigl(m_k,\;n_k\bigr)\Bigr)=d\left(\Bigl\lVert\mathbf{R}\bigl(P,\;m_k\bigr)-\mathbf{R}\bigl(P,\;n_k\bigr)\Bigr\rVert^2\right)=$$

$$=d\left(\bigl(x_{m_k}-x_{n_k}\bigr)^2+\bigl(y_{m_k}-y_{n_k}\bigr)^2+\bigl(z_{m_k}-z_{n_k}\bigr)^2\right)=$$

$$=2\bigl(x_{m_k}-x_{n_k}\bigr)\bigl(dx_{m_k}-dx_{n_k}\bigr)+2\bigl(y_{m_k}-y_{n_k}\bigr)\bigl(dy_{m_k}-dy_{n_k}\bigr)+2\bigl(z_{m_k}-z_{n_k}\bigr)\bigl(dz_{m_k}-dz_{n_k}\bigr)$$

Опять же, этот набор равенств можно записать через тензор третьего ранга или, как я делаю, в матричной форме:

$$Q_E\bigl(P,\;H\bigr)dP=dL$$

Здесь индекс E означает, что матрица Q являются следствием соотношений на жёсткие рёбра (Edges), а H — это множество K пар чисел, являющимися индексами вершин, соединёнными этими жёсткими рёбрами:

$$H=\Bigl\langle\;\left(m_k,\;n_k\right)\;\Bigl|\Bigr{}\;1\le k\le K\Bigr\rangle$$

Матрица Q содержит K строк и 3N столбцов, и, как и предыдущая, в большинстве своём наполнена нулями. В каждой строке имеется по 6 чисел, являющихся разностью пар координат точек в соответствии с равенствами выше. Подробно расписывать не буду, потому что это довольно муторно.

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

$$dL=0_K$$

Вместе со сферическим ограничением это приведёт к такой системе уравнений:

$$\left(\begin{matrix}Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)dP=0_{N+K}$$

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

Одну другую вещь можно сказать наверняка. Даже если система переопределена:

$$N+K>3N$$

ранг её матрицы не может быть быть больше чем

$$\operatorname{rank}\left(\begin{matrix}Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)\le 3N-3$$

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

$$d\mathbf{M}=\sum\limits_{n=1}^N\Bigl[\mathbf{R}\bigl(P,\;n\bigr),\;d\mathbf{R}\bigl(dP,\;n\bigr)\Bigr]$$

Приравнять нулю и записать это выражение в матричной форме:

$$Q_R\bigl(P\bigr)dP=0_3$$

Здесь индекс R у матрицы Q означает, что она призвана убрать вращение (Rotation) из системы. Она имеет три строки и 3N столбца, все элементы заполнены координатами точек с плюсами и минусами в соответствии с суммой векторных произведений. Надо заметить, что я не знаю, как обобщить этот шаг на размерности пространства больше 3, потому что я не знаю, каков аналог будет там у этой векторной операции, и что вообще представляют из себя момент импульса и инерции в пространствах высшей размерности. Все остальные этапы такое обобщение допускают без проблем.

Теперь, если дополненное уравнение

$$\left(\begin{matrix}Q_R\bigl(P\bigr)\\Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)dP=0_{N+K+3}$$

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

$$\operatorname{rank}\left(\begin{matrix}Q_R\bigl(P\bigr)\\Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)<3N$$

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

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

$$\left(\begin{matrix}Q_R\bigl(P\bigr)\\Q_S\bigl(P\bigr)\\Q_E\bigl(P,\;H\bigr)\end{matrix}\right)dP=\left(\begin{matrix}0_{N+3}\\dL\end{matrix}\right)$$

при условии, что все компоненты дифференциала dL неотрицательны:

$$dl_k\le 0,\qquad 1\le k\le K$$

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

Так же меня смущает мой поход к "удалению" вращательных степеней свободы. Фактически, при записи соответствующих уравнений я назначил всем точкам единичные массы. Что будет если сделать массы не единичными и отличающимися друг от друга? Будут ли три новые строчки матрицы системы линейно зависимы с ранее полученными тремя? Если нет, то что-то мутное творится с этим подходом: можно добавлять строчки до бесконечности и в конце-концов переопределить систему на столько, что из неё нельзя будет чего-то путного извлечь. И мой поход будет неверен. Если же три новых уравнения (с произвольными массами) будут линейно зависимы с уравнениями с единичными массами, то это ещё надо доказать (для доказательства корректности подхода), потому что это совсем не очевидно.

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

Есть ещё проблема определения ранга матрицы в случае, когда на промежуточных этапах поиска экстремума точность значений длин жёстких рёбер равна лишь с некоторой погрешностью. То есть они даже не являются строго говоря жёсткими, только с определённой натяжкой. Можно ли ранг матрицы вычислять с произвольной наперёд заданной точностью? Что считать критерием отсечки? Как действовать в такой ситуации мне пока непонятно.

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

(Оффтоп)

А ещё я что-то раскатался в кнопкожмяканье и перевалил за половину величины потолка на длину сообщения в 20 тысяч символов. Ну хотя бы для себя прояснил решения, проблемы и вопросы, даже если никто до конца не дочитал.

[zgemm]

Я сейчас скажу одну глупость, но, возможно Вам поможет. У меня лет 20 назад была такая же задача. Для случаев меньше икосаэдра мне не было нужно. Все, что больше делал очень инженерно и бездоказательно, но, на удивление, всегда работало.

Исходим из того, что каждая грань икосаэдра заполняется равносторонними треугольниками, а потом полученную фигуру можно раздуть до того, что она ляжет на сферу. Назовем это дополненно-раздудый икосаэдр. Для заданного $N$ ищем ближайший сверху такой дополненно-раздутый, и ищем на сколько он отличается от того, что надо нам. Далее находим две параллели близкие по значению одну в северном, другую в южном полушариях, на которых в сумме будет столько точек, как если отнять числа точек в этом раздуто-дополненного икосаэдре значение $N$. Удаляем точки на этих двух параллелях. Разрешаем каждой точке двигаться только вдоль ее меридиана, и пододвигаем так, чтобы дырки от удаленных точек схлопнулись. Запускаем минимизацию, уже без отжига обычным BFGSом. Было предпринято много раз улучшить результат отжигом, но не улучшалось.

PS: не ругаться на меня, что написал не строго и очень инженерно. Надеюсь ход мыслей донес.

PPS: ой, вспомнил, что в моем случае разрешено было сделать так, что число точек в фигуре могло быть на несколько единиц меньше, чем $N$.

[Евгений Машеров]

Dmitriy40 в сообщении #1680930 писал(а):

Разве это не Задача Томсона
?

В задаче Томсона - минимизация суммы обратных величин расстояний. В данной задаче - максимизация минимального расстояния.

[B@R5uk]

Покрутил тут ещё немного 15-вершинник и обнаружил следующую замечательную вещь. Я выше писал, что у ранее изображённой фигуры имеется три точки, которые надо удалить, чтобы эта фигура стала симметричной, а так же три точки, которые надо добавить, игнорируя минимальные расстояния, чтобы она так же стала симметричной. Так вот, эти шесть точек разбиваются в три пары взаимозаменяемых точек. Любую точку из пары можно взять, и в результате получится оптимальная фигура с одним и тем же минимаксным расстоянием между точками. Причём положение всех остальных 12 вершин не меняется. Это, в частности, означает, что можно выбрать точки так (7 равносторонних треугольников сверху и 4 снизу, или наоборот), чтобы результирующий многогранник был симметричным с осью симметрии OZ третьего порядка. На картинке ниже этот многогранник имеет синие рёбра. Красными рёбрами помечены альтернативные положения вершин:

Результаты расчёта координат:

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

=============================================================
Number of points:     15
Number of hard edges: 30
Minimal   distance:   0.902656188229966
Spherical distance:   53.6578501299326 deg
Coordinates:
   0.260574396630126   0.451328094114983   0.853465837209306
   0.260574396630126  -0.451328094114983   0.853465837209306
  -0.521148793260253   0.000000000000000   0.853465837209306
   0.072963527992517   0.516015879890826  -0.853465837209306
   0.410401096745376  -0.321196208736669  -0.853465837209306
  -0.483364624737893  -0.194819671154156  -0.853465837209306
   0.908985923112752                   0   0.416826812456754
  -0.454492961556376   0.787204901098092   0.416826812456754
  -0.454492961556376  -0.787204901098092   0.416826812456754
   0.843082907030753   0.339803796755889  -0.416826812456754
  -0.715820173808380   0.560229316607122  -0.416826812456754
  -0.127262733222372  -0.900033113363010  -0.416826812456754
   0.376835267797373   0.921711972287309   0.091881560099513  or   0.271317492581992   0.958094252739625  -0.091881560099513
   0.609808349074382  -0.787204901098091   0.091881560099513  or   0.694075215801388  -0.714014967436913  -0.091881560099513
  -0.986643616871755  -0.134507071189218   0.091881560099513  or  -0.965392708383380  -0.244079285302712  -0.091881560099513
=============================================================
 

Файл calc_15.m с программой MATLAB для расчёта координат и отображения фигуры, если кто-то захочет покрутить этот многогранник и рассмотреть его со всех сторон. Причём, "покрутить" я имею в виду буквально: MATLAB позволяет вращать мышкой 3D-график во все стороны.

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = [0.55 1.1 1.45 0.1 0.3];

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

tp = [
    w(1)        pi / 3
    w(1)       -pi / 3
    w(1)        pi
    pi - w(1)   pi / 3     + w(5)
    pi - w(1)  -pi / 3     + w(5)
    pi - w(1)   pi         + w(5)
    w(2)        0
    w(2)        2 * pi / 3
    w(2)       -2 * pi / 3
    pi - w(2)   0          + w(5)
    pi - w(2)   2 * pi / 3 + w(5)
    pi - w(2)  -2 * pi / 3 + w(5)
    w(3)        pi / 3 + w(4)
    w(3)       -pi / 3 + w(4)
    w(3)        pi     + w(4)
    pi - w(3)   pi / 3 - w(4) + w(5)
    pi - w(3)  -pi / 3 - w(4) + w(5)
    pi - w(3)   pi     - w(4) + w(5)
];
num = size (tp, 1) - 3;
coords_all = calc_sphere_to_cartesian (tp);
coords = coords_all (1 : num, :);
indices_1 = [
    16      4
    16      8
    16      10
    17      5
    17      7
    17      12
    18      6
    18      9
    18      11
];
indices_2 = [
    16      7
    16      11
    17      9
    17      10
    18      8
    18      12
];

calc_plot_func (coords)
hold on
for k = 1 : 9
    plot3 (coords_all (indices_1 (k, :), 1), coords_all (indices_1 (k, :), 2), coords_all (indices_1 (k, :), 3), 'r', 'LineWidth', 3)
end
for k = 1 : 6
    plot3 (coords_all (indices_2 (k, :), 1), coords_all (indices_2 (k, :), 2), coords_all (indices_2 (k, :), 3), 'r')
end
hold off
view ([73 15])
color_map = [
    255 255 255
    0   0   0
    0   0   255
    128 144 160
    208 240 255
    128 144 255
    255 0   0
    208 160 192
] / 255;
calc_save_image (num, color_map)
 

Вспомогательный код для сохранения картинки в файл calc_save_image.m. Модифицирован (по сравнению с версией, приведённой ранее) для добавления возможности использования расширенной палитры цветов. Имеет обратную совместимость.

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

%   topic160089.html
%   B@R5uk, 2025, April

function calc_save_image (indx, color_map)

if 1 == nargin
    color_map = [
        255 255 255
        0   0   0
        0   0   255
        128 144 160
        208 240 255
        128 144 255
    ] / 255;
end
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
 

Под катом тексты программ для расчёта координат всех остальных многогранников, изображённых в посте выше. Кат удалил, потому что тег off конфликтует с тегом syntax, превращая круглые скобки в последовательности кодов символов типа "&#40;w&#41;". Не знаю, куда на этот глюк жаловаться.

Файл calc_6_a.m для шести точек:

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = 0.8;

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

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

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

Файл calc_7.m для семи точек:

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = [0.8 1.8];

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

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          0
];
num = size (tp, 1);
coords = calc_sphere_to_cartesian (tp);

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

Файл calc_8.m для восьми точек:

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = 1;

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

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

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

Файл calc_9.m для девяти точек:

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = 0.5;

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

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

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

Файл calc_10.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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-31, 'TolX', 1e-17, ...
    'MaxIter', 20000, 'MaxFunEvals', 10000);
w = fminsearch (@(w) target_func (w), w0, fmso);
disp (target_func (w))

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)
 

Файл calc_11.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = 1;

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

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

calc_plot_func (coords)
view ([-18 15])
%calc_save_image (num)
 

Файл calc_12.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = 1;

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

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

calc_plot_func (coords)
view ([-18 15])
%calc_save_image (num)
 

Файл calc_13.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

points_func = @ (w) [
    w(1)        0
    w(1)        pi / 2
    w(2)        pi / 4
    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.7 1.4 2.1];

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

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

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

Файл calc_14.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = [1.0 0.6 1.4];

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

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

calc_plot_func (coords)
view ([20 15])
%calc_save_image (num)
 

Файл calc_16.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = [0.7 1.3];

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

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

calc_plot_func (coords)
view ([-20 15])
%calc_save_image (num)
 

Файл calc_17.m

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

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

w0 = [0.58 0.72 1.29 0.84 1.61 1.82 2.25 0.68];

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

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

calc_plot_func (coords)
view ([80 15])
%calc_save_image (num)
 

[Евгений Машеров]

Ещё попытка в эмпирический алгоритм. Набрасываем точки случайно. Затем вычисляем расстояния между точками. Находим минимальное расстояние. Затем для каждой из двух точек, расстояние между которыми минимально, находим три ближайшие и строим окружности через них. Точку, для которой радиус окружности больше, перемещаем в центр окружности. Пересчитываем расстояния. И т.д. Условие остановки непонятно, возможно, будет зацикливание.

[B@R5uk]

Евгений Машеров в сообщении #1681388 писал(а):

находим три ближайшие и строим окружности через них.

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

Но в целом идея хорошая. Ближайшие точки надо расталкивать. А перемещаемую точку (назовём её точкой A) нужно двигать (каким-либо образом) до тех пор, пока она не будет "зажата" со всех сторон. Математически "зажатость со всех сторон" можно формализовать следующим образом. Для точек, которые "зажимают" точку A (что может выражаться, например, тем, что расстояния до них меньше минимального плюс некоторый порог) строятся радиус-векторы, выходящие из точки A. Эти векторы проецируются на касательную плоскость для точки A, и проекции нормируются на единицу. Дальше проверяется, что не существует вектора, все скалярные произведения на который всех остальных векторов будут положительны. В противном случае, концы всех векторов будут лежать по одну сторону от прямой, перпендикулярной заданному вектору. Другими словами, что точка A на самом деле не будет "зажата". Технически матрица скалярных произведений находится как произведение матрицы координат этих векторов на себя транспонированную (со свёрткой по размерности векторов), поскольку вектора уже отнормированы на единицу. Получается квадратная матрица с размером, равным количеству векторов. Заменяем все отрицательные числа нулями, а все неотрицательные — единицами. Если в каком-то столбце (или строке) все числа единицы (произведение элементов столбца или строки prod в MATLAB), то точка не является зажатой.

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

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

%   topic160089.html
%   B@R5uk, 2025, April

clc
clearvars
format compact
format long

num = 15;
step = 0.2;
thsh_mult = 2.3;
rand_mult = 0.3;
step_mult = 0.85;
step_end = 1e-12;
skip_frames = 50;

coords = randn (num, 3);
coords = coords ./ repmat (sqrt (sum (coords .^ 2, 2)), 1, 3);
coords_new = zeros (num, 3);

flags = zeros (num);
n = 0;

subplot (111)
plot3 (0, 0, 0)
while step_end < step
    n = n + 1;
    dists = repmat (permute (coords, [3, 1, 2]), num, 1);
    dists = sqrt (sum ((permute (dists, [2, 1, 3]) - dists) .^ 2, 3));
    dists = dists + diag (inf (num, 1));
    min_dist = min (dists (:));
    
    flags = zeros (num);
    flag_step_down = 1;
    for k = 1 : num
        cur_dists = dists (k, :);
        cur_min_dist = min (cur_dists);
        addr = min_dist + thsh_mult * step > cur_dists;
        num_nb = sum (addr);
        if 2 > num_nb
            flag_step_down = 0;
            min_pos = find (cur_min_dist == cur_dists, 1, 'first');
            vect = coords (k, :) - coords (min_pos, :);
            vect = vect / sqrt (sum (vect .^ 2));
            vect = vect - coords (k, :) * vect' * coords (k, :);
            vect = vect / sqrt (sum (vect .^ 2));
            vect = 0.5 * vect * step + coords (k, :);
            coords_new (k, :) = vect / sqrt (sum (vect .^ 2));
        elseif 2 == sum (addr)
            flag_step_down = 0;
            vals = cur_dists (addr);
            vects = coords ([k k], :) - coords (addr, :);
            vects = vects ./ repmat (vals', 1, 3);
            vects = vects - repmat (vects * coords (k, :)', 1, 3) .* coords ([k k], :);
            vects = vects ./ repmat (sqrt (sum (vects .^ 2, 2)), 1, 3);
            vect = (vals - min_dist) * vects;
            vect = vect / sqrt (sum (vect .^ 2));
            vect = 0.5 * vect * step + coords (k, :);
            coords_new (k, :) = vect / sqrt (sum (vect .^ 2));
        else
            vals = cur_dists (addr);
            vects = repmat (coords (k, :), num_nb, 1) - coords (addr, :);
            vects = vects ./ repmat (vals', 1, 3);
            vects = vects - repmat (vects * coords (k, :)', 1, 3) .* repmat (coords (k, :), num_nb, 1);
            vects = vects ./ repmat (sqrt (sum (vects .^ 2, 2)), 1, 3);
            if sum (prod (double (0 < vects * vects')))
                flag_step_down = 0;
            end
            vect = (vals - min_dist) * vects;
            vect = vect * step + coords (k, :);
            coords_new (k, :) = vect / sqrt (sum (vect .^ 2));
        end
        flags (k, cur_min_dist + min (5 * thsh_mult * step, 0.05 * min_dist) > cur_dists) = 1;
    end
    flags = flags + flags';
    if flag_step_down
        step = step_mult * step;
    else
        coords_rand = randn (num, 3);
        coords_rand = coords_rand ./ repmat (sqrt (sum (coords_rand .^ 2, 2)), 1, 3);
        coords = coords_new + coords_rand * step * rand_mult;
        coords = coords ./ repmat (sqrt (sum (coords .^ 2, 2)), 1, 3);
    end
    
    if ~mod (n, skip_frames)
        [az, el] = view;
        plot3 (coords (:, 1), coords (:, 2), coords (:, 3), 'o')
        hold on
        for k = 1 : num
            addr = find (flags (:, k) == 1);
            addr = [k * ones(numel (addr), 1), addr]'; %#ok<AGROW>
            plot3 (coords (addr, 1), coords (addr, 2), coords (addr, 3))
            addr = find (flags (:, k) > 1);
            addr = [k * ones(numel (addr), 1), addr]'; %#ok<AGROW>
            plot3 (coords (addr, 1), coords (addr, 2), coords (addr, 3), 'LineWidth', 3)
        end
        hold off
        axis vis3d
        xlim ([-1, 1])
        ylim ([-1, 1])
        zlim ([-1, 1])
        title ([n, step, min_dist, 2 * asind(min_dist / 2)])
        view (az, el)
        drawnow
    end
end
disp ('========')
disp (num)
disp ('--------')
disp (coords)
disp ('--------')
disp (min_dist)
disp ('--------')
disp (2 * asind (min_dist / 2))
 

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

Ну и coup de grace этому подходу: критерий "зажатости" точек не является критерием локальной оптимальности. Да, иногда он работает, но не всегда. Именно по этой причине я ищу аналитический метод проверки оптимальности и, если это возможно, вычисления аналога "градиента" в направлении которого целевая функция будет увеличиваться.

[B@R5uk]

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

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

Забыл сказать, что контр-пример с десятью точками приведёт на первой странице.

[мат-ламер]

B@R5uk
Работу Sloane, ссылку на которую вы выложили, не нашёл. Ну, да ладно. Буду ориентироваться на ваши таблицы. Может вас заинтересует препринт Hars https://www.hars.us/Papers/Numerical_Tammes.pdf . У него есть рисунки, таблицы, алгоритмы. Но оптимальные решения он как-то странно характеризует - не через минимальное расстояние между точками, а через какой-то угол. Думаю, что ничего страшного.

[B@R5uk]

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

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

Это угол, под которым видно ребро минимальной длины из центра сферы. Ну, или, что тоже самое, длина дуги, которая получается проекцией этого ребра на сферу (расстояние между точками на сфере). Длина ребра и угол связаны формулой:

$$l=2\sin\frac{\alpha}{2}$$

Картинки в статье кривые (не отображены все рёбра многогранника), но в целом очень информативно. Спасибо.

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

[B@R5uk]

ОК, я был не прав.

Размещение точек на рисунках выше является экстремальным (все три фигуры — это одна фигура, наклонённая и повёрнутая). Если взять глобально оптимальное расположение 10 точек на левом рисунке ниже:

и привязать длину самого верхнего ребра к длине оптимального ребра по формуле

$$l_{top}=l_{min}\cdot\Bigl(1 + \bigl(\sqrt{2}-1\bigr)\lambda\Bigr)$$

а затем построить график зависимости длины оптимального ребра от параметра лямбда, то на нём будет два максимума:

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

Если двигаться на графике от вертикальной красной линии вправо, то это будет соответствовать удлинению верхнего ребра фигуры. При этом разорвётся пара горизонтальных жёстких рёбер в центре фигуры, что изображено на картинке снизу посередине. В конце-концов растягивание ребра приведёт к единичному значению параметра. Это будет соответствовать тому, что ребро будет в ${}_{{}_{{}_.}}\sqrt{2}$ раз длиннее оптимального, превратившись в диагональ квадратной грани, которая получается слиянием двух прилегающих треугольных граней. Это изображено на нижнем правом рисунке. Эта конфигурация не является не в каком смысле оптимальной (даже седловой), но она имеет одну замечательную особенность: пара самых нижних точек может быть расположена двояко и это не изменит целевую функцию задачи. На рисунке альтернативное расположение точек отмечено выходящими из них красными жёсткими рёбрами.

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

Программа для расчёта графика calc_10_dist_graph.m:

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

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

clc
clearvars
format compact
format long

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

w0 = [0.7 1.0 1.7 0.7 1.7 1.0 2.6 2.6];
a = sqrt (2) - 1;

fmso = optimset ('fminsearch');
fmso = optimset (fmso, 'Display', 'final', 'TolFun', 1e-29, 'TolX', 1e-17, ...
    'MaxIter', 20000, 'MaxFunEvals', 10000);
%{
w = fminsearch (@(w) target_func (w, a), w0, fmso);
disp (target_func (w, a))

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

calc_plot_func (coords)
view ([-24 15])
%}

dist_10_max      = 1.09142629140382;
dist_10_max_next = 1.219891420230531;
xx_mult = sqrt (2) - 1;
xx_max = (dist_10_max_next / dist_10_max - 1) / xx_mult;
xx = [0 : 0.05 : 0.25, xx_max, 0.3 : 0.05 : 1];
yy = zeros (size (xx));
w = w0;
for k = 1 : numel (xx)
    a = xx_mult * xx (k);
    w = fminsearch (@(w) target_func (w, a), w, fmso);
    yy (k) = 2 * sin (pi - (w (7) + w (8)) / 2);
    disp (target_func (w, a))
    disp ([xx(k), yy(k)])
end

plot (xx_max * [1 1], [1.04 1.1], 'r')
hold on
plot ([0 1], dist_10_max * [1 1], 'r')
plot (xx, yy, 'LineWidth', 2)
hold off
ylim ([1.04, 1.1])
 

-- 09.04.2025, 19:33 --

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

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

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

Забыл сказать, что контр-пример с десятью точками приведён на первой странице.

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

[мат-ламер]

B@R5uk
А так уж и важна нам визуализация? В исходной постановке ведь речь идёт о системе точек на сфере. И ни о каком многограннике не говорится. В картинках, что я вижу у вас и у Haars много граней - четырёхугольники и пятиугольники. И не факт, что они всегда плоские.

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