Python环境下高效计算多面体解析中心的可行方案

解析中心的数学定义是严格满足所有约束下，对数障碍函数$\sum_{i=1}^{m_{ineq}} \log(b_i - a_i^T x)$的最大值点，对应中心路径上障碍参数趋于无穷的极限点，不是普通线性规划内点法迭代过程中的任意内点。

Python原生直接落地实现（适配numpy数组输入、大规模MIPLIB级问题） #

• 放弃纯稠密NumPy手写教学版原始对偶牛顿法的思路：叶荫宇教材中的算法是原理演示版本，没有集成稀疏线性代数、预条件、自适应线搜索、海森修正等工业级优化，处理千条约束以上的稀疏问题时效率会低一个数量级以上。
• 最低代码量实现：直接用凸优化建模框架构造解析中心问题，调用工业级稀疏内点求解器计算，输入可直接接numpy格式的约束矩阵，代码示例如下：

import cvxpy as cp
import numpy as np
from scipy import sparse

# 输入参数格式：
# A_ub: 不等式约束矩阵，shape=(m_ineq, n)，numpy数组
# b_ub: 不等式约束右端项，shape=(m_ineq,)，numpy数组
# A_eq: 等式约束矩阵，shape=(m_eq, n)，numpy数组，无等式约束时传空数组
# b_eq: 等式约束右端项，shape=(m_eq,)，numpy数组，无等式约束时传空数组

x = cp.Variable(A_ub.shape[1])
# 转稀疏压缩格式，适配MIPLIB问题99%以上的稀疏度，大幅降低计算开销
A_ub_sp = sparse.csr_matrix(A_ub)
constraints = [A_ub_sp @ x <= b_ub]
if A_eq.size > 0:
    A_eq_sp = sparse.csr_matrix(A_eq)
    constraints.append(A_eq_sp @ x == b_eq)

# 解析中心对应最大化所有不等式约束的松弛量对数和
prob = cp.Problem(
    cp.Maximize(cp.sum(cp.log(b_ub - A_ub_sp @ x))),
    constraints
)
# 选CLARABEL求解器：开源、支持稀疏矩阵、原生支持对数目标，关闭预处理避免修改原多面体结构
prob.solve(
    solver=cp.CLARABEL,
    presolve=False,
    scale=False,
    verbose=False,
    max_iter=200,
    tol_gap_abs=1e-9,
    tol_feas=1e-9
)
analytic_center = x.value

注意必须关闭求解器的预处理、自动缩放、交叉步骤，否则求解器会自动约简冗余约束、缩放变量，导致结果和原多面体的真实解析中心出现偏差。

• 若需要自己实现迭代逻辑进一步提速：核心优化点是全程用scipy.sparse存储约束矩阵，牛顿步求解用sksparse.cholmod的稀疏Cholesky分解做线性方程求解，不要构造稠密海森矩阵做求逆，单步迭代速度可以比纯NumPy稠密实现快10~100倍，收敛判据直接卡对数障碍函数的梯度无穷范数小于1e-8即可，不要用线性规划的最优间隙作为判停标准。

导出MPS调用第三方工具的方案 #

• 模型导出：用scipy.io或者pulp、cvxpy的模型导出接口，将numpy格式的约束转成标准MPS文件，导出时关闭所有自动预处理、约束约简选项，保证多面体描述和原始输入完全一致。
• 求解选项：
  • 用Pyomo读入MPS文件，构造对数障碍最大化目标，调用IPOPT求解，IPOPT的稀疏牛顿迭代对十万级约束规模的问题处理效率很高，收敛精度可以稳定到1e-10以上。
  • 调用专门的多面体内点计算工具LOPOCS，读入MPS后可直接输出解析中心，不需要额外建模，对大规模LP松弛问题的适配性优于通用LP求解器。

此前测试方案失效的核心原因 #

• scipy linprog的interior-point方法返回结果有偏差：该方法是为线性规划最优解设计的，内点路径跟踪只是寻找最优解的中间过程，不会收敛到障碍函数的全局极小点，判停条件也没有针对解析中心的精度要求设置。
• scipy linprog的highs-ipm方法返回顶点解：HiGHS内点法默认开启crossover（交叉）步骤，会将收敛到最优面的内点投影到多面体顶点，即使目标函数为零，也会返回一个可行顶点而非内部的解析中心。

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