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?
I don't know exactly the nature of your particular problem. But maybe we can work around it.
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]
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]
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]