我是一个编程新手，尤其是python，但我正在努力学习它，我发现到目前为止它非常迷人。
我有一个30个固定坐标x和y的列表。
x = np.array([[13,10,12,13,11,12,11,13,12,13,14,15,15,16,18,2,3,4,6,9,1,3,6,7,8,10,12,11,10,30]])
y = np.array([[12,11,10,9,8,7,6,6,7,8,11,12,13,15,14,18,12,11,10,13,15,16,18,17,16,15,14,13,12,3]])
我想找到一个优化的（集中的）位置，通过找到最小距离，可以连接到10个固定坐标的最大值。
我厌倦了使用优化算法，但我只能得到一个坐标列表的最佳位置（即，我不能将位置限制为10个坐标）。
任何帮助都将不胜感激。
def objective(x):
px = x[0]
py = x[1]
wdistance = sum(((px - wx)**2 + (py - wy)**2)**0.5)
tdistance = (((px - tx)**2 + (py - ty)**2)**0.5)
cost = wdistance * 2.5 + tdistance * 5
return cost
# initial estimate of the centralized location
x0 = [10,10]
# boundary values for the centralized locations
bnds = ((0,15), (0,15))
sol = minimize(objective, x0, bounds=bnds, method='SLSQP')
print(sol)
最佳答案
正如我在评论中指出的，我认为这个问题是NP-hard的，因此不容易解决。
这是一个具有某种基数约束的问题，这意味着，通常很难/不可能使用连续优化（假定为scipy.optimize中的几乎所有内容）。
此外，最小化最小距离并不是特别平滑，scipy中的几乎所有内容都假设了这一点。但是可以重新格式化（至少在支持约束的情况下！）.
如果没有基数约束（这个问题的离散性质的核心），它就会陷入Smallest Enclosing Circle Problem问题，很容易在线性时间内求解。更一般的wiki-article。但这对给定的问题没有多大帮助。
这是一个演示，可能不是针对心脏病患者（取决于背景知识）：
基于Mixed Integer Second Order Cone Programming
将求解无限时间（忽略数值问题；FP数学）的最优性（ε逼近）
不超过NP硬度；但可能有助于开发数据中的某种结构
对基于欧氏距离（范数）的问题感觉很自然
使用cvxpy作为建模工具
使用ECOS_BB作为解算器
警告：这是概念证明/玩具求解器（我们将在结果中看到）
MISOCP/MICP是一个活跃的研究领域，但非商业专业解决方案的数量是稀少的！
代替ECOS_BB，人们可能应该使用Pajarito，可能是Bonmin或广告（Gurobi，CPLEX，Mosek和co.）
Pajarito看起来很有趣，但是（imho）与python一起使用并不有趣；因此请将我的代码看作是一个带有演示求解器的演示！
下面的代码将解决20个问题中的10个（前20个点；其中一个副本被丢弃以获得更好的可视化！）
它将失败的完全问题（相当好的近似将被返回；但不是最佳的！）.I责怪ECOS_BB（但同时非常感谢韩旺的代码！）你的问题肯定会很容易通过更好的选择来解决。但不要指望任何解算器能得到1000点。
代码：
import numpy as np
import cvxpy as cvx
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist
""" (Reduced) Data """
# removed a duplicate point (for less strange visualization)
x = np.array([13,10,12,13,11,12,11,13,13,14,15,15,16,18,2,3,4,6,9,1])#,3,6,7,8,10,12,11,10,30])
y = np.array([12,11,10,9,8,7,6,6,8,11,12,13,15,14,18,12,11,10,13,15])#,16,18,17,16,15,14,13,12,3])
N = x.shape[0]
M = 10
""" Mixed-integer Second-order cone problem """
c = cvx.Variable(2) # center position to optimize
d = cvx.Variable(N) # helper-var: lower-bound on distance to center
d_prod = cvx.Variable(N) # helper-var: lower-bound on distance * binary/selected
b = cvx.Bool(N) # binary-variables used for >= M points
dists = pdist(np.vstack((x,y)).T)
U = np.amax(dists) # upper distance-bound (not tight!) for bigM-linearization
helper_expr_c_0 = cvx.vec(c[0] - x)
helper_expr_c_1 = cvx.vec(c[1] - y)
helper_expr_norm = cvx.norm(cvx.hstack(helper_expr_c_0, helper_expr_c_1), axis=1)
constraints = []
constraints.append(cvx.sum_entries(b) >= M) # M or more points covered
constraints.append(d >= helper_expr_norm) # lower-bound of distances
constraints.append(d_prod <= U * b) # linearization of product
constraints.append(d_prod >= 0) # """
constraints.append(d_prod <= d) # """
constraints.append(d_prod >= d - U * (1-b)) # """
objective = cvx.Minimize(cvx.max_entries(d_prod))
problem = cvx.Problem(objective, constraints)
problem.solve(verbose=False)
print(problem.status)
# Visualization
center = np.array(c.value.flat)
tmp = np.array(np.round(b.value.T.flat), dtype=int)
selected_points = np.where(tmp == 1)[0]
non_selected_points = np.where(tmp == 0)[0]
ax = plt.subplot(aspect='equal')
ax.set_title('Optimal solution radius: {0:.2f}'.format(problem.value))
ax.scatter(x[non_selected_points], y[non_selected_points], c='blue', s=70, alpha=0.9)
ax.scatter(x[selected_points], y[selected_points], c='magenta', s=70, alpha=0.9)
ax.scatter(c[0].value, c[1].value, c='red', s=70, alpha=0.9)
circ = plt.Circle((center[0], center[1]), radius=problem.value, color='g', fill=True, alpha=0.1)
ax.add_patch(circ)
plt.tight_layout()
plt.show()
输出：
optimal
绘图：
julia based
编辑
我将上述代码移植到julia，以便能够使用上述解算器Pajarito。尽管过去茱莉亚为我做的节目不多，但现在已经开始运作了：
使用
上面基于cvxpy的代码几乎1:1的映射！
使用3个开源解决方案的组合：
Convex.jl作为MISOCP/MICP解算器
使用Pajarito作为内部MIP解算器（Coin-CBC通常要好得多；但由于Pajarito中可能存在的错误，不适合我）
使用GLPK作为内凸求解器
从上面重用基于matplotlib的绘图代码；手工复制的结果值：-）
这种方法是解决你的全部问题（再次我放弃了复制点）的最佳化！
它输出一些数值不稳定警告，但声称结果是最优的！（Pajarito仍然是相当新的研究软件，我没有调整众多的选项；我们使用的开源解决方案与商业替代方案相比没有那么强大）
我对结果非常满意，我很高兴我给了Julia/凸面.jl/Pajarito一个尝试！
代码：
using Convex, Pajarito, ECOS, GLPKMathProgInterface
# DATA
x = [13 10 12 13 11 12 11 13 13 14 15 15 16 18 2 3 4 6 9 1 3 6 7 8 10 12 11 10 30]
y = [12 11 10 9 8 7 6 6 8 11 12 13 15 14 18 12 11 10 13 15 16 18 17 16 15 14 13 12 3]
N = size(x)[2]
M = 10
# MISOCP
c = Variable(2)
d = Variable(N)
d_prod = Variable(N)
b = Variable(N, :Bin)
U = 100 # TODO
# Objective
problem = minimize(maximum(d_prod))
# Constraints
problem.constraints += sum(b) >= M
problem.constraints += d_prod <= U*b
problem.constraints += d_prod >= 0
problem.constraints += d_prod <= d
problem.constraints += d_prod >= d - U * (1-b)
for i = 1:N
problem.constraints += d[i] >= norm([c[1] - x[i], c[2] - y[i]]) # ugly
end
# Solve
mip_solver_drives = false
log_level = 2
rel_gap = 1e-8
mip_solver = GLPKSolverMIP()
cont_solver = ECOSSolver(verbose=false)
solver = PajaritoSolver(
mip_solver_drives=mip_solver_drives,
log_level=log_level,
rel_gap=rel_gap,
mip_solver=mip_solver,
cont_solver=cont_solver,
init_exp=true,
)
# option taken from https://github.com/JuliaOpt/Pajarito.jl/blob/master/examples/gatesizing.jl
# not necessarily optimal
solve!(problem, solver)
# Print out results
println(problem.status)
println(problem.optval)
println(b.value)
println(c.value)
输出：
Problem dimensions:
variables | 120
constraints | 263
nonzeros in A | 466
Cones summary:
Cone | Count | Min dim. | Max dim.
Second order | 29 | 3 | 3
Variable types:
continuous | 91
binary | 29
Transforming data... 1.10s
Creating conic subproblem... 0.52s
Building MIP model... 2.35s
Solving conic relaxation... 1.38s
Starting iterative algorithm
0 | +Inf | +3.174473e-11 | Inf | 6.434e+00
Warning: numerical instability (primal simplex, phase I)
Iter. | Best feasible | Best bound | Rel. gap | Time (s)
1 | +2.713537e+00 | +2.648653e+00 | 2.391e-02 | 1.192e+01
Warning: numerical instability (primal simplex, phase I)
Warning: numerical instability (primal simplex, phase I)
Warning: numerical instability (dual simplex, phase II)
... some more ...
2 | +2.713537e+00 | +2.713537e+00 | 5.134e-12 | 1.341e+01
Iterative algorithm summary:
- Status = Optimal
- Best feasible = +2.713537e+00
- Best bound = +2.713537e+00
- Relative opt. gap = 5.134e-12
- Total time (s) = 1.34e+01
Outer-approximation cuts added:
Cone | Relax. | Violated | Nonviol.
Second order | 58 | 32 | 24
0 numerically unstable cone duals encountered
Distance to feasibility (negative indicates strict feasibility):
Cone | Variable | Constraint
Linear | NA | 5.26e-13
Second order | NA | 2.20e-11
Distance to integrality of integer/binary variables:
binary | 0.00e+00
Optimal
2.7135366682753825
[1.0; 1.0; 1.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 1.0; 1.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 1.0; 1.0; 0.0]
[12.625; 11.6875]
可视化：
ECOS
在所有点中再绘制一个select-23：