Plotting Region defined by Inequality Constraints in Python (Part 2)

Author: Tran Thu Le
Date: 5/2/2022

In the previous post, ChatGPT and me introduced a method for plotting the 2D convex region defined by a set of linear inequalities $Ax\leq b$. While the method is simple, it had some limitations: the image generated is of low resolution, and the vertices of the polytope are not known precisely. In this post, we provide new methods for solving these limitations.

The following figures illustrate the difference between the new technique (left) and the previous one.

First new method (recommend)

In this method, we use pypoman to find the vertices of the convex set and then use ConvexHull in scipy to rearange these vertices. The code can be downloaded here.

You can install pypoman and scipy from the terminal as follows:

pip install scipy pypoman

Here is a minimal demo.

import numpy as  np
import  matplotlib.pyplot  as  plt
import pypoman
from scipy.spatial import ConvexHull

A = np.array([[1, 2], [1, -1], [-1, 1], [-1, -1]])
b = np.array([2, 2, 2, 1])
vertices = pypoman.compute_polytope_vertices(A, b)
vertices = np.vstack(vertices) # shape (n, 2)
indices = ConvexHull(vertices).vertices
vertices = vertices[indices, :]

plt.fill(vertices[:, 0], vertices[:, 1], color="green",
         alpha=0.5, edgecolor='black', linewidth=3)
plt.xlim([-0.5, 2.5])
plt.ylim([-0.5, 1.5])
plt.show()

Second new method (not recommend)

The Python code for this second method can be downloaded here.

Let $A\in \mathbb R^{n\times 2}$ and $b\in \mathbb R^n$. We assume that $S = \{x\in \mathbb R^2: Ax\leq b\}$ is nonempty, nondegenerate, and bounded. If it is unbounded, we can always restrict it to a box region $[p_1,q_1]\times[p_2, q_2]$. We further assume that the rows of $A$ (the normal vectors) are normalized to $1$.

Let $L_i$ be the $i$-th line, associated with the normal vector $A_i\in \mathbb R^2$ and the intercept $b_i\in \mathbb R$. To find the vertices of the convex set, we use the following pseudo-code:

S = empty set of vertices of the convex set
for i = 1, ..., n do:
  Find the angles αj of Aj and Ai for j ≠ i
  Let j_max be the index of max of negative angles αj
  Let j_min be the index of min of positive angles αj
  Let:
    I = L_j_max ∩ L_j_min
    J_max = L_j_max ∩ L_i
    J_min = L_i ∩ L_j_min
  If J_max, J_min, I forms a counter clock-wise triangle then:
    S = S ∪ {j_min}

Note that if the lines are perpendicular to $L_i$, then the intersection $I = L_{j_{\max}}\cap L_{j_{\min}}$ is $+\infty$, and in this case $J_{\max}, J_{\min}, \infty$ are considered counter clock-wise. The set of indices $S$ is the order of vertices of the convex set.

Note: The correctness of this idea is not proved yet. However, If you are looking for a more general problem: Finding the vertices of a convex set in $\mathbb R^d$ defined by a linear system of inequality, this problem is known as Vertex enumeration problem. In the next posts, we will discuss these methods in more detail. It is also important to distinguish our algorithms and Convex hull algorithms.

Conclusion

With this new method, we can plot the 2D convex region defined by a set of inequalities with higher precision and accuracy.