Search code examples
pythonsympy

Get the matrix-vector representation of a set of symbolic linear inequalities

With Sympy, it is possible to define symbolic linear inequalities, like:

enter image description here

This system is equivalent to:

enter image description here

And then it can be represented by a matrix A of coefficients and a vector b of upper bounds:

A = np.array([
  [-1, 0, 0], # -x
  [ 1, 0, 0], # x
  [ 0,-1, 0], # -y
  [ 1, 1, 0], # x+y
  [ 0, 0,-1], # -z
  [ 1, 1, 1]  # x+y+z
])
b = np.array([5, 4, 5, 3, 10, 6])

I'm wondering whether it is possible to directly get A and b from the symbolic linear inequalities, without manual gymnastics?

Solution

  • As far as I know, SymPy doesn't have a function that does what you want. But we can code one. A very basic attempt is the following, which has a lot of assumptions for it to work properly.

    from sympy import *
import numpy as np
var("x:z")
i1 = (x >= -5) & (x <= 4)
i2 = (y >= -5) & (y <= 3 - x)
i3 = (z >= -10) & (z <= 6 - 2*x - y)

# relational types: StrictGreaterThan, StrictLessThan, GreaterThan, LessThan

def func(inequalities, symbols, required_type):
    # assumptions:
    # 1. all inequalities are written with the same relational: < or <= or > or >=
    # 2. For each inequality, LHS and RHS are linear in `symbols`
    
    def get_ineq_with_correct_type(i, required_type):
        if type(i) != required_type:
            i = i.reversed
        return i
    
    # extract all inequalities, process them so they are all of the same type and
    # terms containing `symbols` are on the LHS.
    ineq = []
    for i in inequalities:
        if isinstance(i, And):
            ineq.extend([get_ineq_with_correct_type(a, required_type) for a in i.args])
        else:
            ineq.append(get_ineq_with_correct_type(i, required_type))
    
    # at this point, all inequalities should be of the same type.
    # rewrite them as expressions: LHS - RHS
    equations = [i.lhs - i.rhs for i in ineq]
    return linear_eq_to_matrix(equations, symbols)

# rewrite all inequalities as LessThan and extract A, b
A, b = func([i1, i2, i3], [x, y, z], LessThan)

print("SymPy A:", A)
print("SymPy b:", b)
A = np.array(A, dtype=float)
b = np.array(b, dtype=float)
print("Numpy A:\n", A)
print("Numpy b:\n", b)

# SymPy A: Matrix([[-1, 0, 0], [1, 0, 0], [0, -1, 0], [1, 1, 0], [0, 0, -1], [2, 1, 1]])
# SymPy b: Matrix([[5], [4], [5], [3], [10], [6]])
# Numpy A:
#  [[-1.  0.  0.]
#  [ 1.  0.  0.]
#  [ 0. -1.  0.]
#  [ 1.  1.  0.]
#  [ 0.  0. -1.]
#  [ 2.  1.  1.]]
# Numpy b:
#  [[ 5.]
#  [ 4.]
#  [ 5.]
#  [ 3.]
#  [10.]
#  [ 6.]]

    Note that this has only been tested on your example, where I only changed one coefficient in the last inequality just for fun.