如何在Python中从线性不等式系统生成随机可行点？

问题背景 #

需要从无限域线性不等式约束系统中生成n个随机可行点，现有两种尝试均失败：

• 使用scipy.optimize.linprog仅能得到单个最优解，添加噪声也无法生成不同可行点；
• 尝试通过半空间交集求凸包再生成凸组合点时，触发HalfspaceIntersectionError报错。

可行实现方法 #

方法一：随机采样+约束验证（简单直接） #

核心思路：在变量预设范围内生成随机点，验证是否满足所有约束条件，筛选出可行点。适合约束复杂度较低的场景。

import numpy as np

def generate_random_feasible_points(A, b, n_points, var_bounds=None):
    # 默认变量范围设为0到1，可根据实际需求调整
    if var_bounds is None:
        var_bounds = [(0, 1) for _ in range(A.shape[1])]
    
    feasible_points = []
    while len(feasible_points) < n_points:
        # 在变量范围内生成随机点
        point = np.array([np.random.uniform(low, high) for low, high in var_bounds])
        # 验证所有约束（加1e-8容差处理浮点计算误差）
        if np.all(A @ point <= b + 1e-8):
            feasible_points.append(point)
    return np.array(feasible_points)

# 测试用例
A_test = np.array([[1, 1, 0, 0, 0], [-1, -1, 0, 0, 0], [0, 0, 1, 1, 0], [0, 0, -1, -1, 0], [0, 0, 0, 0, 1], [0, 0, 0, 0, -1], [1, 1, 0, 0, 0], [-1, -1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, -1, 0, 0], [0, 0, 0, 1, 1], [0, 0, 0, -1, -1], [1, 1, 0, 0, 0], [-1, -1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, -1, 0, 0], [0, 0, 0, 1, 1], [0, 0, 0, -1, -1], [1, 1, 1, 1, 1], [-1, -1, -1, -1, -1]])
b_test = np.array([0.42000000000000004, -0.42000000000000004, 0.35, -0.35, 0.23, -0.23, 0.42000000000000004, -0.42000000000000004, 0.18, -0.18, 0.4, -0.4, 0.42000000000000004, -0.42000000000000004, 0.18, -0.18, 0.4, -0.4, 1, -1])

# 生成5个随机可行点
points = generate_random_feasible_points(A_test, b_test, 5)
print(points)

注意：若可行域占变量范围比例极小，采样效率会降低。可先通过linprog找到几个可行点，将变量范围限定在这些点的极值区间内，缩小采样范围提升效率。

方法二：凸包顶点凸组合采样（分布更均匀） #

针对之前半空间交集报错的问题，先处理约束冗余并确保初始点严格在可行域内部，再通过凸包顶点的凸组合生成随机点，适合有界凸多面体可行域场景。

import numpy as np
from scipy.spatial import HalfspaceIntersection, ConvexHull
from scipy.optimize import linprog

def remove_redundant_constraints(A, b, tol=1e-8):
    # 先去除重复约束
    unique_rows = np.unique(np.hstack((A, b[:, None])), axis=0)
    A_unique = unique_rows[:, :-1]
    b_unique = unique_rows[:, -1]
    
    # 进一步去除冗余约束（判断约束是否被其他约束隐含）
    redundant_indices = []
    for i in range(len(A_unique)):
        c = -A_unique[i]
        res = linprog(c, A_ub=A_unique[np.arange(len(A_unique)) != i], b_ub=b_unique[np.arange(len(A_unique)) != i])
        if res.success and np.dot(A_unique[i], res.x) <= b_unique[i] + tol:
            redundant_indices.append(i)
    A_pruned = np.delete(A_unique, redundant_indices, axis=0)
    b_pruned = np.delete(b_unique, redundant_indices)
    return A_pruned, b_pruned

def generate_convex_hull_points(A, b, n_points):
    # 去除冗余约束，优化半空间计算
    A_pruned, b_pruned = remove_redundant_constraints(A, b)
    
    # 找到可行域内的初始点，确保严格在内部
    c = np.zeros(A_pruned.shape[1])
    res_feas = linprog(c, A_ub=A_pruned, b_ub=b_pruned)
    if not res_feas.success:
        raise ValueError("该约束系统无可行域")
    
    # 计算初始点到各约束的距离，确保严格内部
    distances = b_pruned - A_pruned @ res_feas.x
    min_dist = np.min(distances)
    if min_dist < 1e-8:
        # 向可行域中心微调，避免在边界上
        center = res_feas.x + 1e-6 * np.mean(A_pruned, axis=0)
    else:
        center = res_feas.x
    
    # 转换为scipy半空间要求的格式：Ax + b <= 0（原约束A@x <=b 转为A@x -b <=0）
    halfspaces = np.hstack((A_pruned, -b_pruned[:, None]))
    hs = HalfspaceIntersection(halfspaces, center)
    
    # 获取凸包顶点
    hull = ConvexHull(hs.intersections)
    vertices = hull.points
    
    # 生成随机凸组合点（权重服从Dirichlet分布，保证和为1且非负）
    feasible_points = []
    for _ in range(n_points):
        weights = np.random.dirichlet(np.ones(len(vertices)))
        point = weights @ vertices
        feasible_points.append(point)
    return np.array(feasible_points)

# 测试
points = generate_convex_hull_points(A_test, b_test, 5)
print(points)

原代码报错原因 #

1. 冗余约束过多：输入中存在大量重复的线性约束，导致半空间交集计算时出现数值不稳定；
2. 初始点在边界上：linprog得到的初始点刚好处于约束边界，而HalfspaceIntersection要求初始点必须严格在可行域内部，这是触发报错的直接原因。

内容的提问来源于stack exchange，提问作者Simon K.