如何用Sympy解决不等式？

Question

目标：我想使用Sympy实现一个傅里叶-Motzkin消除。这方面的第一步是解决一些不平等问题。

问题:使用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的表达式。我想用傅里叶-Motzkin消去所有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)]

我不明白如何处理这些ConditionSets。他们肯定不是我问题的解决办法。

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