Блог

Я хочу решить систему неравенств 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])
 

Кажется, что fillplots — это надмножество того, что вам нужно. Это должно очень легко обрабатывать линейные неравенства.

Обновить

Я снова думал об этом, и я подумал, что попробую посмотреть, без чего можно обойтись fillplots , просто используя стандартные библиотеки, такие как scipy и numpy .

В такой системе неравенств каждое уравнение определяет полупространство. Система является пересечением всех этих полупространств и представляет собой выпуклое множество.

Поиск вершин этого множества (например, для их построения) называется задачей перечисления вершин. К счастью, существуют мощные алгоритмы для управления выпуклыми оболочками, вычисления полупространственных пересечений (и многих других замечательных вещей) в n измерениях. Примером реализации является библиотека Qhull.

Еще более удачно для нас, что мы можем получить доступ к аспектам этой библиотеки напрямую через scipy.spacial , в частности: HalfspaceIntersection и ConvexHull .

В следующем:

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

Существует отличная библиотека pypoman, которая решает проблему перечисления вершин и может помочь с вашей проблемой, но, к сожалению, она выводит только вершины множества, а не визуализацию. Вершины могут быть неупорядоченными, и без дополнительных действий визуализация будет некорректной. Чтобы преодолеть эту проблему, вы можете использовать алгоритмы с этого сайта https://habr.com/ru/post/144921 / (Сканирование Грэма или алгоритм Джарвиса).

Вот пример кода:

 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)
 

 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)