如何使用 Sympy 解决不平等问题？
How do I solve inequalities with Sympy?
Objective: 我想使用 Sympy 实现 Fourier-Motzkin Elimination。第一步是解决一些不平等问题。
问题：使用 Sympy 解决不等式的工具，如 solveset、solve_poly_inequality 或 reduce_inequalities 似乎不起作用。
数据：
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]
这些都是>=0的表达式。
我想用 Fourier-Motzkin Elimination 去除所有 y 变量。因此，第一步我想从 y1.
开始
所需的解决方案：
例如 list_of_inequalites[8] 即 50*y1 - y4 我应该得到 y1>=y4/50 或类似的。
最后我想要两个列表。一种表达式小于 y1，它将包含 y4/50，另一种表达式大于 y1。
我将需要这些列表用于 Fourier-Motzkin 消除的下一步。
我的尝试：
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))
这样我得到一个这样的列表：
[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)]
我不明白如何处理这些条件集。它们肯定不是我问题的解决方案。
另一种方法是使用 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),'>='))
这样我得到一个
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]
导致这个错误的不等式是50*u1 - x1 - 50*y1 + y4 >= 0。
我找到的最后一个解决不等式的方法是 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]))
但是，这次我得到以下错误：
---------------------------------------------------------------------------
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.
你有什么办法可以解决这个问题吗？
我不完全了解您的特定问题的性质。但也许我们可以解决它。
第一个解决方案
solve 也可以解决不等式问题，尽管提取正确的解决方案可能并不容易：
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]
第二种解决方案
由于你的不等式似乎是线性的，也许我们可以将它们作为方程求解，从而节省一些 post-processing。例如：
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]
第三种解决方案
我认为 solveset 返回了 ConditionSet，因为它不知道您符号的性质。如果您的符号代表真实变量，您可以对它们设置假设：
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]
Objective: 我想使用 Sympy 实现 Fourier-Motzkin Elimination。第一步是解决一些不平等问题。
问题：使用 Sympy 解决不等式的工具，如 solveset、solve_poly_inequality 或 reduce_inequalities 似乎不起作用。
数据：
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]
这些都是>=0的表达式。
我想用 Fourier-Motzkin Elimination 去除所有 y 变量。因此，第一步我想从 y1.
所需的解决方案：
例如 list_of_inequalites[8] 即 50*y1 - y4 我应该得到 y1>=y4/50 或类似的。
最后我想要两个列表。一种表达式小于 y1，它将包含 y4/50，另一种表达式大于 y1。
我将需要这些列表用于 Fourier-Motzkin 消除的下一步。
我的尝试：
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))
这样我得到一个这样的列表：
[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)]
我不明白如何处理这些条件集。它们肯定不是我问题的解决方案。
另一种方法是使用 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),'>='))
这样我得到一个
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]
导致这个错误的不等式是50*u1 - x1 - 50*y1 + y4 >= 0。
我找到的最后一个解决不等式的方法是 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]))
但是，这次我得到以下错误：
---------------------------------------------------------------------------
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.
你有什么办法可以解决这个问题吗？
我不完全了解您的特定问题的性质。但也许我们可以解决它。
第一个解决方案
solve 也可以解决不等式问题，尽管提取正确的解决方案可能并不容易：
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]
第二种解决方案
由于你的不等式似乎是线性的，也许我们可以将它们作为方程求解，从而节省一些 post-processing。例如：
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]
第三种解决方案
我认为 solveset 返回了 ConditionSet，因为它不知道您符号的性质。如果您的符号代表真实变量，您可以对它们设置假设：
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]