ConstrainedZonotope¶

This python notebook goes over various functionalities of pycvxset for manipulating and plotting constrained zonotopes.

License information: pycvxset code is released under AGPL-3.0-or-later license, as found in the LICENSE.md file. The documentation for pycvxset is released under CC-BY-4.0 license, as found in the LICENSES/CC-BY-4.0.md.

Constrained zonotope characterizes a convex and compact set in the working dimension $\mathbb{R}^n$, and is an affine transformation of a constrained zonotope defined in latent space $B_\infty(A_e, b_e)\subset \mathbb{R}^{N_C}$. Here, $B_\infty(A_e, b_e)$ is defined as the intersection of a unit $\ell_\infty$-norm ball and a collection of $M_C$ linear constraints. Formally, $$\mathcal{C} = \{G \xi + c | \xi \in B_\infty(A_e, b_e)\} \subset \mathbb{R}^n,$$ where $$B_\infty(A_e, b_e) = \{\xi | \| \xi \|_\infty \leq 1, A_e \xi = b_e\} \subset \mathbb{R}^{N_C}.$$

Note that every constrained zonotope is a polytope, and vice-versa. Consequently, pycvxset uses Python's duck typing to treat Polytope and ConstrainedZonotope objects (almost) equally.

We organize this notebook on pycvxset operations on ConstrainedZonotopes as follows:

The links below are for the HTML page in the documentation website, and they will not work in the Jupyter notebook.

  1. Representations and plotting
    1. Constrained zonotope from polytopes
      1. Removing redundancies
    2. Using polytopes
      1. Plotting a constrained zonotope
      2. Generating direction vectors for 3D plotting
    3. Axis-aligned cuboids
    4. Empty constrained zonotopes (and checking for emptiness)
    5. Embedded constrained zonotopes
  2. Operations involving a single constrained zonotope
    1. Affine transformation
      1. Negation and scaling by a scalar
    2. Inverse affine transformation with an invertible linear map
    3. Projection of a constrained zonotope
    4. Project a point on to a constrained zonotope
    5. Containment: Check if a point is in a constrained zonotope
    6. Support function computation
    7. Chebyshev centering
    8. Maximum volume inscribing ellipsoid
    9. Minimum volume circumscribing rectangle
    10. Interior point computation
    11. Volume computation
    12. Cartesian product with itself
  3. Operations involving two constrained zonotopes or a constrained zonotope and another set
    1. Minkowski sum with another constrained zonotope
    2. Pontryagin difference with a zonotope or an ellipsoid
    3. Intersection with another constrained zonotope, polytope, a halfspace, or an affine set
      1. Intersection under inverse affine map
      2. Intersection with a halfspace
      3. Intersection with an affine set
      4. Slicing
    4. Containment: Check if a constrained zonotope is in another constrained zonotope or polytope

Representations and plotting¶

In preparation to run the notebook, we import the necessary packages.

In [1]:
from pycvxset import (
    Ellipsoid,
    Polytope,
    ConstrainedZonotope,
    approximate_volume_from_grid,
    spread_points_on_a_unit_sphere,
)
import numpy as np
import matplotlib.pyplot as plt

We also set up matplotlib to have a standard figure size, and call matplotlib widget to enable 3D plotting in the notebook.

In [2]:
plt.rcParams["figure.figsize"] = [5, 3]
plt.rcParams["figure.dpi"] = 150
%matplotlib widget

Constrained zonotopes and zonotopes¶

We can define constrained zonotopes (and zonotopes) directly in pycvxset.

In [3]:
C1 = ConstrainedZonotope(G=np.eye(2), c=[0, 0], Ae=[1, 0], be=[1])
print("C1 is a", repr(C1))
C2 = ConstrainedZonotope(G=np.eye(2), c=[0, 0])
print("C2 is a", repr(C2))

C1 is a Constrained Zonotope in R^2
	with latent dimension 2 and 1 equality constraint
C2 is a Constrained Zonotope in R^2
	that is a zonotope with latent dimension 2

Removing redundancies¶

pycvxset can detect and remove redundancies in representation.

In [4]:
C1a = ConstrainedZonotope(
    G=np.hstack((np.eye(2), 2 * np.eye(2))),
    c=[0, 0],
    Ae=np.eye(4),
    be=[1, 1, 1, 1],
)
print("C1a is a", repr(C1a))
print("C1a.G.shape", C1a.G.shape)
print("C1a.Ae.shape", C1a.Ae.shape)
C1a.remove_redundancies()
print("C1a is a", repr(C1a))
print("C1a.G.shape", C1a.G.shape)
print("C1a.Ae.shape", C1a.Ae.shape)

C1a is a Constrained Zonotope in R^2
	with latent dimension 4 and 4 equality constraints
C1a.G.shape (2, 4)
C1a.Ae.shape (4, 4)
C1a is a Constrained Zonotope in R^4
	that is a zonotope representing a single point
C1a.G.shape (4, 0)
C1a.Ae.shape (0, 0)

Using polytopes¶

We create a $\mathbb{R}^3$-dimensional constrained zonotope C_3D from the corresponding polytope P_hrep. We confirm that the constrained zonotope is three-dimensional. You could also define a polytope in V-rep and construct a constrained zonotope.

In [5]:
A = np.array(
    [[1, 0, 0], [0, 1, 0], [0, 0, 1], [-1, 0, 0], [0, -1, 0], [0, 0, -1]]
)
b = [2, 3, 1, 2, 3, 1]
P_3D = Polytope(A=A, b=b)
C_3D = ConstrainedZonotope(polytope=P_3D)
print(repr(C_3D))
print(f"The dimension of C_3D is {C_3D.dim:d}.")
print("G", np.array2string(C_3D.G, precision=2, suppress_small=True))
print("c", np.array2string(C_3D.c, precision=2, suppress_small=True))
print("Ae", np.array2string(C_3D.Ae, precision=2, suppress_small=True))
print("be", np.array2string(C_3D.be, precision=2, suppress_small=True))
print("He", np.array2string(C_3D.He, precision=2, suppress_small=True))
print("is_bounded?", C_3D.is_bounded)
print("is_full_dimensional?", C_3D.is_full_dimensional)

Constrained Zonotope in R^3
	with latent dimension 9 and 6 equality constraints
The dimension of C_3D is 3.
G [[2. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 3. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0.]]
c [-0. -0.  0.]
Ae [[ 2.  0.  0. -2.  0.  0.  0.  0.  0.]
 [ 0.  3.  0.  0. -3.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0. -1.  0.  0.  0.]
 [-2.  0.  0.  0.  0.  0. -2.  0.  0.]
 [ 0. -3.  0.  0.  0.  0.  0. -3.  0.]
 [ 0.  0. -1.  0.  0.  0.  0.  0. -1.]]
be [0. 0. 0. 0. 0. 0.]
He [[ 2.  0.  0. -2.  0.  0.  0.  0.  0.  0.]
 [ 0.  3.  0.  0. -3.  0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0. -1.  0.  0.  0.  0.]
 [-2.  0.  0.  0.  0.  0. -2.  0.  0.  0.]
 [ 0. -3.  0.  0.  0.  0.  0. -3.  0.  0.]
 [ 0.  0. -1.  0.  0.  0.  0.  0. -1.  0.]]
is_bounded? True
is_full_dimensional? True

Plotting a constrained zonotope¶

Next, we plot the constrained zonotope. Note that pycvxset internally performs a polytopic approximation (by default inner, but an outer-approximation can be requested as well) in order to plot the constrained zonotope C_3D.

In [6]:
ax, _, _ = C_3D.plot(patch_args={"label": "C_3D"})
ax.legend()
ax.set_title("Plotting C_3D")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z");

Generating direction vectors for 3D plotting¶

Internally, pycvxset spreads points on a 3-dimensional unit sphere, and then performs ray-shooting. This is accomplished by using spread_points_on_a_unit_sphere of pycvxset and extreme method of ConstrainedZonotope.

The plot below shows a set of vectors spread on the unit sphere.

In [7]:
dir_vectors = spread_points_on_a_unit_sphere(3)[0]
plt.figure()
ax = plt.subplot(111, projection="3d")
ax.scatter(dir_vectors[:, 0], dir_vectors[:, 1], dir_vectors[:, 2])
ax.set_aspect("equal")

Since the computation of these points takes few seconds, we will make use of dir_vectors in all future 3D plots.

Axis-aligned cuboids¶

Next, we define an axis-aligned cuboid using the bounds. Internally, pycvxset stores the constrained zonotope in H-representation. One can also define an axis-aligned cuboid based on their center and half-sides.

In [8]:
axis_aligned_cuboid_from_bounds = ConstrainedZonotope(
    lb=[-1, 2], ub=[5, 4]
)
print(repr(axis_aligned_cuboid_from_bounds))
# Plot axis_aligned_cuboid_from_bounds
ax, _, _ = axis_aligned_cuboid_from_bounds.plot(
    patch_args={"label": "Rectangle"}, vertex_args={"visible": True}
)

axis_aligned_cuboid_from_center_and_sides = ConstrainedZonotope(
    c=[0, 0], h=1
)
print(repr(axis_aligned_cuboid_from_center_and_sides))
# Plot axis_aligned_cuboid_from_center_and_sides
axis_aligned_cuboid_from_center_and_sides.plot(
    ax=ax,
    patch_args={"label": "Cube", "facecolor": "m"},
    vertex_args={"visible": True},
)
ax.set_aspect("equal")
ax.legend(loc="best")
ax.set_title("Axis-aligned cuboids");

Constrained Zonotope in R^2
	that is a zonotope with latent dimension 2

Constrained Zonotope in R^2
	that is a zonotope with latent dimension 2

Empty constrained zonotope (and checking for emptiness)¶

We can also define an empty constrained zonotope of desired dimension as well. pycvxset can also check if a constrained zonotope is empty or nonempty.

In [9]:
empty_constrained_zonotope = ConstrainedZonotope(dim=3)
print(empty_constrained_zonotope)
print(
    "Is empty_constrained_zonotope empty?",
    empty_constrained_zonotope.is_empty,
)

Constrained Zonotope (empty) in R^3
Is empty_constrained_zonotope empty? True

We also confirm if previously defined constrained zonotopes were not empty.

In [10]:
print("Was C_3D empty?", C_3D.is_empty)
print(
    "Was axis_aligned_cuboid_from_bounds empty?",
    axis_aligned_cuboid_from_bounds.is_empty,
)
print(
    "Was axis_aligned_cuboid_from_center_and_sides empty?",
    axis_aligned_cuboid_from_center_and_sides.is_empty,
)

Was C_3D empty? False
Was axis_aligned_cuboid_from_bounds empty? False
Was axis_aligned_cuboid_from_center_and_sides empty? False

Embedded constrained zonotopes¶

Similarly to Polytope class, we can embed a lower-dimensional constrained zonotope in a higher-dimensional space.

In [11]:
unit_box_in_3D = ConstrainedZonotope(c=[0, 0, 0], h=1)
print(repr(unit_box_in_3D))
low_dim_P_in_high_dim_space = unit_box_in_3D.intersection_with_affine_set(
    Ae=[0.03, 0, 0.2], be=[0]
)
print(repr(low_dim_P_in_high_dim_space))

Constrained Zonotope in R^3
	that is a zonotope with latent dimension 3
Constrained Zonotope in R^3
	with latent dimension 3 and 1 equality constraint

However, plotting low-dimensional constrained zonotopes embedded in high-dimensional volumes is difficult since ray shooting approach may completely miss the affine set in which the constrained zonotope lives.

Operations involving a single constrained zonotope¶

We provide examples for the following operations supported by pycvxset on the ConstrainedZonotope objects.

  1. Affine transformation
  2. Inverse affine transformation with an invertible linear map
  3. Projection of a constrained zonotope
  4. Project a point on to a constrained zonotope
  5. Support function computation
  6. Chebyshev centering
  7. Maximum volume inscribing ellipsoid
  8. Minimum volume circumscribing ellipsoid
  9. Minimum volume circumscribing rectangle
  10. Containment: Check if a point is in a constrained zonotope
  11. Volume computation
  12. Interior point computation

For constrained zonotope transformations, we will use OLD_POLYTOPE_COLOR_1 to denote the old constrained zonotope and NEW_POLYTOPE_COLOR to denote the new constrained zonotope.

In [12]:
OLD_CONSTRAINED_ZONOTOPE_COLOR_1 = "skyblue"
NEW_CONSTRAINED_ZONOTOPE_COLOR = "lightgray"

Affine transformation¶

Affine transformations of constrained zonotopes can rotate, shear, and translate constrained zonotopes.

Recall that rotation matrix $R(\theta)=\left[\begin{array}{cc}\cos(\theta)&-\sin(\theta)\\\sin(\theta)&\cos(\theta)\end{array}\right]$ rotates a 2D vector by $\theta\in(-\pi,\pi]$. Thus, given a 2D constrained zonotope $\mathcal{P}$, we are interested in computing the constrained zonotope $$\mathcal{Q}=R(\theta)\mathcal{P} + \left[\begin{array}{c} 2\\-2\end{array}\right].$$

Similarly to Polytope class, the ConstrainedZonotope class overloads the operator @ and + to accomplish these transformations. Alternatively, you can use affine_map and plus methods.

In [13]:
P_from_V = Polytope(V=[[0, 0], [0, 1], [2, 0]])
CZ = ConstrainedZonotope(polytope=P_from_V)
rotation_degree = 45
rotation_radians = np.deg2rad(rotation_degree)
rotation_matrix = np.array(
    [
        [np.cos(rotation_radians), -np.sin(rotation_radians)],
        [np.sin(rotation_radians), np.cos(rotation_radians)],
    ]
)
translation = np.array([[2, -2]]).T
# Computation via operator overloading
new_constrained_zonotope_from_affine_transformation = (
    rotation_matrix @ CZ + translation
)

patch_args_old = {
    "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
    "label": "Original constrained zonotope",
}
patch_args_new = {
    "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
    "label": "Affine transformed constrained zonotope",
}
ax, _, _ = CZ.plot(patch_args=patch_args_old)
new_constrained_zonotope_from_affine_transformation.plot(
    ax=ax, patch_args=patch_args_new
)
ax.legend()
ax.set_title("Affine transformation")

# Explicit computation via methods
new_constrained_zonotope_from_affine_transformation_explicit = (
    CZ.affine_map(rotation_matrix).plus(translation)
)
print(
    "Are constrained zonotopes obtained from these methods the same?",
    new_constrained_zonotope_from_affine_transformation
    == new_constrained_zonotope_from_affine_transformation_explicit,
)

Set parameter ServerPassword

Set parameter TokenServer to value "atlas3.merl.com"

Are constrained zonotopes obtained from these methods the same? True
In [14]:
patch_args_old = {
    "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
    "label": "Original constrained zonotope",
}
patch_args_new = {
    "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
    "label": "Affine transformed constrained zonotope",
}
C = Polytope(V=spread_points_on_a_unit_sphere(2, 7)[0])
C = ConstrainedZonotope(polytope=C)
new_C = [[2, 0], [0, 0.5]] @ C
ax, _, _ = C.plot(patch_args=patch_args_old)
new_C.plot(ax=ax, patch_args=dict(patch_args_new, **{"alpha": 0.5}))
ax.legend(loc="lower right", bbox_to_anchor=(1, 0))
ax.grid()
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Grow in x, but shrink in y");

Negation and scaling by a scalar¶

In [15]:
negated_C = -C
scaled_C = 0.2 * C
ax = C.plot(
    patch_args={
        "alpha": 0.4,
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
        "label": "Original",
    }
)[0]
negated_C.plot(
    ax=ax,
    patch_args={
        "alpha": 0.4,
        "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
        "label": "-C",
    },
)
scaled_C.plot(
    ax=ax,
    patch_args={"alpha": 0.4, "facecolor": "lightgreen", "label": "0.2C"},
)
ax.legend()
Out[15]:
<matplotlib.legend.Legend at 0x7bb73239cd90>

Inverse affine transformation¶

We can also undo an affine transformation when it is invertible.

Note that the ConstrainedZonotope class overloads the operator @ and - to accomplish these transformations. Alternatively, you can use inverse_affine_map and minus methods.

We use == to check if constrained zonotopes are equal. For more details, see containment.

In [16]:
old_constrained_zonotope_from_new_constrained_zonotope = (
    new_constrained_zonotope_from_affine_transformation - translation
) @ rotation_matrix
old_constrained_zonotope_from_new_constrained_zonotope_explicit = (
    new_constrained_zonotope_from_affine_transformation.minus(
        translation
    ).inverse_affine_map_under_invertible_matrix(rotation_matrix)
)
print(
    "Are constrained zonotopes obtained from these methods the same?",
    old_constrained_zonotope_from_new_constrained_zonotope_explicit
    == old_constrained_zonotope_from_new_constrained_zonotope,
)
print(
    "Did we get the old constrained zonotope back?",
    old_constrained_zonotope_from_new_constrained_zonotope_explicit == CZ,
)

patch_args_new = {
    "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
    "label": "Affine transformed constrained zonotope",
}
patch_args_new_undo = {
    "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
    "label": "Undo the affine transformation",
}
ax, _, _ = new_constrained_zonotope_from_affine_transformation.plot(
    patch_args=patch_args_new
)
old_constrained_zonotope_from_new_constrained_zonotope_explicit.plot(
    ax=ax, patch_args=patch_args_new_undo
)
ax.legend()
ax.set_title("Inverse affine transformation");

Are constrained zonotopes obtained from these methods the same? True

Did we get the old constrained zonotope back? True

Projection of a constrained zonotope¶

Given a constrained zonotope $\mathcal{P}\subset\mathbb{R}^n$, the (orthogonal) projection of a constrained zonotopes is defined as $$\mathcal{R} = \{r \in \mathbb{R}^{m}: \exists v \in \mathbb{R}^{n - m},\ \mathrm{Lift}(r,v)\in \mathcal{P}\}$$ Here, $m = \mathcal{P}.dim - \text{length}(\text{project\_away\_dimensions})$, and $\mathrm{Lift}(r,v)$ lifts (undoes the projection) by combining $r$ and $v$ to produce a $n$-dimensional vector.

You can rotate the 3D plots using your mouse.

WARNING: You will see issues with matplotlib's rendering when you rotate the figure (see the second plot). See https://matplotlib.org/stable/api/toolkits/mplot3d/faq.html#my-3d-plot-doesn-t-look-right-at-certain-viewing-angles for more details. When a visual comparison of the plots is desired, it may be better to just plot the frame by setting patch_args['facecolor]=None. This issue does not appear for 2D plotting.

In [17]:
L1_norm_ball_polytope = Polytope(V=np.vstack((np.eye(3), -np.eye(3))))
L1_norm_ball = ConstrainedZonotope(polytope=L1_norm_ball_polytope)
simplex_in_R_pos_polytope = Polytope(V=np.vstack((np.eye(3), [0, 0, 0])))
simplex_in_R_pos = ConstrainedZonotope(polytope=simplex_in_R_pos_polytope)
constrained_zonotope_list = [
    L1_norm_ball,
    simplex_in_R_pos,
    simplex_in_R_pos,
]
view_dict = [
    {"elev": 21, "azim": -120},
    {"elev": 21, "azim": -35},
    {"elev": 21, "azim": -35},
]
facecolor_3d_list = ["lightblue", "lightblue", None]
for C, view_dict, facecolor_3d in zip(
    constrained_zonotope_list, view_dict, facecolor_3d_list
):
    fig = plt.figure()
    ax = fig.add_subplot(1, 2, 1, projection="3d")
    ax.view_init(**view_dict)
    C.plot(
        ax=ax,
        direction_vectors=dir_vectors,
        patch_args={"facecolor": facecolor_3d},
    )
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.set_aspect("equal")
    ax.set_title("3D constrained zonotope")
    P_projection = C.projection(project_away_dims=2)
    ax2d = fig.add_subplot(1, 2, 2)
    P_projection.plot(ax=ax2d)
    ax2d.set_xlabel("x")
    ax2d.set_ylabel("y")
    ax2d.set_aspect("equal")
    ax2d.set_title("Projection on x-y")
    plt.subplots_adjust(wspace=1)

Project a point on to a constrained zonotope¶

Given a constrained zonotope $\mathcal{P}\subset\mathbb{R}^n$ and a point $y\in\mathbb{R}^n$, the projection of the point on to $\mathcal{P}$ is defined as the optimal solution to the optimization problem $x^\ast = \arg\inf_{x\in\mathcal{P}} \|x - y\|_2$.

The ConstrainedZonotope class also provides closest_point and distance methods that are wrapper to project method. Projection can also be computed using 1 and inf norms.

In [18]:
L1_norm_ball_polytope = Polytope(V=np.vstack((np.eye(2), -np.eye(2))))
L1_norm_ball = ConstrainedZonotope(polytope=L1_norm_ball_polytope)
point_to_project = np.array([1, 0.5])
projected_point, distance_to_constrained_zonotope = L1_norm_ball.project(
    point_to_project
)
print(
    f"Distance of point to project {np.array2string(point_to_project)} "
    f"to the constrained zonotope: {distance_to_constrained_zonotope[0]:1.2f}"
)
print(
    "Did closest_point return the same point?",
    np.allclose(
        projected_point, L1_norm_ball.closest_point(point_to_project)
    ),
)
print(
    "Did distance return the same distance?",
    np.allclose(
        distance_to_constrained_zonotope,
        L1_norm_ball.distance(point_to_project),
    ),
)

patch_args = {
    "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
    "label": "L1 norm ball",
}
ax, _, _ = L1_norm_ball.plot(patch_args=patch_args)
ax.scatter(
    point_to_project[0],
    point_to_project[1],
    color="red",
    label="Point to project",
)
ax.scatter(
    projected_point[0, 0],
    projected_point[0, 1],
    color="k",
    label="Projected point",
)
ax.plot(
    [point_to_project[0], projected_point[0, 0]],
    [point_to_project[1], projected_point[0, 1]],
    "k--",
    label="Distance",
)
ax.legend(bbox_to_anchor=(1, 1))
ax.set_aspect("equal")
ax.grid()
plt.subplots_adjust(right=0.7)
ax.set_title("Projection of a point on a constrained zonotope");

Distance of point to project [1.  0.5] to the constrained zonotope: 0.35
Did closest_point return the same point? True
Did distance return the same distance? True

Containment: Check if a point is in a constrained zonotope¶

Given a constrained zonotope $\mathcal{P}\subset\mathbb{R}^n$ and a point $y\in\mathbb{R}^n$, we can check if $y\in\mathcal{P}$.

The ConstrainedZonotope class also overloads in, <, <= operator for this purpose.

In [19]:
L1_norm_ball_polytope = Polytope(V=np.vstack((np.eye(2), -np.eye(2))))
L1_norm_ball = ConstrainedZonotope(polytope=L1_norm_ball_polytope)
point_to_project = np.array([1, 0.5])
projected_point = L1_norm_ball.closest_point(point_to_project)
print(
    "Does the unit l1-norm ball contain point_to_project? "
    f"{L1_norm_ball.contains(point_to_project)}"
)
print(
    "Is point_to_project in the unit l1-norm ball? "
    f"{point_to_project in L1_norm_ball}"
)
print(
    "Is point_to_project <= the unit l1-norm ball? "
    f"{point_to_project <= L1_norm_ball}"
)
print(
    "Is point_to_project < the unit l1-norm ball? "
    f"{point_to_project < L1_norm_ball}"
)

print(
    "\nDoes the unit l1-norm ball contain projected_point? "
    f"{L1_norm_ball.contains(projected_point)}"
)
print(
    "Is projected_point in the unit l1-norm ball? "
    f"{projected_point in L1_norm_ball}"
)
print(
    "Is projected_point <= the unit l1-norm ball? "
    f"{projected_point <= L1_norm_ball}"
)
print(
    "Is projected_point < the unit l1-norm ball? "
    f"{projected_point < L1_norm_ball}"
)

Does the unit l1-norm ball contain point_to_project? False
Is point_to_project in the unit l1-norm ball? False
Is point_to_project <= the unit l1-norm ball? False
Is point_to_project < the unit l1-norm ball? False

Does the unit l1-norm ball contain projected_point? True
Is projected_point in the unit l1-norm ball? True
Is projected_point <= the unit l1-norm ball? True
Is projected_point < the unit l1-norm ball? True

Support function¶

Given a constrained zonotope $\mathcal{P}\subset\mathbb{R}^n$ and a vector $\eta\in\mathbb{R}^n$, the support function of the constrained zonotope is defined as $\rho_{\mathcal{P}}=\sup_{x\in\mathcal{P}} \eta^\top x$, and its optimal solution as the support vector. Recall that $$\mathcal{P}\subset\{x:\eta^\top x \leq \rho_{\mathcal{P}}(\eta)\}$$ for any direction $\eta\in\mathbb{R}^n$.

The ConstrainedZonotope class also provides extreme as a wrapper method to the support method to directly compute the support vectors.

In [20]:
P = Polytope(V=spread_points_on_a_unit_sphere(2, 7)[0])
C = np.diag([0.5, 0.75]) @ ConstrainedZonotope(polytope=P)
direction_vectors = np.array([[0, -1], [0, 1], [1, 0.5], [-2.5, -1]])
support_evaluations, support_points = C.support(direction_vectors)

outer_approximating_polytope = Polytope(
    A=direction_vectors, b=support_evaluations
)
ax, _, _ = C.plot(
    patch_args={"facecolor": "lightgreen", "label": "ConstrainedZonotope"}
)
outer_approx_patch_args = {
    "facecolor": "lightblue",
    "alpha": 0.4,
    "label": "Outer-approx. polytope",
}
outer_approximating_polytope.plot(
    ax=ax, patch_args=outer_approx_patch_args
)
ax.scatter(
    support_points[:, 0],
    support_points[:, 1],
    color="red",
    label="Support points",
)
ax.set_aspect("equal")
ax.grid()
ax.legend(bbox_to_anchor=(0.8, 1))
plt.subplots_adjust(right=0.6)
print(
    "Did the extreme function return the same support points?",
    np.allclose(C.extreme(direction_vectors), support_points),
)
ax.set_title("Support function/vector evaluation");

Did the extreme function return the same support points? True

Chebyshev centering¶

We can compute the Chebyshev center (approximately) of a constrained zonotope.

In [21]:
CZ1 = ConstrainedZonotope(
    polytope=np.diag([2, 0.75])
    @ Polytope(V=spread_points_on_a_unit_sphere(2, 7)[0])
)
CZ2 = ConstrainedZonotope(
    G=np.hstack((np.eye(2), 3 * np.eye(2))),
    c=[10, 10],
    Ae=[1, 1, 1, 1],
    be=2.85,
)
for P in [CZ1, CZ2]:
    c, r = P.chebyshev_centering()
    ax, _, _ = P.plot(
        patch_args={"facecolor": "lightgreen", "label": "Polytope"}
    )
    cheby_ball = Ellipsoid(c=c, r=r)
    cheby_ball.plot(
        ax=ax,
        patch_args={"facecolor": None, "label": "Chebyshev ball"},
        center_args={"color": "k", "label": "Chebyshev center"},
    )
    ax.set_aspect("equal")
    ax.grid()
    ax.legend(loc="lower right", bbox_to_anchor=(1.2, 0))
    plt.subplots_adjust(right=0.8)
    ax.set_title("Chebyshev centering");

/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_unary.py:78: UserWarning: This function returns a sub-optimal but feasible solution for Chebyshev centering.
  warnings.warn("This function returns a sub-optimal but feasible solution for Chebyshev centering.", UserWarning)
In [22]:
P = Polytope(V=spread_points_on_a_unit_sphere(2, 7)[0])
C = np.diag([0.5, 0.75]) @ ConstrainedZonotope(polytope=P)
direction_vectors = np.array([[0, -1], [0, 1], [1, 0.5], [-2.5, -1]])
support_evaluations, support_points = C.support(direction_vectors)

outer_approximating_polytope = Polytope(
    A=direction_vectors, b=support_evaluations
)
ax, _, _ = C.plot(
    patch_args={"facecolor": "lightgreen", "label": "ConstrainedZonotope"}
)
outer_approx_patch_args = {
    "facecolor": "lightblue",
    "alpha": 0.4,
    "label": "Outer-approx. polytope",
}
outer_approximating_polytope.plot(
    ax=ax, patch_args=outer_approx_patch_args
)
ax.scatter(
    support_points[:, 0],
    support_points[:, 1],
    color="red",
    label="Support points",
)
ax.set_aspect("equal")
ax.grid()
ax.legend(bbox_to_anchor=(0.8, 1))
plt.subplots_adjust(right=0.6)
print(
    "Did the extreme function return the same support points?",
    np.allclose(C.extreme(direction_vectors), support_points),
)
ax.set_title("Support function/vector evaluation");

Did the extreme function return the same support points? True

Maximum volume inscribing ellipsoid¶

We can compute the maximum volume inscribing ellipsoid center (approximately) of a constrained zonotope.

In [23]:
CZ1 = ConstrainedZonotope(
    polytope=np.diag([2, 0.75])
    @ Polytope(V=spread_points_on_a_unit_sphere(2, 7)[0])
)
CZ2 = ConstrainedZonotope(
    G=np.hstack((np.eye(2), 3 * np.eye(2))),
    c=[10, 10],
    Ae=[1, 1, 1, 1],
    be=2.85,
)
for P in [CZ1, CZ2]:
    c, Q, _ = P.maximum_volume_inscribing_ellipsoid()
    ax, _, _ = P.plot(
        patch_args={"facecolor": "lightgreen", "label": "Polytope"}
    )
    mvi_ellipsoid = Ellipsoid(c=c, Q=Q)
    mvi_ellipsoid.plot(
        ax=ax,
        patch_args={"facecolor": None, "label": "MVIE"},
        center_args={"color": "k", "label": "MVIE center"},
    )
    ax.set_aspect("equal")
    ax.grid()
    ax.legend(loc="lower right", bbox_to_anchor=(1.2, 0))
    plt.subplots_adjust(right=0.8)
    ax.set_title("Maximum volume inscribing ellipsoid (MVIE)");

/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_unary.py:146: UserWarning: This function returns a sub-optimal but feasible solution for maximum volume inscribed ellipsoid-based centering.
  warnings.warn(

Minimum volume circumscribing rectangle¶

Given a constrained zonotope $\mathcal{P}\subset\mathbb{R}^n$, we now compute the minimum volume rectangle that completely covers the constrained zonotope. Recall that rectangles are parameterized by their lower and upper bounds $(l, u)$, i.e., $\mathcal{R}(l, u)=\{x| l \leq x \leq u \}$.

We compute the ellipsoid by solving the following optimization problem $\inf_{\mathcal{P}\subseteq\mathcal{R}(l, u)} \mathrm{Vol}(\mathcal{R}(l,u))$. This optimization problem can be solved via linear programming.

In [24]:
P = Polytope(V=[[0, 0], [1, 0], [0, 2]])
C = (
    rotation_matrix
    @ np.diag([1.5, 1.25])
    @ ConstrainedZonotope(polytope=P)
)
l, u = C.minimum_volume_circumscribing_rectangle()
min_vol_circ_rect = Polytope(lb=l, ub=u)
print(f"Rectangle lower bound is: {np.array2string(l)}")
print(f"Rectangle upper bound is: {np.array2string(u)}")
print(f"Rectangle volume is: {min_vol_circ_rect.volume():1.2f}")

P_patch_args = {"facecolor": "lightgreen"}
min_vol_circ_rect_patch_args = {"facecolor": "lightblue"}
ax, hE, _ = min_vol_circ_rect.plot(patch_args=min_vol_circ_rect_patch_args)
_, hP, _ = C.plot(ax=ax, patch_args=P_patch_args)
ax.set_aspect("equal")
ax.grid()
ax.legend(
    [hP, hE],
    ["ConstrainedZonotope", "Min. vol. rectangle"],
    bbox_to_anchor=(1, 0.75),
)
plt.subplots_adjust(right=0.6)
ax.set_title("Minimum volume circumscribing rectangle");

Rectangle lower bound is: [-1.76776695e+00  5.97257454e-10]
Rectangle upper bound is: [1.06066017 1.76776695]
Rectangle volume is: 5.00

Interior point¶

Unlike zonotopes, the constrained zonotope c need be an interior point.

In [25]:
CZ1 = ConstrainedZonotope(
    polytope=np.diag([2, 0.75])
    @ Polytope(V=spread_points_on_a_unit_sphere(2, 7)[0])
)
CZ2 = ConstrainedZonotope(
    G=np.hstack((np.eye(2), 3 * np.eye(2))),
    c=[10, 10],
    Ae=[1, 1, 1, 1],
    be=2.85,
)
for P in [CZ1, CZ2]:
    chebyshev_center = P.interior_point()
    ax, _, _ = P.plot(
        patch_args={"facecolor": "lightgreen", "label": "Polytope"}
    )
    ax.scatter(
        chebyshev_center[0],
        chebyshev_center[1],
        color="red",
        label="Chebyshev center",
    )
    ax.scatter(
        P.c[0],
        P.c[1],
        color="k",
        label="CZ.c",
    )
    ax.set_aspect("equal")
    ax.grid()
    ax.legend(loc="lower right", bbox_to_anchor=(1.2, 0))
    plt.subplots_adjust(right=0.8)
    ax.set_title("Interior point");

Volume computation¶

We can scipy.spatial.ConvexHull to compute the volume of a constrained zonotope.

We perform a simple sanity check where we compute the volume a unit l1-norm ball. From the figures above, it is evident that the "rhombus" consists of four isosceles right angled triangles with base as $1$.

In [26]:
L1_norm_ball_polytope = Polytope(V=np.vstack((np.eye(2), -np.eye(2))))
L1_norm_ball = ConstrainedZonotope(polytope=L1_norm_ball_polytope)
volume_of_an_isosceles_triangle_with_unit_base = 0.5
volume_of_unit_l1_norm_ball = approximate_volume_from_grid(
    L1_norm_ball, 0.1
)
print(volume_of_unit_l1_norm_ball)
print(
    "Is the volume of the unit l1-norm ball (approximately) 2?",
    abs(
        volume_of_unit_l1_norm_ball
        - 4 * volume_of_an_isosceles_triangle_with_unit_base
    )
    <= 0.25,
)

2.1999999987300662
Is the volume of the unit l1-norm ball (approximately) 2? True

Cartesian product with itself¶

In [27]:
cartesian_Z = ConstrainedZonotope(lb=[-1, -1], ub=[1, 1])
cartesian_Z_4D = cartesian_Z**2
print(
    "Is cartesian_Z_4D a cuboid?",
    cartesian_Z_4D == ConstrainedZonotope(c=[0, 0, 0, 0], h=1),
)

Is cartesian_Z_4D a cuboid? True

Operations involving two constrained zonotopes¶

In [28]:
OLD_CONSTRAINED_ZONOTOPE_2 = ConstrainedZonotope(c=[0, 0], h=2) + [-2, 4]
OLD_POLYTOPE_2 = Polytope(c=[0, 0], h=2) + [-2, 4]
OLD_CONSTRAINED_ZONOTOPE_COLOR_2 = "lightsalmon"

Minkowski sum¶

We perform a Minkowski sum operation on two constrained zonotopes. Specifically, given two constrained zonotopes $\mathcal{P}_1$ and $\mathcal{P}_2$, we compute $$\mathcal{Q}=\{p_1+p_2\ |\ p_1\in\mathcal{P}_1,p_2\in\mathcal{P}_2\}.$$

Sanity check: The convex hull of the constrained zonotopes obtained by translating $\mathcal{P}_1$ to $\mathcal{P}_2$'s vertices coincides with $\mathcal{Q}$.

Note: This code snippet also shows the advantage of obtaining handles. In line 15, if label is passed to patch_args, we would have as many legend entries as the number of vertices in OLD_POLYTOPE_2.V.

In [29]:
new_constrained_zonotope = CZ + OLD_CONSTRAINED_ZONOTOPE_2

fig, ax = plt.subplots(1, 2)
CZ.plot(
    ax=ax[0],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
        "label": "P1",
    },
)
OLD_CONSTRAINED_ZONOTOPE_2.plot(
    ax=ax[0],
    patch_args={
        "zorder": 10,
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_2,
        "label": "P2",
    },
)
new_constrained_zonotope.plot(
    ax=ax[0],
    patch_args={"facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR, "label": "Q"},
)
ax[0].legend(loc="best")
ax[0].set_title("Minkowski sum Q = P1 $\oplus$ P2")

CZ.plot(
    ax=ax[1], patch_args={"facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1}
)
new_constrained_zonotope.plot(
    ax=ax[1],
    patch_args={
        "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
        "alpha": 0.9,
        "label": r"Q=P1 $\oplus$ P2.v",
    },
)
OLD_CONSTRAINED_ZONOTOPE_2.plot(
    ax=ax[1],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_2,
        "alpha": 0.2,
    },
    vertex_args={"s": 50, "color": "r"},
)
for v in OLD_POLYTOPE_2.V:
    temp_constrained_zonotope = v + CZ
    _, h, _ = temp_constrained_zonotope.plot(
        ax=ax[1], patch_args={"facecolor": "k", "alpha": 0.8}
    )
handles, labels = ax[1].get_legend_handles_labels()
ax[1].legend([*handles, h], [*labels, "P2.v"], loc="best")
ax[1].set_title("Sanity check");

/tmp/ipykernel_236822/1984422472.py:3: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
  fig, ax = plt.subplots(1, 2)

Pontryagin difference¶

We perform a Pontryagin difference operation on two constrained zonotopes. Specifically, given two constrained zonotopes $\mathcal{P}_1$ and $\mathcal{P}_2$, we compute $$\mathcal{Q}=\mathcal{P}_2 \ominus \mathcal{P}_1 = \{p_2\ |\ \{p_2\} \oplus \mathcal{P}_1 \subseteq\mathcal{P}_2\} = \{p_2\ |\ \forall p_1 \in\mathcal{P}_1,\ p_1 + p_2 \in\mathcal{P}_2\}.$$

Sanity check: 1) $\mathcal{Q}\subseteq \mathcal{P}_2$, and 2) Translating $\mathcal{P}_1$ by $\mathcal{Q}$'s vertices still keeps it inside $\mathcal{P}_2$.

In [30]:
P_from_V = Polytope(V=[[0, 0], [0, 1], [2, 0]])
CZ = ConstrainedZonotope(polytope=P_from_V)
Z = 0.1 * ConstrainedZonotope(c=[0, 0], G=np.eye(2))
CZ_minus_Z = CZ - Z
CZ_minus_Z_polytope = CZ_minus_Z.polytopic_inner_approximation()
CZ_minus_Z_polytope.minimize_V_rep()

fig, ax = plt.subplots(1, 2)
CZ.plot(
    ax=ax[0],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
        "label": "C",
    },
)
Z.plot(
    ax=ax[0],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_2,
        "label": "Z",
    },
)
CZ_minus_Z.plot(
    ax=ax[0],
    patch_args={
        "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
        "label": r"Q=C $\ominus$ Z",
    },
)
ax[0].legend(loc="best")
ax[0].set_title("Pontryagin difference")

CZ.plot(
    ax=ax[1],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_2,
        "label": "C",
    },
)
CZ_minus_Z.plot(
    ax=ax[1],
    patch_args={
        "facecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
        "label": "Q = C $\ominus$ Z",
    },
)
for v in CZ_minus_Z_polytope.V:
    temp_polytope = v + Z
    _, h, _ = temp_polytope.plot(
        ax=ax[1], patch_args={"facecolor": "k", "alpha": 0.8}
    )
handles, labels = ax[1].get_legend_handles_labels()
ax[1].legend([*handles, h], [*labels, r"Z $\oplus$ Q.v"], loc="best")
ax[1].set_title("Sanity check");

/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_binary.py:644: UserWarning: This function returns an inner-approximation of the Pontryagin difference with a zonotope.
  warnings.warn(

Intersection¶

We perform an intersection operation on two constrained zonotopes. Specifically, given two constrained zonotopes $\mathcal{P}_1$ and $\mathcal{P}_2$, we compute $$\mathcal{Q}=\{p\ |\ p\in\mathcal{P}_1,p\in\mathcal{P}_2\}.$$

To obtain a non-empty intersection, we first translate $\mathcal{P}_1$ to obtain $\mathcal{P}_1^\dagger$.

Sanity check: Defining $\mathcal{Q}=\mathcal{P}_1^\dagger\cap\mathcal{P}_2$, we have 1) $\mathcal{Q}\subseteq\mathcal{P}_1^\dagger$ and 2) $\mathcal{Q}\subseteq\mathcal{P}_2$.

You can also compute intersection with affine sets, halfspaces, as well as intersections with other constrained zonotopes under an inverse-affine maps.

In [31]:
translate_by_vector = (-1, 3)
P4 = CZ + translate_by_vector
new_constrained_zonotope = OLD_CONSTRAINED_ZONOTOPE_2.intersection(P4)

fig = plt.figure()
ax = plt.gca()
CZ.plot(
    ax=ax,
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
        "label": "P1",
    },
)
OLD_CONSTRAINED_ZONOTOPE_2.plot(
    ax=ax,
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_2,
        "label": "P2",
    },
)
P4.plot(
    ax=ax,
    patch_args={
        "facecolor": "yellow",
        "linewidth": 3,
        "label": r"P1$^\dagger$",
    },
)
new_constrained_zonotope.plot(
    ax=ax,
    patch_args={
        "edgecolor": NEW_CONSTRAINED_ZONOTOPE_COLOR,
        "fill": False,
        "linewidth": 1,
        "label": "Q",
    },
)
ax.legend(loc="best")
ax.set_title(r"Intersection Q=P1$^\dagger\cap$ P2");

Intersection under inverse affine map¶

Given constrained zonotopes $\mathcal{P}=\{x|-1\leq x \leq 1\}\subset\mathbb{R}^3$ and $\mathcal{Q}=\{x|-1\leq x - 1\leq 1\}\subset\mathbb{R}^2$ and matrix $R=[I_2, 0_{2\times 1}]$, we compute $\mathcal{P}\cap_R\mathcal{Q}=\{x\in\mathcal{P}|Rx\in\mathcal{Q}\}=\{x|-1\leq x \leq 1, [I_2, 0_{2\times 1}]x\in\mathcal{Q}\}=\{x|-1\leq x \leq 1, 0\leq x_i \leq 2\text{ for }i\in\{1,2\}\}$. Unlike $\mathcal{P}\cap(\mathcal{Q}@R)$, this function does not require $R$ to be square.

In [32]:
C_3D = ConstrainedZonotope(c=[0, 0, 0], h=1)
ax, _, _ = C_3D.plot(
    direction_vectors=dir_vectors, patch_args={"label": "C_3D"}
)
C_3D.intersection_under_inverse_affine_map(
    ConstrainedZonotope(c=[1, 1], h=1), [[1, 0, 0], [0, 1, 0]]
).plot(
    ax=ax,
    direction_vectors=dir_vectors,
    patch_args={"facecolor": "yellow", "label": "C_3D after intersection"},
)
ax.legend()
ax.view_init(elev=15, azim=-24)
ax.set_title("Intersection under inverse affine map")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z");

Intersection with a halfspace¶

In [33]:
ax, _, _ = C_3D.plot(
    direction_vectors=dir_vectors, patch_args={"label": "C_3D"}
)
C_3D.intersection_with_halfspaces(A=[1, 1, 1], b=0).plot(
    ax=ax,
    direction_vectors=dir_vectors,
    patch_args={"facecolor": "yellow", "label": "C_3D after intersection"},
)
ax.legend()
ax.view_init(elev=15, azim=-24)
ax.set_title("Intersection with a halfspace")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z");

Intersection with an affine set¶

In [34]:
ax, _, _ = C_3D.plot(
    direction_vectors=dir_vectors, patch_args={"label": "C_3D"}
)
C_3D.intersection_with_affine_set(Ae=[1, 1, 1], be=0).plot(
    ax=ax,
    direction_vectors=dir_vectors,
    patch_args={"facecolor": "yellow", "label": "C_3D after intersection"},
)
ax.legend()
ax.view_init(elev=15, azim=-24)
ax.set_title("Intersection with affine set")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z");

Slice¶

In [35]:
ax, _, _ = C_3D.plot(
    direction_vectors=dir_vectors, patch_args={"label": "C_3D"}
)
C_3D.slice(dims=2, constants=0.5).plot(
    ax=ax,
    direction_vectors=dir_vectors,
    patch_args={"facecolor": "yellow", "label": "C_3D after slicing"},
)
ax.legend()
ax.view_init(elev=15, azim=-24)
ax.set_title("Slicing")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z");

Containment and equality check¶

We can see if a constrained zonotope covers another --- given two constrained zonotopes $\mathcal{P}_1$ and $\mathcal{P}_2$, we can check if $\mathcal{P}_1\overset{?}{\subseteq}\mathcal{P}_2$. pycvxset provides a method contains to do so, but also overloads binary operators <=, <, >=, >, in for the same. Finally, we can also check if constrained zonotopes are equal using ==.

You can also check the constrained zonotope contains another ellipsoid or a constrained zonotope as well.

In [36]:
P5 = 0.5 * (CZ - CZ.c) + CZ.c
fig, ax = plt.subplots(1, 2)
CZ.plot(
    ax=ax[0],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_1,
        "label": "P1",
    },
)
P5.plot(ax=ax[0], patch_args={"facecolor": "black", "label": "P5"})
ax[0].legend(loc="best")
print(f"Does P1 cover P5? {CZ.contains(P5)}")
print(f"Is P1 >= P5? {CZ >= P5}")
print(f"Is P1 > P5? {CZ > P5}")
print(f"Is P5 in P1? {P5 in CZ}")
print(f"Is P1 <= P5? {CZ <= P5}")
print(f"Is P1 < P5? {CZ < P5}")
print(f"Is P1 in P5? {CZ in P5}")
OLD_CONSTRAINED_ZONOTOPE_2.plot(
    ax=ax[1],
    patch_args={
        "facecolor": OLD_CONSTRAINED_ZONOTOPE_COLOR_2,
        "label": "P2",
    },
)
P4.plot(ax=ax[1], patch_args={"facecolor": "grey", "label": "P4"})
ax[1].legend(loc="best")
print(f"\n\nDoes P2 cover P4? {OLD_CONSTRAINED_ZONOTOPE_2.contains(P4)}")
print(f"Is P2 >= P4? {OLD_CONSTRAINED_ZONOTOPE_2 >= P4}")
print(f"Is P2 > P4? {OLD_CONSTRAINED_ZONOTOPE_2 > P4}")
print(f"Is P4 in P2? {P4 in OLD_CONSTRAINED_ZONOTOPE_2}")
print(f"Is P2 <= P4? {OLD_CONSTRAINED_ZONOTOPE_2 <= P4}")
print(f"Is P2 < P4? {OLD_CONSTRAINED_ZONOTOPE_2 < P4}")
print(f"Is P2 in P4? {OLD_CONSTRAINED_ZONOTOPE_2 in P4}")
plt.suptitle("Containment check for two sets of constrained zonotopes");

Does P1 cover P5? True
Is P1 >= P5? True

Is P1 > P5? True
Is P5 in P1? True
Is P1 <= P5? False

Is P1 < P5? False
Is P1 in P5? False


Does P2 cover P4? False
Is P2 >= P4? False
Is P2 > P4? False
Is P4 in P2? False
Is P2 <= P4? False
Is P2 < P4? False
Is P2 in P4? False