我正在寻找一种算法来优化球体表面上的 N 个点分布(不是随机的,我说的是 Tammes problem 意义上的最优分布)。
我想在 C++ 中这样做,但问题不是语言问题,而是算法问题。这是我目前使用的算法(在 python 中):
# Function to distribute N points on the surface of a sphere
# (source: http://www.softimageblog.com/archives/115)
def uniform_spherical_distribution(N):
pts = []
inc = math.pi * (3 - math.sqrt(5))
off = 2 / float(N)
for k in range(0, int(N)):
y = k * off - 1 + (off / 2)
r = math.sqrt(1 - y*y)
phi = k * inc
pts.append([math.cos(phi)*r, y, math.sin(phi)*r])
return pts
对于大量的点,视觉结果似乎相当不错,但对于少量的点(2 到 10 之间),结果似乎一点也不好(半球中有更多点)。有没有一种方法可以增强算法,使其在点数较少的情况下也可以...(答案可以是 C、C++ 或 Python)。
如果您只需要较少的点数,蛮力可能会奏效。我所说的蛮力是指预先计算 2..10
的解决方案。如果 N<=10
,则在数组中查找它们.
对于小N
就像那样,应该可以解决数学优化问题,以最佳方式分配您的点数。你需要 multistart 但对于小 N
这应该不是问题。这是我在 AMPL 中的尝试:
param N;
var x{1..N};
var y{1..N};
var z{1..N};
var d;
param xbest{1..N};
param ybest{1..N};
param zbest{1..N};
param dbest;
maximize obj: d;
all_points_on_the_sphere{i in 1..N}:
x[i]^2 + y[i]^2 + z[i]^2 = 1;
all_pairs{i in 1..N, j in 1..N : i<j}:
(x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2 >= d;
fix_first_x: x[1] = 1;
fix_first_y: y[1] = 0;
fix_first_z: z[1] = 0;
fix_second_z: z[2] = 0;
#############################################
option show_stats 1;
option presolve 10;
option substout 1;
option var_bounds 2;
#option nl_comments 0;
#option nl_permute 0;
option display_precision 0;
option solver "/home/ali/ampl/ipopt";
for {k in 2..10} {
let N := k;
solexpand _con;
let dbest := -1.0;
# Multistart
for {1..2000} {
for {j in 1.._snvars}
let _svar[j] := Uniform(-1, 1);
let d := N;
solve;
if (solve_result_num < 200 and d > dbest) then {
let dbest := d;
for {i in 1..N} {
let xbest[i] := x[i];
let ybest[i] := y[i];
let zbest[i] := z[i];
}
}
}
print "@@@";
printf "N=%d, min distance %6f\n", N, sqrt(dbest);
for {i in 1..N}
printf "(%9.6f, %9.6f, %9.6f)\n", xbest[i], ybest[i], zbest[i];
}
在我的机器上运行这个脚本花了 5 分钟,解决方案是:
N=2, min distance 2.000000
( 1.000000, 0.000000, 0.000000)
(-1.000000, 0.000000, 0.000000)
N=3, min distance 1.732051
( 1.000000, 0.000000, 0.000000)
(-0.500000, 0.866025, 0.000000)
(-0.500000, -0.866025, 0.000000)
N=4, min distance 1.632993
( 1.000000, 0.000000, 0.000000)
(-0.333333, -0.942809, 0.000000)
(-0.333333, 0.471405, -0.816497)
(-0.333333, 0.471405, 0.816497)
N=5, min distance 1.414214
( 1.000000, 0.000000, 0.000000)
(-0.208840, 0.977950, 0.000000)
(-0.000000, 0.000000, 1.000000)
(-0.212683, -0.977121, 0.000000)
( 0.000000, 0.000000, -1.000000)
N=6, min distance 1.414214
( 1.000000, 0.000000, 0.000000)
(-1.000000, 0.000000, 0.000000)
( 0.000000, -0.752754, -0.658302)
( 0.000000, 0.752754, 0.658302)
( 0.000000, 0.658302, -0.752754)
( 0.000000, -0.658302, 0.752754)
N=7, min distance 1.256870
( 1.000000, 0.000000, 0.000000)
(-0.688059, -0.725655, 0.000000)
( 0.210138, -0.488836, -0.846689)
( 0.210138, -0.488836, 0.846688)
(-0.688059, 0.362827, 0.628435)
(-0.688059, 0.362827, -0.628435)
( 0.210138, 0.977672, 0.000000)
N=8, min distance 1.215563
( 1.000000, 0.000000, 0.000000)
( 0.261204, -0.965284, 0.000000)
( 0.261204, 0.565450, 0.782329)
(-0.783612, -0.482642, -0.391165)
( 0.261204, -0.199917, -0.944355)
( 0.261204, 0.882475, -0.391165)
(-0.783612, 0.599750, 0.162026)
(-0.477592, -0.399834, 0.782329)
N=9, min distance 1.154701
( 1.000000, 0.000000, 0.000000)
(-0.500000, -0.866025, 0.000000)
( 0.333333, -0.577350, -0.745356)
(-0.500000, 0.866025, 0.000000)
(-0.666667, -0.000000, 0.745356)
(-0.666667, 0.000000, -0.745356)
( 0.333333, -0.577350, 0.745356)
( 0.333333, 0.577350, -0.745356)
( 0.333333, 0.577350, 0.745356)
N=10, min distance 1.091426
( 1.000000, 0.000000, 0.000000)
(-0.605995, 0.795469, 0.000000)
( 0.404394, 0.816443, 0.412172)
(-0.664045, -0.746251, -0.046407)
( 0.404394, -0.363508, -0.839242)
(-0.664045, 0.002497, 0.747688)
(-0.605995, 0.046721, -0.794096)
( 0.404394, -0.908368, 0.106452)
( 0.255933, 0.703344, -0.663179)
( 0.404394, -0.159620, 0.900548)
通过查看数字,很明显可以通过分析计算某些解决方案(我认识 sqrt(2) 和 sqrt(3) 等)。我相信 N=2、4 和 6 的解是直线 ([-1,0,0], [1,0,0])、四面体、八面体。
没有strong保证以上是点的最佳分布。非线性求解器可能陷入局部最优;作为N
增长,局部最优的数量也增长。
您可以将上述解决方案放在一个数组中,然后在 Python、C++ 或您使用的任何语言中使用它们。