Solve linear Inequalities

I want to solve a system of inequalities A x <= b, namely to visualize the set of solutions of this system. Are there any ways to do this in Python? The solutions I've found using the scipy library give only one vertex.

A = np.array([[-1, 1],
          [0, 1],
          [0.5, 1],
          [1.5, 1],
          [-1, 0],
          [0, -1]])
 b = np.array([1, 2, 3, 6, 0, 0])

It seems that fillplots is a superset of what you need. That should handle linear inequations very easily.

Update

I was thinking about this again, and I thought I would try to see what can be done without fillplots, just using standard libraries such as scipy and numpy.

In such a system of inequalities, each equation defines a half-space. The system is the intersection of all these half-spaces and is a convex set.

Finding the vertices of that set (for example, to plot them) is called the Vertex enumeration problem. Fortunately, there are powerful algorithms to manipulate convex hulls, compute half-space intersections (and do many other wonderful things) in n dimensions. An example implementation is the Qhull library.

Even more fortunate for us, we can access aspects of that library directly via scipy.spacial, specifically: HalfspaceIntersection and ConvexHull.

In the following:

1. We find a suitable feasible point, or interior_point, needed by HalfspaceIntersection.
2. In order to avoid warnings (and Inf, nan in the result) when the convex set is open, we augment the original system Ax <= b with constraints that define a bounding box (to be supplied by the caller, and also used as plotting boundaries).
3. We get the half-space intersections, and reorder them into a convex hull (a bit wasteful, but I didn't quite follow the order returned by HalfspaceIntersection, and in 2D the vertices of the hull are guaranteed to be in counterclockwise order).
4. We plot the convex hull (in red), and all the lines corresponding to the equations.

Here we go:

import matplotlib.pyplot as plt

import numpy as np
from scipy.spatial import HalfspaceIntersection, ConvexHull
from scipy.optimize import linprog

def feasible_point(A, b):
    # finds the center of the largest sphere fitting in the convex hull
    norm_vector = np.linalg.norm(A, axis=1)
    A_ = np.hstack((A, norm_vector[:, None]))
    b_ = b[:, None]
    c = np.zeros((A.shape[1] + 1,))
    c[-1] = -1
    res = linprog(c, A_ub=A_, b_ub=b[:, None], bounds=(None, None))
    return res.x[:-1]

def hs_intersection(A, b):
    interior_point = feasible_point(A, b)
    halfspaces = np.hstack((A, -b[:, None]))
    hs = HalfspaceIntersection(halfspaces, interior_point)
    return hs

def plt_halfspace(a, b, bbox, ax):
    if a[1] == 0:
        ax.axvline(b / a[0])
    else:
        x = np.linspace(bbox[0][0], bbox[0][1], 100)
        ax.plot(x, (b - a[0]*x) / a[1])

def add_bbox(A, b, xrange, yrange):
    A = np.vstack((A, [
        [-1,  0],
        [ 1,  0],
        [ 0, -1],
        [ 0,  1],
    ]))
    b = np.hstack((b, [-xrange[0], xrange[1], -yrange[0], yrange[1]]))
    return A, b

def solve_convex_set(A, b, bbox, ax=None):
    A_, b_ = add_bbox(A, b, *bbox)
    interior_point = feasible_point(A_, b_)
    hs = hs_intersection(A_, b_)
    points = hs.intersections
    hull = ConvexHull(points)
    return points[hull.vertices], interior_point, hs

def plot_convex_set(A, b, bbox, ax=None):
    # solve and plot just the convex set (no lines for the inequations)
    points, interior_point, hs = solve_convex_set(A, b, bbox, ax=ax)
    if ax is None:
        _, ax = plt.subplots()
    ax.set_aspect('equal')
    ax.set_xlim(bbox[0])
    ax.set_ylim(bbox[1])
    ax.fill(points[:, 0], points[:, 1], 'r')
    return points, interior_point, hs

def plot_inequalities(A, b, bbox, ax=None):
    # solve and plot the convex set,
    # the inequation lines, and
    # the interior point that was used for the halfspace intersections
    points, interior_point, hs = plot_convex_set(A, b, bbox, ax=ax)
    ax.plot(*interior_point, 'o')
    for a_k, b_k in zip(A, b):
        plt_halfspace(a_k, b_k, bbox, ax)
    return points, interior_point, hs

Tests

(Your original system):

plt.rcParams['figure.figsize'] = (6, 3)

A = np.array([[-1, 1],
          [0, 1],
          [0.5, 1],
          [1.5, 1],
          [-1, 0],
          [0, -1]])
b = np.array([1, 2, 3, 6, 0, 0])

bbox = [(-1, 5), (-1, 4)]
fig, ax = plt.subplots(ncols=2)
plot_convex_set(A, b, bbox, ax=ax[0])
plot_inequalities(A, b, bbox, ax=ax[1]);

A modified system that results in an open set:

A = np.array([
    [-1, 1],
    [0, 1],
    [-1, 0],
    [0, -1],
])
b = np.array([1, 2, 0, 0])

fig, ax = plt.subplots(ncols=2)
plot_convex_set(A, b, bbox, ax=ax[0])
plot_inequalities(A, b, bbox, ax=ax[1]);