<p>似乎<a href="https://fillplots.readthedocs.io/en/latest/examples.html#two-regions-two-py" rel="nofollow noreferrer">fillplots</a>是您所需内容的超集。这应该很容易处理线性不等式</p>
<p><strong>更新</strong></p>
<p>我又在考虑这个问题,我想我会尝试看看没有<code>fillplots</code>可以做些什么,只使用标准库,比如<code>scipy</code>和<code>numpy</code></p>
<p>在这样一个不等式系统中,每个方程定义了一个半空间。该系统是所有这些半空间的交集,是一个凸集</p>
<p>查找该集合的顶点(例如,绘制它们)称为<a href="https://en.wikipedia.org/wiki/Vertex_enumeration_problem" rel="nofollow noreferrer">Vertex enumeration problem</a>。幸运的是,在<em>n</em>维度上,有强大的算法来处理凸包、计算半空间交点(以及做许多其他奇妙的事情)。一个示例实现是<a href="http://www.qhull.org/html/index.htm" rel="nofollow noreferrer">Qhull library</a></p>
<p>更幸运的是,我们可以通过<code>scipy.spacial</code>直接访问该库的各个方面,特别是:<a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.HalfspaceIntersection.html" rel="nofollow noreferrer">^{<cd5>}</a>和<a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.ConvexHull.html" rel="nofollow noreferrer">^{<cd6>}</a></p>
<p>在以下方面:</p>
<ol>
<li>我们找到一个合适的可行点,或<code>interior_point</code>,它是<code>HalfspaceIntersection</code>所需要的</李>
<li>为了避免凸集打开时出现警告(以及结果中的<code>Inf</code>、<code>nan</code>),我们使用定义边界框(由调用方提供,也用作打印边界)的约束来扩充原始系统<code>Ax <= b</code></李>
<li>我们得到半空间交点,并将它们重新排序为凸壳(这有点浪费,但我没有完全遵循<code>HalfspaceIntersection</code>返回的顺序,在2D中,壳的顶点保证为逆时针顺序)</李>
<li>我们绘制凸壳(红色),以及与方程对应的所有线</李>
</ol>
<p>我们开始:</p>
<pre><code>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
</code></pre>
<p><strong>测试</strong></p>
<p>(您的原始系统):</p>
<pre><code>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]);
</code></pre>
<p><a href="https://i.stack.imgur.com/F2IPx.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/F2IPx.png" alt="enter image description here"/></a></p>
<p>产生开放集的修改系统:</p>
<pre><code>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]);
</code></pre>
<p><a href="https://i.stack.imgur.com/sMqXT.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/sMqXT.png" alt="enter image description here"/></a></p>