How do I solve inequalities with Sympy?

Question

Objective: I want to implement a Fourier-Motzkin Elimination using Sympy. The first step for this would be to solve a number of inequalities.

Problem: The tools for solving inequalities with Sympy like solveset, solve_poly_inequality or reduce_inequalities do no seem to work.

Data:

from sympy import *
u1, u2, x1, x2 = symbols('u1 u2 x1 x2')
y1, y2, y3, y4, y5 = symbols('y1 y2 y3 y4 y5')

list_of_inequalities = [50*u2 - y5, -x2 + y5, 0, -u2 + 1, -35*u2 + y4 + y5, 35*u2 + x1 + x2 - y4 - y5 - 35, 35*u2 - y4 - y5, -35*u2 - x1 - x2 + y4 + y5 + 35, 50*y1 - y4, 50*u1 - x1 - 50*y1 + y4, -50*y1 + y4, -50*u1 + x1 + 50*y1 - y4, u2 - y1, -u1 - u2 + y1 + 1, 50*u2 - y5, -50*u2 - x2 + y5 + 50, 65*y1, 65*u1 - 65*y1, 35*u2 + 65*y1 - y4 - y5]

These are all expressions that are >=0. I want to get rid of all the y-variables with Fourier-Motzkin Elimination. So, in a first step I would like to start with y1.

Desired Solution:

For example for list_of_inequalites[8] which is 50*y1 - y4 I should get y1>=y4/50 or similar. In the end I want to have two lists. One with expressions that are smaller than y1 which would contain y4/50 and one with expressions larger than y1. I will need these lists for the next step in the Fourier-Motzkin Elimination.

My Try:

y_1=[]
for eq in list_of_equations:
    expr= eq>=0
    if y1 in eq.free_symbols:
        y_1.append(solveset(expr.lhs>=0,y1,domain=S.Reals))

This way I get a list like this:

[ConditionSet(y1, 50*y1 - y4 >= 0, Reals), ConditionSet(y1, 50*u1 - x1 - 50*y1 + y4 >= 0, Reals), ConditionSet(y1, -50*y1 + y4 >= 0, Reals), ConditionSet(y1, -50*u1 + x1 + 50*y1 - y4 >= 0, Reals), ConditionSet(y1, u2 - y1 >= 0, Reals), ConditionSet(y1, -u1 - u2 + y1 + 1 >= 0, Reals), Interval(0, oo), ConditionSet(y1, 65*u1 - 65*y1 >= 0, Reals), ConditionSet(y1, 35*u2 + 65*y1 - y4 - y5 >= 0, Reals)]

I do not understand how to deal with these ConditionSets. They are certainly not the solution of my problem.

Another way would be to work with solve_poly_inequality:

for eq in list_of_equations:   
    expr= eq>=0
    if y1 in eq.free_symbols:
        y_1.append(solve_poly_inequality(Poly(expr.lhs,y1),'>='))

This way I get a

NotImplementedError                       Traceback (most recent call last)
<ipython-input-269-686426e9455b> in <module>
      9     expr= eq>=0
     10     if y1 in eq.free_symbols:
---> 11         y_1.append(solve_poly_inequality(Poly(expr.lhs,y1),'>='))

~\Anaconda3\lib\site-packages\sympy\solvers\inequalities.py in solve_poly_inequality(poly, rel)
     56                 "could not determine truth value of %s" % t)
     57 
---> 58     reals, intervals = poly.real_roots(multiple=False), []
     59 
     60     if rel == '==':

~\Anaconda3\lib\site-packages\sympy\polys\polytools.py in real_roots(f, multiple, radicals)
   3498 
   3499         """
-> 3500         reals = sympy.polys.rootoftools.CRootOf.real_roots(f, radicals=radicals)
   3501 
   3502         if multiple:

~\Anaconda3\lib\site-packages\sympy\polys\rootoftools.py in real_roots(cls, poly, radicals)
    385     def real_roots(cls, poly, radicals=True):
    386         """Get real roots of a polynomial. """
--> 387         return cls._get_roots("_real_roots", poly, radicals)
    388 
    389     @classmethod

~\Anaconda3\lib\site-packages\sympy\polys\rootoftools.py in _get_roots(cls, method, poly, radicals)
    717             raise PolynomialError("only univariate polynomials are allowed")
    718 
--> 719         coeff, poly = cls._preprocess_roots(poly)
    720         roots = []
    721 

~\Anaconda3\lib\site-packages\sympy\polys\rootoftools.py in _preprocess_roots(cls, poly)
    696         if not dom.is_ZZ:
    697             raise NotImplementedError(
--> 698                 "sorted roots not supported over %s" % dom)
    699 
    700         return coeff, poly

NotImplementedError: sorted roots not supported over ZZ[x1,y4,u1]

The inequality that causes this error is 50*u1 - x1 - 50*y1 + y4 >= 0.

The last way I found for solving inequalities is reduce_inequalities:

for eq in list_of_equations:   
    expr= eq>=0
    if y1 in eq.free_symbols:
        y_1.append(reduce_inequalities(expr>=0,[y1]))

But, this time I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-266-50cdf532f9fa> in <module>
     22     expr= eq>=0
     23     if y1 in eq.free_symbols:
---> 24         y_1.append(reduce_inequalities(expr>=0,[y1]))
     25 #     print(y_1[-1])
     26 

~\Anaconda3\lib\site-packages\sympy\solvers\inequalities.py in reduce_inequalities(inequalities, symbols)
    985 
    986     # solve system
--> 987     rv = _reduce_inequalities(inequalities, symbols)
    988 
    989     # restore original symbols and return

~\Anaconda3\lib\site-packages\sympy\solvers\inequalities.py in _reduce_inequalities(inequalities, symbols)
    900             if len(common) == 1:
    901                 gen = common.pop()
--> 902                 other.append(_solve_inequality(Relational(expr, 0, rel), gen))
    903                 continue
    904             else:

~\Anaconda3\lib\site-packages\sympy\core\relational.py in __new__(cls, lhs, rhs, rop, **assumptions)
     87                             Eq and Ne; all other relationals expect
     88                             real expressions.
---> 89                         '''))
     90             # \\\
     91             return rv
TypeError: 
A Boolean argument can only be used in Eq and Ne; all other
relationals expect real expressions.

Do you have any ideas how I can solve this problem?

When you ask about an error, show the **whole thing**! For example what actually produces the `truth value` error? `if y>0: print('yes')` produces that error, since `y>0` is a symbolic relation, which can't be used in a context like `if` with expects a simple True/False. — hpaulj, Feb 22 '22 at 17:53
Hmm, actually I couldn't reproduce the ````truth value```` error, but got ````NotImplementedError: sorted roots not supported over ZZ[x1,y4,u1]````. Otherwise, I now included the whole error. — R.Deilke, Feb 22 '22 at 18:03
There is some work [here](https://github.com/sympy/sympy/pull/22389) with a reference to another approach there, too. — smichr, Feb 22 '22 at 23:32
`u1, u2, x1, x2 = smp.symbols('u1 u2 x1 x2')` gives `NameError: name 'smp' is not defined` — Nasser, Feb 23 '22 at 02:28
This looks like it's using iPython and Anaconda - infrequently restarting the kernel and re-running all the cells or copying the relevant parts to a simple script may help a great amount with debugging (for example as @Nasser shows, there's probably at least a line with `import sympy as smp` not shown amongst any amount of others..) — ti7, Feb 23 '22 at 02:46
I restarted the kernel and re-ran all the cells in a simple script. But, the errors are exactly the same. I also got rid of the 'smp'. — R.Deilke, Feb 23 '22 at 10:08
Beyond some basics it is hard to generalize `sympy` solutions. You have to study the docs. I haven't done much with inequalities, or what I think is being called `relationals`. And if your `solve_set` produces `ConditionSets`, then you need to study that part of the documentation. Also the coverage of subpackages varies widely; this is a open source project, not something that was developed and maintained by a company with paid staff. — hpaulj, Feb 23 '22 at 16:35
Is there any other way to solve linear inequalities with Python? I actually solved the example problem by hand. So I could also start to work with string manipulations to solve it. I just thought Sympy should have a way to solve this easily. — R.Deilke, Feb 24 '22 at 14:40

score 2 · Accepted Answer · answered Apr 01 '22 at 18:57

I don't know exactly the nature of your particular problem. But maybe we can work around it.

First solution

solve is also capable of solving inequalities, though it might not be easy to extract the proper solution:

sols = [solve(t >= 0, y1) for t in list_of_inequalities if y1 in S(t).free_symbols]
sols

# [(y1 < oo) & (y1 >= y4/50),
#  (-oo < y1) & (y1 <= u1 - x1/50 + y4/50),
#  (-oo < y1) & (y1 <= y4/50),
#  (y1 < oo) & (y1 >= u1 - x1/50 + y4/50),
#  (y1 <= u2) & (-oo < y1),
#  (y1 < oo) & (y1 >= u1 + u2 - 1),
#  (0 <= y1) & (y1 < oo),
#  (y1 <= u1) & (-oo < y1),
#  (y1 < oo) & (y1 >= -7*u2/13 + y4/65 + y5/65)]

# now a bit of post-processing:
from sympy.core.numbers import Infinity, NegativeInfinity
test_inf = lambda x: x.has(Infinity) or x.has(NegativeInfinity)
correct_arg = lambda x, y: x.args[0] if not x.args[0].has(y) else x.args[1]
final = [
    correct_arg(u, y1) for u in # extract the arguments from Relational that doesn't contain y1
    [t[0] if not test_inf(t[0]) else t[1] for t in # exclude arguments containing infinity
    [s.args[:2] for s in sols]] # extract args from Boolean And
]
final

# [y4/50,
#  u1 - x1/50 + y4/50,
#  y4/50,
#  u1 - x1/50 + y4/50,
#  u2,
#  u1 + u2 - 1,
#  0,
#  u1,
#  -7*u2/13 + y4/65 + y5/65]

Second solution

Since your inequalities appear to be linear, maybe we can solve them as equations, thus sparing us some post-processing. For example:

sols = [solve(t, y1)[0] for t in list_of_inequalities if y1 in S(t).free_symbols]
sols

# [y4/50,
#  u1 - x1/50 + y4/50,
#  y4/50,
#  u1 - x1/50 + y4/50,
#  u2,
#  u1 + u2 - 1,
#  0,
#  u1,
#  -7*u2/13 + y4/65 + y5/65]

Third solution

I believe that solveset returned ConditionSet because it doesn't know the nature of your symbols. If your symbols represent real variables, you can set assumptions on them:

u1, u2, x1, x2 = symbols('u1 u2 x1 x2', real=True)
y1, y2, y3, y4, y5 = symbols('y1 y2 y3 y4 y5', real=True)

list_of_inequalities = [50*u2 - y5, -x2 + y5, 0, -u2 + 1, -35*u2 + y4 + y5, 35*u2 + x1 + x2 - y4 - y5 - 35, 35*u2 - y4 - y5, -35*u2 - x1 - x2 + y4 + y5 + 35, 50*y1 - y4, 50*u1 - x1 - 50*y1 + y4, -50*y1 + y4, -50*u1 + x1 + 50*y1 - y4, u2 - y1, -u1 - u2 + y1 + 1, 50*u2 - y5, -50*u2 - x2 + y5 + 50, 65*y1, 65*u1 - 65*y1, 35*u2 + 65*y1 - y4 - y5]

sols = [solveset(t >= 0, y1, S.Reals) for t in list_of_inequalities if y1 in S(t).free_symbols]
sols

# [Interval(y4/50, oo),
#  Interval(-oo, u1 - x1/50 + y4/50),
#  Interval(-oo, y4/50),
#  Interval(u1 - x1/50 + y4/50, oo),
#  Interval(-oo, u2),
#  Interval(u1 + u2 - 1, oo),
#  Interval(0, oo),
#  Interval(-oo, u1),
#  Interval(-7*u2/13 + y4/65 + y5/65, oo)]

# extract the expressions, disregarding the infinity symbols
from sympy.core.numbers import Infinity, NegativeInfinity
test_inf = lambda x: x.has(Infinity) or x.has(NegativeInfinity)
[t[0] if not test_inf(t[0]) else t[1] for t in [s.args[:2] for s in sols]]
# [y4/50,
#  u1 - x1/50 + y4/50,
#  y4/50,
#  u1 - x1/50 + y4/50,
#  u2,
#  u1 + u2 - 1,
#  0,
#  u1,
#  -7*u2/13 + y4/65 + y5/65]

I have solved it differently, but your answer is much more flexible. Thanks! — R.Deilke, Apr 20 '22 at 12:35

How do I solve inequalities with Sympy?

1 Answers1

First solution

Second solution

Third solution