Getting any feasible point from a system of linear inequalities in Python

Question

I have a problem with solving a system of linear inequalities in Python. What I need is basically draw n random feasible points from the inequality system. You could also call it a constraint system over an infinite domain.

So far, I tried two approaches which both not seem to work. First, I tried solving the system with scipy.optimize.linprog, having the two arrays A_lb, and b_lb. This, of course, only gives me one solution. I read online that adding noise to the system changes the solution (e.g. adding variables without relations to the initial system), as the objective function is still 0. However, it did not work. I always get the same solution.

The next Idea I had was creating the halfspace intersection, getting the convex hull and creating convex combinations of the generator points. But I did not manage to get this running in Python.

Maybe someone has an idea or knows how to implement the issue,

Thank you all very much!

Edit:

The code I am trying to use for halfspaceintersection is the following: (related to Solve linear Inequalities)

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

A_test: [[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: [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]
A_test = np.array(A_test)
b_test = np.array(b_test)

# Get a feasible point 
norm_vector = np.linalg.norm(A_test, axis=1)
A_ = np.hstack((A_test, norm_vector[:, None]))
b_ = b_test[:, None]
c = np.zeros((A_test.shape[1] + 1,))
c[-1] = -1
res = linprog(c, A_ub=A_, b_ub=b_test[:, None], bounds=(0, 1))
print(res["x"])

# Returns [2.10000000e-01 2.10000000e-01 1.80000000e-01 1.70000000e-01
# 2.30000000e-01 1.92880526e-13]

stacked_inequalities = np.hstack((A_test, b_test[:, None]))
x = HalfspaceIntersection(np.array(stacked_inequalities), res["x"][:-1])

The error I am getting is the following:

HalfspaceIntersectionError

But the underlying problem is not to get this specific program running but just find a way to randomly draw from a system of inequalities. I thought that would be the easiest approach (As the system is convex).