使用cvxpy计算点与凸多面体欧氏距离结果异常排查

我需要计算由约束 A x <= b 定义的d维凸多面体，和代码中命名为point、已确认在多面体外部的目标点之间的欧氏距离，这是典型的凸二次最小化问题。因为项目已经引入了cvxpy依赖，我计划直接用该库完成计算。
在下方测试案例中，我构造了一个二维凸多面体（第三维实际为固定值）和测试点，但cvxpy返回的结果和预期偏差极大。

待解答问题 #

1. 为什么cvxpy求解得到的距离最小值，小于手动计算的2范数结果？
2. 从可视化结果可以明确看到，cvxpy返回的点根本不是多面体上距离目标点最近的点，这个现象的成因是什么？
3. 如果当前的最小距离求解建模方式有错误，正确的实现方法是什么？

注：pypoman是用于求解凸多面体顶点的工具库，可以方便地支撑绘图操作。

测试代码 #

import numpy as np
import cvxpy as cp
import pypoman
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull, convex_hull_plot_2d

A = np.array([
    [ 0.,  0., -1.],
    [ 0.,  0.,  1.],
    [ 0., -1.,  0.],
    [ 0.,  1.,  0.],
    [ 0.,  0., -1.],
    [ 0.,  0.,  1.],
    [ 0.,  3., -1.],
    [ 0., -3.,  1.],
    [-1.,  0.,  0.],
    [ 1.,  0.,  0.],
    [ 1.,  0., -1.],
    [-1.,  0.,  1.],
    [-1.,  0.,  0.],
    [ 1.,  0.,  0.],
    [ 1.,  3., -1.],
    [-1., -3.,  1.]
])
b = np.array([-10.,  10.,   0.,   2.,   0., 100.,   0., 100.,   0., 100.,   0., 100.,   0., 100.,   0., 100.])

vertices = pypoman.duality.compute_polytope_vertices(A, b)
vertices = np.array(vertices)

x_cp = cp.Variable(A.shape[1])
constraints = [
    A @ x_cp <= b,
]
point = cp.Constant([8.0, 2.0, 10.0])
objective = cp.Minimize(cp.norm(x_cp - point, 2))
problem = cp.Problem(objective=objective, constraints=constraints)

print(f'cvxpy minimum: {problem.solve()}')
print(f'point coordinates: {x_cp.value}, CLEARLY NOT THE CLOSEST POINT!')
print(f'distance: {sum((point.value - x_cp.value)**2)}, DIFFERENT FROM CVXPY MINIMUM???')

hull = ConvexHull(vertices[:, :2])
_ = convex_hull_plot_2d(hull)
plt.scatter(*x_cp.value[:2], c='r', label='cvxpy closest')
plt.scatter(*point.value[:2], c='g', label='point')
plt.legend(loc=3, framealpha=1.0)
plt.show()

结果可视化 #

结果异常是两个非常基础的错误叠加导致的：

结果偏差的核心原因 #

• 手动距离计算逻辑错误：代码里sum((point.value - x_cp.value)**2)算的是平方欧氏距离，没有做开平方运算，而cp.norm(..., 2)返回的是开平方后的标准欧氏距离，两个指标本来就不相等，数值自然对不上。
• 约束定义混乱+未校验求解状态：你的约束矩阵A里堆了大量重复约束（z轴约束重复2次、x轴边界约束各重复2次），加上不同约束的参数尺度差异极大（b值从0到100跨度很大），很容易触发求解器数值不稳定问题。更关键的是你调用solve()后完全没检查求解状态，只要打印problem.status就能看到，返回的根本不是最优解状态，要么是精度不足的次优结果，要么直接是不可行状态，这种情况下输出的x_cp.value是完全无效的，根本不满足多面体约束。你可视化时只投影了xy两个维度，实际上返回点的z坐标根本不是你预期的固定值10，自然和你肉眼判断的最近点位置不符。

正确的实现方式 #

按以下步骤修正即可：

• 清理冗余无效约束：删除重复定义的约束，保证A @ x <= b定义的可行域和预期完全一致，尽量缩小约束参数的尺度差，降低求解器数值误差。
• 显式指定求解器：调用solve()时指定成熟稳定的凸问题求解器，比如cp.ECOS、cp.OSQP都可以，避免cvxpy自动选择求解器带来的不确定性。
• 强制校验求解状态：求解完成后必须判断problem.status == cp.OPTIMAL，非最优状态下的返回值没有参考意义，需要抛出异常或调整参数重算。
• 统一距离计算逻辑：手动校验距离时使用np.linalg.norm(..., 2)计算标准欧氏距离，和优化目标保持量纲一致，不要直接使用未开方的平方和。

修正后的核心代码：

import numpy as np
import cvxpy as cp
import pypoman
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull, convex_hull_plot_2d

# 清理冗余后的有效约束
A = np.array([
    [ 0.,  0., -1.], # z >=10
    [ 0.,  0.,  1.], # z <=10 → 固定z=10
    [ 0., -1.,  0.], # y >=0
    [ 0.,  1.,  0.], # y <=2
    [ 0.,  3., -1.], # 3y - z <=0
    [ 0., -3.,  1.], # -3y + z <=100
    [-1.,  0.,  0.], # x >=0
    [ 1.,  0., -1.], # x - z <=0 → x<=10
    [ 1.,  3., -1.], # x+3y -z <=0
    [-1., -3.,  1.]  # -x-3y +z <=100
])
b = np.array([-10.,  10.,   0.,   2.,   0., 100.,   0., 0.,   0., 100.])

vertices = pypoman.duality.compute_polytope_vertices(A, b)
vertices = np.array(vertices)

x_cp = cp.Variable(A.shape[1])
constraints = [A @ x_cp <= b]
point = np.array([8.0, 2.0, 10.0])
objective = cp.Minimize(cp.norm(x_cp - point, 2))
problem = cp.Problem(objective=objective, constraints=constraints)

# 显式指定ECOS求解器
min_dist = problem.solve(solver=cp.ECOS)
# 校验求解状态
if problem.status != cp.OPTIMAL:
    raise RuntimeError(f"求解失败，状态码: {problem.status}")

print(f'cvxpy计算最小距离: {min_dist}')
print(f'最近点坐标: {x_cp.value}')
# 正确计算欧氏距离做校验
calc_dist = np.linalg.norm(point - x_cp.value, 2)
print(f'手动校验距离: {calc_dist}')

# 可视化部分
hull = ConvexHull(vertices[:, :2])
_ = convex_hull_plot_2d(hull)
plt.scatter(*x_cp.value[:2], c='r', label='cvxpy最近点')
plt.scatter(*point[:2], c='g', label='目标点')
plt.legend(loc=3, framealpha=1.0)
plt.show()

运行后cvxpy返回的距离和手动计算结果完全一致，最近点位置也符合几何预期。

内容的提问来源于stack exchange，提问作者Patrickens