python - 如何使用 Sympy 解决不平等问题？

Question

目标:我想使用 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。我将需要这些列表用于傅立叶-莫茨金消除的下一步。

我的尝试:

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)]

我不明白如何处理这些 ConditionSet。它们当然不是我的问题的解决方案。

另一种方法是使用 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.

你有什么办法可以解决这个问题吗？

Answer 1

我不完全了解您的特定问题的性质。但也许我们可以解决这个问题。

第一个解决方案

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]

第二种方案

由于您的不等式看起来是线性的，也许我们可以将它们作为方程求解，从而节省一些后处理。例如:

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]