Как исправить текущий решатель неравенства и есть ли лучший подход/библиотека?

Я не уверен, как лучше это решить, но с помощью кода, который уже доступен здесь, вы можете сделать следующее:

>>> from symyp import *
>>> from sympy.abc import x,y,z
>>> eqs = [31*x + 2.9*y + 25*z <= 160, 23*y + 1.3*z <= 160, 3.6*x + 0.9*y + 33*z <= 50, 165*x + 110.9*y + 402*z <= 1730]
>>> eqs = [nsimplify(i, rational=True) for i in eqs]
>>> sol = solve_linear_inequalities(eqs,(x,y,z))

Вы получите решение для y, указывающее, что оно находится в Interval(-oo, min(...functions of x and z...)). Выполнение того же для функций x и z в min дает решение для z, подобное Interval(-oo, min(...const or functions of x...)). Наконец, то же самое для функций x в min показывает, что x находится в интервале (-oo, 5.16). Таким образом, для x в диапазоне (1,6) вы можете вычислить соответствующие диапазоны для z, а затем y. Вот просто цикл грубой силы:

>>> yy = sol[y]
for xx in range(6):
    for zz in range(6):
        yi = yy.subs(x,xx).subs(z,zz).intersection(Range(1,oo))
        if not yi:continue
        print(xx,yi,zz)

который дает

0 Range(1, 7, 1) 0
0 Range(1, 7, 1) 1
1 Range(1, 7, 1) 0
1 Range(1, 7, 1) 1
2 Range(1, 7, 1) 0
2 Range(1, 7, 1) 1
3 Range(1, 7, 1) 0
3 Range(1, 7, 1) 1
4 Range(1, 7, 1) 0
4 Range(1, 3, 1) 1
5 Range(1, 2, 1) 0

Без целевой функции, которую нужно максимизировать, похоже, что существует конечное число целочисленных решений ваших уравнений. Возможно, этого достаточно, чтобы подтолкнуть вас к более общему решению вашей проблемы? Кроме того, если вы укажете диапазон интересов для x, y и z, это поможет, например. если все они положительные, то дано решение для каждой переменной

>>> sol = solve_linear_inequalities(eqs + [x>0,y>0,z>0], (x, y, z))
>>> for i in sol:
...     print(i, sol[i])
z Interval(0, 50/33)
x Interval(0, Min(160/31 - 25*z/31, 346/33 - 134*z/55, 125/9 - 55*z/6))
y Interval(0, Min(160/23 - 13*z/230, -310*x/29 - 250*z/29 + 1600/29, -4*x - 110*z/3 + 500/9, -1650*x/1109 - 4020*z/1109 + 17300/1109))

Поэтому выберите значение из диапазона для z, рассчитайте допустимые диапазоны для x и выберите значение x и, наконец, вычислите допустимый диапазон для y.