1

我想解决一个不等式系统 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])
4

2 回答 2

1

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

更新

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

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

找到该集合的顶点（例如，绘制它们）称为顶点枚举问题。幸运的是，有强大的算法可以在n维中操作凸包、计算半空间交点（以及做许多其他奇妙的事情）。一个示例实现是Qhull 库。

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

在下面的：

  1. 我们找到了一个合适的可行点，或者interior_point，需要HalfspaceIntersection。
  2. 为了避免在凸集打开时出现警告（以及结果），我们使用定义边界框（由调用者提供，也用作绘制边界）的约束来扩充原始Inf系统。nanAx <= 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]);

在此处输入图像描述

于 2020-12-17T16:36:25.113 回答
0

有一个优秀的库 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 上的示例）。函数 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)
于 2020-12-18T08:04:06.257 回答