求解线性不等式

Question

我想求解一个不等式系统A x <= b，即可视化这个系统的解集。在 Python 中有什么方法可以做到这一点吗？我使用 scipy 库找到的解决方案仅提供一个顶点。

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])

Answer 1

fillplots 似乎是您所需要的超集。那应该很容易处理线性不等式。

更新

我又在考虑这个问题，我想我会尝试看看没有 fillplots 可以做什么，只使用 scipy 和 numpy 等标准库。

在这样的不等式系统中，每个方程定义了一半-space。系统是所有这些half-space的交集，是一个凸集

查找该集合的顶点（例如，绘制它们）称为 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。

对我们来说更幸运的是，我们可以直接通过 scipy.spacial 访问该库的各个方面，特别是：HalfspaceIntersection and ConvexHull.

在下面：

我们找到一个合适的可行点，或 interior_point，HalfspaceIntersection 需要。
为了避免在凸集打开时出现警告（以及结果中的 Inf、nan），我们使用定义边界的约束扩充原始系统 Ax <= b框（由调用者提供，也用作绘图边界）。
我们得到一半-space个交集，并将它们重新排序成一个凸包（有点浪费，但我没有完全按照HalfspaceIntersection返回的顺序，而在2D中船体的顶点保证按逆时针顺序排列。
我们绘制凸包（红色），以及与方程对应的所有直线。

我们开始：

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

测试

（您原来的系统）：

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 = 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]);

Answer 2

有一个优秀的库 pypoman 解决了顶点枚举问题，可以帮助解决你的问题，但不幸的是它只输出集合的顶点，而不是可视化。顶点可能是无序的，如果没有额外的操作，可视化将不正确。要解决此问题，您可以使用此站点的算法 https://habr.com/ru/post/144921/（Graham 扫描或算法 Jarvis）。

这是一个示例代码：

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)

我也在为我的硕士论文写一个 intvalpy 库（还没有文档，只有 githab 上的例子）。函数 lineqs 也可以帮助你。它求解系统 A x >= b 并输出有序顶点并将集合可视化。

对于您的问题，代码如下所示：

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)

Answer 3

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)

求解线性不等式

Solve linear Inequalities

python

numpy

matplotlib

inequalities