如何基于多元有理不等式组生成多维网格？（Python工具实现）

有6个有理但不一定线性的多元不等式组，形式示例：

• 1 + x1*x2 > x3
• x4 > (1 + x5)*x6
• x3/x4 > 1 - x2

需要生成满足所有不等式的(x1,x2,x3,x4,x5,x6,x7)多维网格，需满足两个条件：

1. 所有点均符合不等式约束；
2. 合理覆盖约束界定的子空间（例如由x+y≤6、x-y≥1、2≤y≤4定义的网格为((3,2),(4,2),(3,3),(4,3))）。

已知可通过先求解各变量极值、生成笛卡尔积母网格再逐一筛选，但该方法会产生大量无效点，优先寻求SymPy实现方案，其他Python包也可。疑问：SymPy无法直接化简多元不等式，但仅需生成网格是否有间接实现方式？有人提及SciPy的linprog可化简多元不等式，但它仅支持线性不等式，且不清楚如何从优化问题转化到不等式化简问题。

一、SymPy间接实现思路 #

1. 分步约束下的变量采样 #

无需直接化简整个不等式组，而是按变量依赖关系分层采样：

• 先确定无依赖变量的初始合理采样区间（比如无约束的变量可先设为[-10,10]，后续根据其他约束调整）；
• 对有依赖的变量，基于已采样的变量值动态计算可行范围。例如：
  • 确定x1、x2的值后，通过1 + x1*x2 > x3得到x3的上界，再结合其他涉及x3的约束进一步缩小可行区间；
  • 对x4，结合x4 > (1 + x5)*x6和变形后的x3/x4 > 1 - x2（需注意1 - x2的正负性，分情况处理）确定其可行区间。

2. 利用单变量不等式求解缩小范围 #

针对每个变量，固定其他变量的部分值后，用sympy.solve_univariate_inequality求解单变量不等式：

from sympy import symbols, solve_univariate_inequality

x1, x2, x3 = symbols('x1 x2 x3')
# 假设x1=2, x2=3，求解x3的可行范围
expr = 1 + x1*x2 > x3
sol = solve_univariate_inequality(expr.subs({x1:2, x2:3}), x3)
# sol会给出x3 < 7，后续可在该区间内采样x3

这种逐变量、带条件的求解方式，能大幅减少无效采样点，避免生成全量笛卡尔积。

二、其他Python包方案 #

1. SciPy的局限性说明 #

SciPy的linprog仅支持线性不等式约束，且用于求解线性规划最优解，无法处理非线性不等式的化简或网格生成。若你的不等式组包含非线性项（如x1*x2、x3/x4），linprog完全不适用，无需考虑该方向。

2. 拉丁超立方采样+约束验证 #

用拉丁超立方采样(LHS) 替代笛卡尔积，在变量初始大致范围内生成均匀分布的样本点，再筛选符合约束的点：

from scipy.stats.qmc import LatinHypercube
import numpy as np

# 定义变量数量和初始范围（示例：每个变量初始范围[-10,10]）
n_vars = 7
l_bounds = [-10]*n_vars
u_bounds = [10]*n_vars

# 生成拉丁超立方样本
sampler = LatinHypercube(n_vars)
sample = sampler.random(n=1000)  # 生成1000个样本
# 转换到实际范围
sample = l_bounds + sample * (np.array(u_bounds) - np.array(l_bounds))

# 定义约束验证函数
def check_constraints(point):
    x1, x2, x3, x4, x5, x6, x7 = point
    # 处理分母为0的情况
    if x4 == 0:
        return False
    cond1 = 1 + x1*x2 > x3
    cond2 = x4 > (1 + x5)*x6
    cond3 = x3/x4 > 1 - x2
    # 补充其余3个约束判断
    return cond1 and cond2 and cond3 and ...  # 所有约束需同时满足

# 筛选符合条件的点
valid_points = [p for p in sample if check_constraints(p)]

该方法生成的样本点比笛卡尔积更均匀，能以更少的点覆盖可行域，减少无效计算。

3. Pyomo约束建模+采样 #

Pyomo是优化建模工具，可规范定义非线性约束，结合采样生成可行点：

from pyomo.environ import ConcreteModel, Var, Constraint

# 定义模型
model = ConcreteModel()
model.x1 = Var(bounds=(-10,10))
model.x2 = Var(bounds=(-10,10))
# 定义其余变量
model.x7 = Var(bounds=(-10,10))

# 添加约束
model.c1 = Constraint(expr=1 + model.x1*model.x2 > model.x3)
model.c2 = Constraint(expr=model.x4 > (1 + model.x5)*model.x6)
model.c3 = Constraint(expr=model.x3/model.x4 > 1 - model.x2)
# 补充其余约束

# 可结合随机采样验证约束，或调用求解器辅助生成可行点

Pyomo的优势是能更系统地管理约束，方便后续调整或扩展。

三、关键注意事项 #

• 处理分式约束时，必须额外判断分母是否为0，避免计算错误；
• 若变量范围不确定，可通过极值优化确定每个变量的大致边界：对目标变量，在其他变量取可行值的前提下，分别最大化/最小化该变量，得到上下界以缩小初始采样范围；
• 网格的“合理覆盖”可通过控制采样密度实现，在可行域内对变量等间隔采样，或在约束变化剧烈的区域增加采样点。

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