解线性不等式

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

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

查找该集合的顶点（例如，绘制它们）称为Vertex enumeration problem。幸运的是，在n维度上，有强大的算法来处理凸包、计算半空间交点（以及做许多其他奇妙的事情）。一个示例实现是Qhull library

更幸运的是，我们可以通过scipy.spacial直接访问该库的各个方面，特别是：^{}和^{}

在以下方面：

1. 我们找到一个合适的可行点，或interior_point，它是HalfspaceIntersection所需要的
2. 为了避免凸集打开时出现警告（以及结果中的Inf、nan），我们使用定义边界框（由调用方提供，也用作打印边界）的约束来扩充原始系统Ax <= b
3. 我们得到半空间交点，并将它们重新排序为凸壳（这有点浪费，但我没有完全遵循HalfspaceIntersection返回的顺序，在2D中，壳的顶点保证为逆时针顺序）
4. 我们绘制凸壳（红色），以及与方程对应的所有线

我们开始：

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 scan或algo 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上的示例）。功能行也可以帮助您。它解决了系统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)