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:
interior_point, needed by HalfspaceIntersection.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).HalfspaceIntersection, and in 2D the vertices of the hull are guaranteed to be in counterclockwise order).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]);
There is an excellent library pypoman which solves the vertex enumerate problem and can help with your problem, but unfortunately it only outputs the vertices of the set, not the visualization. The vertices may be disordered and without additional actions the visualization will not be correct. To overcome this problem, you can use the algorithms from this site https://habr.com/ru/post/144921/ (Graham scan or algo Jarvis).
Here is a sample code:
import pypoman
import cdd
import matplotlib.pyplot as plt
def grahamscan(A):
def rotate(A,B,C):
return (B[0]-A[0])*(C[1]-B[1])-(B[1]-A[1])*(C[0]-B[0])
n = len(A)
if len(A) == 0:
return A
P = np.arange(n)
for i in range(1,n):
if A[P[i]][0]<A[P[0]][0]:
P[i], P[0] = P[0], P[i]
for i in range(2,n):
j = i
while j>1 and (rotate(A[P[0]],A[P[j-1]],A[P[j]])<0):
P[j], P[j-1] = P[j-1], P[j]
j -= 1
S = [P[0],P[1]]
for i in range(2,n):
while rotate(A[S[-2]],A[S[-1]],A[P[i]])<0:
del S[-1]
S.append(P[i])
return S
def compute_poly_vertices(A, b):
b = b.reshape((b.shape[0], 1))
mat = cdd.Matrix(np.hstack([b, -A]), number_type='float')
mat.rep_type = cdd.RepType.INEQUALITY
P = cdd.Polyhedron(mat)
g = P.get_generators()
V = np.array(g)
vertices = []
for i in range(V.shape[0]):
if V[i, 0] != 1: continue
if i not in g.lin_set:
vertices.append(V[i, 1:])
return vertices
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])
vertices = np.array(compute_poly_vertices(A, b))
print(vertices)
vertices = np.array(vertices[grahamscan(vertices)])
x, y = vertices[:, 0], vertices[:, 1]
fig=plt.figure(figsize=(15,15))
ax = fig.add_subplot(111, title="Solution")
ax.fill(x, y, linestyle = '-', linewidth = 1, color='gray', alpha=0.5)
ax.scatter(x, y, s=10, color='black', alpha=1)
I am also writing an intvalpy library for my master's thesis (no documentation yet, only examples on githab). The function lineqs can also help you. It solves the system A x >= b and outputs ordered vertices and visualizes the set.
For your problem, the code looks like this:
from intvalpy import lineqs
import numpy as np
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])
lineqs(-A, -b)
import numpy as np
import cdd as pcdd
from fractions import Fraction
A = np.array(
[[-1, 1],
[0, 1],
[Fraction(1,2), 1],
[Fraction(3,2), 1],
[-1, 0],
[0, -1]]
)
b = np.array([[1], [2], [3], [6], [0], [0]])
M = np.hstack( (b, -A) )
mat = pcdd.Matrix(M, linear=False, number_type="fraction")
mat.rep_type = pcdd.RepType.INEQUALITY
poly = pcdd.Polyhedron(mat)
ext = poly.get_generators()
print(ext)
