Constrained polynomial zonotopes

We introduce constrained polynomial zonotopes, a novel non-convex set representation that is closed under linear map, Minkowski sum, Cartesian product, convex hull, intersection, union, and quadratic as well as higher-order maps. We show that the computational complexity of the above-mentioned set operations for constrained polynomial zonotopes is at most polynomial in the representation size. The fact that constrained polynomial zonotopes are generalizations of zonotopes, polytopes, polynomial zonotopes, Taylor models, and ellipsoids further substantiates the relevance of this new set representation. In addition, the conversion from other set representations to constrained polynomial zonotopes is at most polynomial with respect to the dimension, and we present efficient methods for representation size reduction and for enclosing constrained polynomial zonotopes by simpler set representations.


Introduction
Many applications like, e.g., controller synthesis, state estimation, and formal verification, are based on algorithms that compute with sets [7,11,12,28].The performance of these algorithms therefore mainly depends on efficient set representations.Ideally, a set representation is not only closed under all relevant set operations, but can also compute these efficiently.We introduce constrained polynomial zonotopes, a novel non-convex set representation that is closed under linear map, Minkowski sum, Cartesian product, convex hull, Table 1 Relation between set representations and set operations.The symbol √ indicates that the set representation is closed under the corresponding set operation and a closed-form expression for the computation exists, where ( √ ) indicates that this holds for linear maps represented by invertible matrices only.The symbol − indicates that the set representation is closed under the corresponding set operation, but no closed-form expression for the computation is known (iterative algorithms, such as Fourier-Motzkin elimination, are not counted as closed-form expressions).The symbol × indicates that the set representation is not closed under the corresponding set operation.An overview for the computational complexity of set operations is provided in [5, Table 1].

Set Representation Lin. Map
Mink.Sum Cart.Prod.

Quad. Map
Intersection Union are zonotope bundles [6] and constrained zonotopes [29], which are both able to represent any bounded polytope.Constrained zonotopes additionally consider linear equality constraints for the zonotope factors, whereas zonotope bundles represent the set implicitly by the intersection of several zonotopes.Two special cases of zonotopes are parallelotopes, which are zonotopes with linearly independent generators, and multi-dimensional intervals.Since intervals are not closed under linear map, algorithms computing with intervals often split them to obtain a desired accuracy [18].Common non-convex set representations are star sets, level sets, Taylor models, and polynomial zonotopes.The concept of star sets [8,13] is similar to the one of constrained zonotopes, but logical predicates instead of linear equality constraints are used to constrain the values of the zonotope factors.Level sets of nonlinear functions [25] can represent any shape.While star sets and level sets are very expressive (see Fig. 1), it is for many of the relevant operations unclear how they are computed (see Table 1).Taylor models [24] consist of a polynomial and an interval remainder part.A set representation that is very similar to Taylor models are polynomial zonotopes, which were first introduced in [2].A computationally efficient sparse representation of polynomial zonotopes was recently proposed in [21].Due to their polynomial nature, Taylor models and polynomial zonotopes are both closed under quadratic and higher-order maps (see Table 1).
In this work we introduce constrained polynomial zonotopes, a novel non-convex set representation that combines the concept of adding equality constraints for the zonotope factors used by constrained zonotopes [29] with the sparse polynomial zonotope representation in [21].Constrained polynomial zonotopes are closed under all relevant set operations (see Table 1) and can represent any set in Fig. 1, except star sets, level sets, and sets defined by their support function.As shown in Table 1, constrained polynomial zonotopes are the only set representation for which closed-form expressions for the calculation of all relevant set operations are known.

Notation and Assumptions
In the remainder of this work, we use the following notations: Sets are denoted by calligraphic letters, matrices by uppercase letters, and vectors by lowercase letters.Moreover, the set of natural numbers is denoted by N = {1, 2, . . .}, the set of natural numbers including zero is denoted by N 0 = {0, 1, 2, . . .}, and the set of real numbers is denoted by R. Given a set H = {h 1 , . . ., h n }, |H| = n denotes the cardinality of the set.Given a vector b ∈ R n , b (i) refers to the i-th entry.Likewise, given a matrix A ∈ R n×w , A (i,•) represents the i-th matrix row, A (•,j) the j-th column, and A (i,j) the j-th entry of matrix row i.Given a set of positive integer indices H = {h 1 , . . ., h |H| } with ∀i ∈ {1, . . ., |H|}, 1 ≤ h i ≤ w, notation A (•,H) is used for [A (•,h1) . . .A (•,h |H| ) ], where [C D] denotes the concatenation of two matrices C and D. The symbols 0 and 1 represent matrices and vectors of zeros and ones of proper dimension, and diag(a) returns a square matrix with a ∈ R n on the diagonal.The empty matrix is denoted by [ ] and the identity matrix of dimension n × n is denoted by I n ∈ R n×n .Moreover, we use the shorthand I = [l, u] for an n-dimensional interval I := {x ∈ R n | l (i) ≤ x (i) ≤ u (i) , i = 1, . . ., n}.For the derivation of computational complexity, we consider all binary operations, except concatenations; initializations are also not considered.

Definitions
Let us first provide some definitions that are important for the remainder of the paper.We begin with zonotopes: Definition 1 (Zonotope) [14,Def. 1] Given a constant offset c ∈ R n and a generator matrix G ∈ R n×p , a zonotope Z ⊂ R n is defined as The scalars α k are called factors and we use the shorthand Z = c, G Z .
Constrained zonotopes [29] can represent arbitrary bounded polytopes: Definition 2 (Constrained Zonotope) [29,Def. 3] Given a constant offset c ∈ R n , a generator matrix G ∈ R n×p , a constraint matrix A ∈ R m×p , and a constraint vector b ∈ R m , a constrained zonotope CZ ⊂ R n is defined as Polynomial zonotopes are a non-convex set representation first introduced in [2].We use the sparse representation of polynomial zonotopes [21]: a generator matrix G ∈ R n×h , and an exponent matrix E ∈ N p×h 0 , a polynomial zonotope PZ ⊂ R n is defined as In contrast to [21, Def.1], we explicitly do not integrate the constant offset c in G, and we do not consider independent generators since each polynomial zonotope with independent generators can be equivalently represented as a polynomial zonotope without independent generators [22,Prop. 1].We use the shorthand PZ = c, G, E P Z .
An ellipsoid is defined as follows: Definition 4 (Ellipsoid) [10, Eq. 2.3] Given a constant offset c ∈ R n and a symmetric and positive definite matrix Q ∈ R n×n , an ellipsoid E ⊂ R n is defined as In this paper we consider the standard set operations listed in Table 1.Given two sets S 1 , S 2 ⊂ R n , a set S 3 ⊂ R w , a matrix M ∈ R w×n , and a discrete set of matrices Q = {Q 1 , . . ., Q w } with Q i ∈ R n×n , i = 1, . . ., w, these operations are defined as follows: Linear map: Minkowski sum: Cartesian prod.: Intersection: Union: Moreover, we consider another set operation that we refer to as the linear combination of two sets: 1 Definition according to Caratheodory's theorem [9].This rather complex definition is required since constrained polynomial zonotopes can represent disjoint sets.
For convex sets, the convex hull and the linear combination are identical.However, for non-convex sets as considered in this paper, the two operations differ 2 .We consider both operations since for many algorithms, such as reachability analysis [1,Eq. (3.4)], it is sufficient to compute the linear combination instead of the convex hull.

Constrained Polynomial Zonotopes
In this section, we introduce constrained polynomial zonotopes (CPZs).A CPZ is constructed by adding polynomial equality constraints to a polynomial zonotope: and a constraint exponent matrix R ∈ N p×q 0 , a constrained polynomial zonotope is defined as The constrained polynomial zonotope is regular if the exponent matrix E and the constrained exponent matrix R do not contain duplicate columns or all-zero columns: The scalars α k are called factors, where the number of factors is p, the number of generators G (•,i) is h, the number of constraints is m, and the number of constraint generators A (•,i) is q.The order ρ = h+q n estimates the complexity of a constrained polynomial zonotope.We use the shorthand CPZ = c, G, E, A, b, R CP Z .All components of a set i have index i, e.g., the parameter p i , h i , m i , and q i as defined in Def. 5 belong to CPZ i .The quantity of scalar numbers µ required to store a CPZ is since c has n entries, G has nh entries, E has ph entries, A has mq entries, b has m entries, and R has pq entries.We call µ the representation size of the CPZ.Moreover, we call the polynomial zonotope PZ = c, G, E P Z corresponding to CPZ = c, G, E, A, b, R CP Z the constructing polynomial zonotope.For the derivation of the computational complexity of set operations with respect to the dimension n, we make the assumption that with a p , a h , a q , a m ∈ R ≥0 .This assumption is justified by the fact that one usually reduces the representation size to a desired upper bound when computing with CPZs. 2 The convex hull that we consider in our previous work on sparse polynomial zonotopes in [21] actually defines a linear combination since polynomial zonotopes are non-convex.
PSfrag replacements We demonstrate the concept of CPZs by an example: defines the set which is visualized in Fig. 2.

Preliminaries
We begin with some preliminary results that are required throughout this paper.

Identities
Let us first establish some identities that are useful for subsequent derivations.According to the definition of CPZs in Def. 5, it holds that and .

Transformation to a Regular Representation
Some set operations result in a CPZ that is not regular.We therefore introduce operations that transform a non-regular CPZ into a regular one.The compactGen operation returns a CPZ with a regular exponent matrix: where the operation uniqueColumns removes identical matrix columns until all columns are unique.
Proof For a CPZ where the exponent matrix E = [e e] consists of two identical columns e ∈ N p 0 , it holds that Summation of the generators for terms α with identical exponents therefore does not change the set, which proves that compactGen(CPZ) = CPZ.In addition, since the operation uniqueColomns removes all identical matrix columns and we add all-zero columns to the constant offset, it holds that the resulting exponent matrix E is regular according to Def. 5.
Complexity: We assume that the operation uniqueColumns in combination with the construction of the sets H j is implemented by first sorting the matrix columns, followed by an identification of identical neighbors.Moreover, we assume that in order to sort the matrix columns one first sorts the entries in the first row.For all columns with identical entries in the first row one then sorts the columns according to the entries in the second row.Since this process is continued for all p matrix rows and the complexity for sorting one row of the matrix The compactCon operation returns a CPZ with a regular constraint exponent matrix: the operation compactCon returns a representation of CPZ with a regular constraint exponent matrix and has complexity O(pq log(q) + mq): where the operation uniqueColumns removes identical matrix columns until all columns are unique.
Proof The proof is analogous to the proof for Prop. 1.

Lifted Polynomial Zonotopes
Finally, we introduce the lifted polynomial zonotope corresponding to a CPZ in the following lemma, which is inspired by [29,Prop. 3]: Proof With the definition of CPZs in Def. 5 we obtain According to Lemma 1, a CPZ can be interpreted as the intersection of the lifted polynomial zonotope PZ + with the subspace {x ∈ R n+m | x (n+1) , . . ., x (n+m) = 0}.Moreover, with the lifted polynomial zonotope we can transfer results for polynomial zonotopes to CPZs, as we demonstrate later.Potential redundancies in the lifted polynomial zonotope due to common columns in the exponent and the constraint exponent matrix can be removed using the compact operation for polynomial zonotopes in [21, Prop.2].

Rescaling
Later, in Sec.6 and Sec. 8, we describe how to enclose CPZs by other set representations and how to reduce the representation size of a CPZ by enclosing it with a simpler CPZ.The tightness of these enclosures mainly depends on the size of the corresponding constructing polynomial zonotope.Since the constraints often intersect only part of the factor hypercube α 1 , . . ., α p ∈ [ 1, 1], we can reduce the size of the constructing polynomial zonotope in advance to obtain tighter results.This can be achieved with a contractor: Definition 6 (Contractor) [18,Chapter 4.1] Given an interval I ⊂ R p and a vector field f : R p → R m which defines the constraint f (x) = 0, the operation contract returns an interval that satisfies contract f (x), I ⊆ I and ∀x ∈ I, f (x) = 0 ⇒ x ∈ contract f (x), I , so that it is guaranteed that all solutions for f (x) = 0 contained in I are also contained in the contracted interval. - PSfrag replacements x 2 α 3 Fig. 3 Visualization of rescaling for CPZ from Example 2 (red, right), where the corresponding constraint is visualized on the left.The constructing polynomial zonotope before rescaling is shown in blue, and the constructing polynomial zonotope after rescaling is shown in green.
There exist many sophisticated approaches for implementing a contractor, an overview of which is provided in [18,Chapter 4].
for the factors by applying a contractor to the polynomial constraint of the CPZ: Using the contracted domain [l, u], the CPZ can be equivalently represented as We show in Appendix B that this set can be represented as a CPZ.Let us demonstrate rescaling by an example: , which is visualized in Fig. 3.As depicted on the left side of Fig. 3, the constraint only intersects a small part of the factor domain α 1 , α 2 , α 3 ∈ [ 1, 1], so that the domain can be contracted to α 1 , α 2 , α 3 ∈ [ 1, 0].Rescaling therefore significantly reduces the size of the constructing polynomial zonotope, as visualized on the right side of Fig. 3.

Conversion from Other Set Representations
This section shows how other set representations can be converted to CPZs.

Taylor Models, Intervals, and Zonotopic Set Representations
Since a polynomial zonotope is simply a CPZ without constraints, the conversion is trivial in this case.For polynomial zonotopes that are defined with additional independent generators as in [21, Def.1], one can first convert the polynomial zonotope to a polynomial zonotope without independent generators using [22,Prop. 1].According to [21,Prop. 4], the set defined by a Taylor model can be equivalently represented as a polynomial zonotope.Moreover, according to [21,Prop. 3] any zonotope can be represented as a polynomial zonotope, and any interval can be represented as a zonotope [1, Prop.2.1].Finally, a constrained zonotope is a special case of a CPZ where all polynomial functions are linear, so the conversion is straightforward.In summary, we therefore obtain the following conversion rules: Zonotope: Constrained zonotope: Polynomial zonotope: The conversion of an interval has complexity O(n) with respect to the dimension n due to the summation and subtraction of the vectors l and u, while all other conversions have constant complexity O(1) since no computations are required.

Polytopes
There are two possibilities to represent a bounded polytope as a CPZ.According to [20] and [21, Theorem 1], every bounded polytope can be represented as a polynomial zonotope.Therefore, any bounded polytope can be converted to a CPZ by first representing it as a polynomial zonotope followed by a conversion of the polynomial zonotope to a CPZ using (17).Moreover, it holds according to [29, Theorem 1] that any bounded polytope can be represented as a constrained zonotope.Consequently, the second possibility for the conversion of a bounded polytope to a CPZ is to first represent the polytope as a constrained zonotope, and then convert the constrained zonotope to a CPZ using (16).Which of the two methods results in a more compact representation depends on the polytope.

Ellipsoids
Any ellipsoid can be converted to a CPZ: be equivalently represented by a CPZ: where the eigenvalues λ 1 , . . ., λ n , the matrix of eigenvalues D, and the matrix of eigenvectors V are obtained by the eigenvalue decomposition The complexity of the conversion is O(n 3 ).
Proof The matrices A, R and the vector b in ( 18) define the constraint Since Using the eigenvalue decomposition of the matrix Q from ( 19) it holds that since V is an orthonormal matrix satisfying V −1 = V T .Inserting (22) into the definition of an ellipsoid in Def. 4 yields We define the factors α k of the CPZ as Inserting ( 24) into ( 23) finally yields = c, G, E, A, b, R CP Z , which concludes the proof.
Complexity: Computation of the eigenvalue decomposition [26].The computation of G in (18) requires n 2 multiplications and the calculation of n square roots and therefore has complexity ).Since all other required operations are concatenations, the overall complexity results by adding the complexity of the eigenvalue decomposition and the complexity of computing G, which yields O(n 2 ) + O(n 3 ) = O(n 3 ).

Enclosure by Other Set Representations
To speed up computations, one often encloses sets by simpler set representations in setbased computing.In this section, we therefore show how to enclose CPZs by constrained zonotopes, polynomial zonotopes, zonotopes, and intervals.The over-approximation error for all enclosures can be reduced by applying rescaling as described in Sec.4.4 in advance.
To demonstrate the tightness of the enclosures, we use the CPZ as a running example throughout this section.

Constrained Zonotopes
We first show how to enclose a CPZ by a constrained zonotope: , the operation conZono returns a constrained zonotope that encloses CPZ: , where the compact operation as defined in [21, Prop.2] returns a regular polynomial zonotope and the zono operation as defined in [21, Prop.5] returns a zonotope that encloses a polynomial zonotope.The computational complexity is O(µ 2 ) with respect to the representation size µ and O(n 2 log(n)) with respect to the dimension n.
Proof To obtain an enclosing constrained zonotope we calculate a zonotope enclosure of the corresponding lifted polynomial zonotope as defined in Lemma 1. Back-transformation of the lifted zonotope to the original state space then yields an enclosing constrained zonotope: where we omitted the compact operation since it only changes the representation of the set, but not the set itself.
Complexity: Let n + = n + m, p + = p, and h + = h + q denote the dimension, the number of factors, and the number of generators of the lifted polynomial zonotope PZ + .According to [21, Prop.2], the compact operation for polynomial zonotopes has complexity [21,Prop. 5].The overall computational complexity is therefore which is O(n 2 log(n)) using (10).
The enclosing constrained zonotope for the CPZ in (25) is shown in Fig. 4.

Polynomial Zonotopes
Clearly, an enclosing polynomial zonotope for a CPZ can simply be obtained by dropping the constraints.However, this might yield large over-approximation errors.Another possibility is to reduce all constraints using Prop.13 introduced later in Sec. 8. Which method results in the tighter enclosure depends on the CPZ.The resulting enclosing polynomial zonotope for the CPZ in (25) obtained by dropping the constraints is visualized in Fig. 4.

Zonotopes and Intervals
An enclosure of a CPZ by a zonotope or interval can be computed using the previously presented enclosures by constrained zonotopes or polynomial zonotopes.For polynomial zonotopes, an enclosing zonotope can be computed using [21,Prop. 5], and an enclosing Fig. 4 Enlosing constrained zonotope (left) and enclosing polynomial zonotope (right) for CPZ in (25).
interval can be computed based on the support function enclosure in [21,Prop. 7].For constrained zonotopes, an enclosing zonotope can be calculated by reducing all constraints as described in [29,Sec. 4.2], and an enclosing interval can be computed using linear programming [27, Prop.1].

Set Operations
In this section, we derive closed-form expressions for all set operations introduced in Sec. 2 on CPZs.We begin with the linear map: Proposition 5 (Linear Map) Given CPZ = c, G, E, A, b, R CP Z ⊂ R n and a matrix M ∈ R w×n , the linear map is which has complexity O(wµ) with respect to the representation size µ and complexity O(wn 2 ) with respect to the dimension n.The resulting CPZ is regular if CPZ is regular.
Proof The result follows directly from inserting the definition of CPZs in Def. 5 into the definition of the operator ⊗ in (1).
Next, we consider the Minkowski sum: , which has complexity O(n) with respect to the dimension n.The resulting CPZ is regular if CPZ 1 and CPZ 2 are regular.
Proof The result is obtained by inserting the definition of CPZs in Def. 5 into the definition of the Minkowski sum in (2): where we used the identities ( 11) and (12).
Complexity: The computation of the new constant offset c 1 + c 2 has complexity O(n).Since all other operations required for the construction of the resulting CPZ are concatenations, it holds that the overall complexity is O(n).Now, we provide a closed-form expression for the Cartesian product: , which has complexity O(1).The resulting CPZ is regular if CPZ 1 and CPZ 2 are regular.
Proof The result is obtained by inserting the definition of CPZs in Def. 5 into the definition of the Cartesian product in (3): where we used the identities in (11) and (12).
Complexity: The construction of the resulting CPZ only involves concatenations and therefore has constant complexity O(1).
Before we examine the convex hull, we first derive a closed-form expression for the linear combination since we can reuse this result for the convex hull: , ,( 12) , where we used the identities in (11) and (12).For the transformation in the last line, we substituted λ with an additional factor α p1+p2+1 .Since λ ∈ [ 1, 1] and α p1+p2+1 ∈ [ 1, 1], the substitution does not change the set.
Complexity: The construction of the constant offset c = 0.5(c 1 + c 2 ) requires n additions and n multiplications.Moreover, the construction of the generator matrix The overall complexity is therefore which is O(n 2 ) using (10).
The convex hull can be computed based on the linear combination: where the linear combination comb(CPZ 1 , CPZ 2 ) is calculated using Prop.8 and the scalars p, h, q, and m denote respectively the number of factors, the number of generators, the number of constraint generators, and the number of constraints of the CPZ c, G, E, A, b, R CP Z .The complexity is O(µ 1 + µ 2 ) with respect to the representation sizes µ 1 and µ 2 and O(n 2 ) with respect to the dimension n.The resulting CPZ is regular if CPZ 1 and CPZ 2 are regular.
Proof According to the definition of the convex hull in (4), the definition of the union in (7), and the definition of the linear combination in (8), it holds that The relation in (28) allows us to substitute the union in the definition of the convex hull in (4) with the linear combination.This yields a resulting CPZ with fewer factors compared to using the union according to Theorem 1, which is often advantageous: (1 + λ j ) = 1, ∀j ∈ {1, . . ., n + 1} : ∀j ∈ {1, . . ., n + 1} : ,( 12) where we used the identities in (11) and (12).For the transformation in the last line, we substituted the scalars λ j by additional factors α ap+j .Since λ j ∈ [ 1, 1] and α ap+j ∈ [ 1, 1], the substitution does not change the set.
Complexity: The calculation of the linear combination comb(CPZ 1 , CPZ 2 ) using Prop.8 has complexity O(n(h 1 + h 2 )) according to (26).Moreover, the construction of the constant offset a c requires n multiplications and therefore has complexity O(n).Since all other operations that are required are initializations and concatenations which have constant complexity O(1), the overall complexity for the computation of the convex hull is which is O(n 2 ) using (10).
For the convex hull conv(CPZ) = conv(CPZ, CPZ) of a single set CPZ, we can exploit that CPZ ∪CPZ = CPZ holds to obtain a more compact representation.Next, we consider the quadratic map: The compactGen operation is applied to obtain a regular CPZ.The complexity is O(µ 2 w)+ O(µ 2 log(µ)) with respect to the representation size µ and O(n 3 (w + log(n))) with respect to the dimension n.
Proof The result is obtained by inserting the definition of CPZs in Def. 5 into the definition of the quadratic map in (5), which yields sq(Q, CPZ) .
Note that only the generator matrix, but not the exponent matrix, is different for each dimension x (i) .
Complexity: The construction of the constant offset c has complexity O(wn 2 ) and the construction of the matrices G 1 and G 2 has complexity O(n 2 hw).Moreover, the construction of the matrices E j has complexity O(h 2 p), and the construction of the matrices G j has complexity O(n 2 hw) + O(nh 2 w) if the results for Q i G are stored and reused.The resulting CPZ has dimension n = w and consists of h = h 2 + h generators.Consequently, subsequent application of the compactGen operation has complexity O(ph log(h) + nh) = O(p(h 2 + h) log(h 2 + h) + w(h 2 + h)) according to Prop. 1.The resulting overall complexity is which is O(n 3 (w + log(n))) using (10).
The extension to cubic or higher-order maps of sets as well as the extension to mixed quadratic maps involving two different CPZs are straightforward and therefore omitted.We continue with the intersection: , which has complexity O (µ 1 + µ 2 ) 2 log(µ 1 + µ 2 ) with respect to the representation sizes µ 1 and µ 2 and complexity O(n 2 log(n)) with respect to the dimension n.The compactCon operation is applied to obtain a regular CPZ.
Proof The outline of the proof is inspired by [29,Prop. 1].We compute the intersection by restricting the factors α k of CPZ 1 to values that belong to points that are located inside CPZ 2 , which is identical to adding the equality constraint to CPZ 1 : where we used the identity in (12).Complexity: Computation of c 2 − c 1 has complexity O(n).The resulting CPZ has p = p 1 + p 2 factors, q = q 1 + q 2 + h 1 + h 2 constraint generators, and m = m 1 + m 2 + n constraints.Since the subsequent application of the compactCon operation has complexity O(pq log(q) + mq) according to Prop. 2, we therefore obtain an overall complexity of which is O(n 2 log(n)) using (10).
As a last operation, we consider the union: ) with respect to the representation sizes µ 1 and µ 2 and O(n 3 log(n)) with respect to the dimension n.The compactCon operation is applied to obtain a regular CPZ.
Proof The proof is provided in Appendix A.
Complexity: We first consider the assembly of the resulting CPZ.The computation of the vectors 0.5(c 1 + c 2 ) and 0.5(c 1 − c 2 ) requires n additions, n subtractions, and 2n multiplications.Moreover, computation of 0.5 b 1 , 0.5 b 1 and 0.5 b 2 requires 2m 1 + m 2 multiplications.Computation of the matrix A requires 3 multiplications and 2 divisions.Since the construction of the remaining matrices only involves concatenations, the resulting complexity for the construction of the CPZ is Next, we consider the subsequent application of the compactCon operation.The constraint generator matrix A for the resulting CPZ has q = 1 + q + q + q 1 + q 2 = 4 + 2p 1 + 2p 2 + 2p 1 p 2 + q 1 + q 2 columns since A has one column ( q = 1), A has q = 2 + 2p 1 + 2p 2 + 2p 1 p 2 columns, A 1 has q 1 columns, and A 2 has q 2 columns.Moreover, the matrix ), A has one row (m = 1), A 1 has m 1 rows, and A 2 has m 2 rows.The number of factors of the resulting CPZ is p = p 1 + p 2 + 2. Since the complexity of the compactCon operation is O(pq log(q) + mq) according to Prop. 2, subsequent application of compactCon has complexity ≤ µ1+µ2 ≤ µ1+µ2 Combining ( 29) and (30) yields for the overall complexity with respect to the representation sizes µ 1 and µ 2 .Using (10) it furthermore holds that the overall complexity resulting from the combination of ( 29) and ( 30) is identical to O(n 3 log(n)).
Table 2 Growth of the number of generators h, the number of factors p, the number of constraints m, and the number of constraint generators q for basic set operations on n-dimensional CPZs.

Set Operation Generators Factors Constraints Constraint Generators
Linear map

Representation Size Reduction
As shown in Table 2, many operations on CPZs significantly increase the number of factors, generators, constraints, and constraint generators, and consequently also the representation size.For computational reasons, an efficient strategy for representation size reduction is therefore crucial when computing with CPZs.Thus, we now introduce the operations reduce and reduceCon for reducing the number of generators and the number of constraints of a CPZ.For both operations the tightness of the result can be improved by applying rescaling as described in Sec.4.4 in advance.

Order Reduction
Our method for reducing the number of generators is inspired by order reduction for constrained zonotopes [29,Sec. 4.3] and applies order reduction for polynomial zonotopes: Proposition 12 (Order Reduction) Given CPZ = c, G, E, A, b, R CP Z ⊂ R n and a desired order ρ d ≥ 2 n+m n , the operation reduce returns a CPZ with an order smaller than or equal to ρ d that encloses CPZ: where the sets H and K defined as store the indices of non-zero generators.The compact operation as defined in [21, Prop.2] returns a regular polynomial zonotope and the reduce operation for polynomial zonotopes as defined in [21,Prop. 16] reduces the order to ρ + d .The resulting CPZ is regular and the complexity is O(µ 2 ) + O(reduce) with respect to the representation size µ and O(n 2 ) + O(reduce) with respect to the dimension n, where O(reduce) is the complexity of order reduction for zonotopes.
Proof To calculate a reduced-order CPZ, we reduce the order of the corresponding lifted polynomial zonotope as defined in Lemma 1 using the reduce operation for polynomial zonotopes in [21,Prop. 16].Back-transformation of the lifted polynomial zonotope to the original state space yields an over-approximative CPZ, which can be proven using Lemma 1: where we omitted the compact operation since it only changes the representation of the set, but not the set itself.It remains to show that the order of the resulting CPZ is smaller than or equal to the desired order ρ d .According to [21,Prop. 16], we have where h denotes the number of columns of the matrix G and n + = n + m is the dimension of the lifted polynomial zonotope PZ + .Solving (32) for ρ d yields 2 h/n ≤ ρ d so that holds since the number of elements in the sets H and K is at most h according to (31).
Complexity: Let n + = n + m, p + = p, and h + = h + q the dimension, the number of factors, and the number of generators of the lifted polynomial zonotope PZ + .According to [21, Prop.2], the compact operation has complexity O(p + h + log(h + )) = O(p + (h + q) log(h + q)) and the complexity of order reduction of a polynomial zonotope using [21,Prop. 16] , where O(reduce) denotes the complexity of order reduction for zonotopes, which depends on the method that is used.Moreover, construction of the sets H and K has complexity O((h + q)(n + m)) in the worst case.The overall computational complexity is therefore PSfrag replacements which is O(n 2 ) + O(reduce) using (10).
Let us demonstrate the tightness of our order reduction method for CPZs by an example: , which has order ρ = 7.The resulting CPZ after order reduction to the desired order ρ d = 6 using Prop.12 is visualized in Fig. 5, where we used principal component analysis for order reduction of zonotopes [23, Sec.III.A], which is required for order reduction of polynomial zonotopes using [21,Prop. 16].

Constraint Reduction
Next, we present an approach for reducing the number of constraints of a CPZ, which is inspired by constraint reduction for constrained zonotopes [29, Sec.
the operation reduceCon removes the constraint with index r and returns a CPZ that encloses CPZ: where and the sets H, K, and N are defined as The compactGen operation is applied to make the resulting CPZ regular and the complexity is O(µ 2 ) with respect to the representation size µ and O(n 2 log(n)) with respect to the dimension n.
Proof To remove the constraint with index r, we solve the corresponding equation for the term that is multiplied with the constraint generator with index s: For the constraint with index r, we obtain with the substitution from (35 the trivial constraint 0 = 0, which can be removed by restricting the indices of the constraints to the set K as defined in (34).Finally, inserting the substitution in (35) into the definition of a CPZ in Def. 5 yields The set N as defined in (34) only removes all-zero rows from the exponent matrix and the constraint exponent matrix, and therefore does not change the set.
Complexity: The construction of the matrix G requires nq multiplications, the construction of the matrix A requires (m − 1)q multiplications and (m − 1)(q − 1) subtractions, and the construction of the vectors c and b requires n + m − 1 multiplications, n additions, and m − 1 subtractions.Let p = p, h = h + q − 2, and q = q − 1 denote the number of factors, the number of generators, and the number of constraint generators of the resulting CPZ.Subsequent application of the compactGen operation has complexity O(ph log(h) + nh) = O(p(h + q − 2) log(h + q − 2) + n(h + q − 2)) according to Prop. 1, and the construction of the set N in (34) has in the worst case complexity which is O(n 2 log(n)) using (10).
The crucial point for Prop.13 is the selection of the constraint with index r that is removed, as well as the selection of suitable indices s, d that satisfy the conditions in (33) and can therefore be used for reduction.Clearly, we want to select the indices r, s, and d such that the over-approximation resulting from constraint reduction is minimized.Since it is computationally infeasible to determine the optimal indices for reduction, we instead present some heuristics on how to choose good values for r, s, and d.When removing a constraint from a CPZ using Prop.13, there are two sources contributing to the resulting over-approximation: 1. Over-approximation due to lost bounds on factors.2. Over-approximation due to a loss of dependency.
We now explain both sources in detail and provide illustrative examples.

Over-Approximation due to Lost Bounds
The over-approximation due to lost bounds results from the fact that due to the replacement of the term  s) , where p k=1 [ 1, 1] R (k,s) denotes interval multiplication and exponentiation.Consequently, constraint reduction results in an over-approximation if the solved constraint in (35) has feasible values outside the domain Let us demonstrate this by an example: Example 4 We consider the CPZ which is visualized in Fig. 6.We first choose the indices s = 1 and d = 1 that correspond to the term α 1 for constraint reduction.In this case, solving the constraint α 1 +2α 2 +0.5α 3 3 = 0 for α 1 yields α 1 = −2α 2 −0.5α 3  3 .As visible in Fig. 6 (left), the solved constraint has feasible values outside the domain α 1 ∈ [ 1, 1]: Constraint reduction using Prop.13 with the indices s = 1 and d = 1 therefore results in an over-approximation reduceCon(CPZ, 1, 1, 1) ⊃ CPZ (see Fig. 6 (right)).Next, we consider the indices s = 2 and d = 2 that correspond to the term α 2 .Solving the constraint for α 2 yields α 2 = −0.5α 1 − 0.25α 3  3 , so that the set of feasible values for the solved constraint is a subset of the domain α 2 ∈ [ 1, 1]: Consequently, constraint reduction using Prop.13 with the indices s = 2 and d = 2 does not result in an over-approximation, so that reduceCon(CPZ, 1, 2, 2) = CPZ (see Fig. 6 (right)).
Computing the exact bounds for the solved constraint in (37) is in general computationally infeasible.Instead, one can use range bounding to compute over-approximations of the bounds.Another heuristic that we observed to perform well in practice is to first enclose the CPZ with a constrained zonotope using Prop.4, and then select the constraint that is removed as well as the indices that are used for reduction based on the constrained zonotope enclosure.For constrained zonotopes, sophisticated methods for selecting the constraints and indices that result in the least over-approximation are available in [29,Appendix].

Over-Approximation due to Loss of Dependency
The second source contributing to the over-approximation during constraint reduction is the loss of dependency.This loss arises since if we substitute the term  In the above example we could prevent the loss of dependency by substituting α 2 2 with α 2 2 = (−0.5α 1 − 0.25α 3 3 ) 2 , which corresponds to the square of the solved constraint α 2 = −0.5α 1 − 0.25α 3  3 .Similarly, it is possible to relax the condition in (33) to ∀i ∈ {1, . . ., p}, R e (i,s) = E (i,d) and A (r,s) = 0, which allows powers of the selected term with arbitrary exponents e ∈ N.While this relaxation often enables us to reduce constraints with less over-approximation, taking powers of the solved constraint significantly increases the number of generators of the resulting CPZ, and therefore also the computational complexity for the subsequent compactGen operation.

Numerical Example
For the numerical experiments, we implemented CPZs in the MATLAB toolbox CORA [3], which is available at https://cora.in.tum.de.All computations are carried out on a 2.9GHz quad-core i7 processor with 32GB memory.One often occurring task in set-based computing is to calculate the image of a given initial set under a nonlinear function.We demonstrate by a numerical example that with CPZs the image can be computed exactly if the nonlinear function is polynomial.Let us consider the nonlinear function x 2 (1) + 2 x 2 (2) , 0.5 x 2 (1) ≤ x (2) and the polytope P with vertices [ 1 1] T , [0 1] T , and [1 0] T .The task is to compute the image of P under the nonlinear function f (x).First, we convert the polytope P to a CPZ according to Sec. 5.2 based on a conversion to a polynomial zonotope [20], which yields Using the linear map in Prop.5, the quadratic map in Prop.10, the intersection in Prop.11, and the union in Theorem 1, we can compute the image exactly as where the parameter M and Q for the linear and quadratic map are and the CPZs represent the regions 0.5 x 2 (1) ≤ x (2) and 0.5 x 2 (1) > x (2) , respectively.The resulting image is visualized in Fig. 8. Computation of the image takes 0.02 seconds, and the resulting CPZ has p = 12 factors, h = 13 generators, m = 8 constraints, and q = 85 constraint generators, so that the representation size is µ = 1892.

Conclusion
We introduced constrained polynomial zonotopes, a novel non-convex set representation that is closed under linear map, Minkowski sum, Cartesian product, convex hull, intersection, union, and quadratic and higher-order maps.We derived closed-form expressions for all relevant set operations and showed that the computational complexity of all relevant set operations is at most polynomial in the representation size.In addition, we derived closedform expressions for the representation of zonotopes, polytopes, polynomial zonotopes, )012 Taylor models, and ellipsoids as constrained polynomial zonotopes.Moreover, we demonstrated how to enclose constrained polynomial zonotopes by simpler set representations.In combination with our efficient techniques for representation size reduction, constrained polynomial zonotopes are well suited for many algorithms that compute with sets.

Appendices A Proof for the Union
We now provide the proof for the closed-form expression for the union of two CPZs as specified in Theorem 1.As a prerequisite for the proof, we introduce the constrDom operation: Definition 7 Given a constraint defined by the constraint generator matrix A ∈ R m×q , the constraint vector b ∈ R m , and the constraint exponent matrix R ∈ N p×q 0 , constrDom returns the set of values satisfying the constraint: where α = [α 1 . . .α p ] T .Using Def. 7, it is straightforward to see that the following identity holds: where As another prerequisite, we introduce the following lemma: Proof Using the definition of CPZs in Def. 5 and the definition of the union in (7), we obtain which concludes the proof.

A.2 Reformulation as Union of Sets
We now prove that the resulting CPZ c, G, E, A, b, R CP Z from Theorem 1 defines the union of two sets S 1 and S 2 .Using the domain D as defined in (46) and introducing , where p = p 1 + p 2 + 2, q L = q 1 + q 2 + 1, and α = [α 1 . . .The proof that S 2 = CPZ 2 is similar to the proof for S 1 and therefore omitted.

Fig. 2
Fig. 2 Visualization of the polynomial constraint (left), the constrained polynomial zonotope (right, red), and the corresponding constructing polynomial zonotope (right, blue) for CPZ from Example 1.
Chapter 5], sorting the matrix columns has a worst-case complexity of O(p|N | log(|N |)), which is O(ph log(h)) since |N | ≤ h.The identification and removal of identical neighbors requires at most p(h − 1) comparison operations and therefore has worst-case complexity O(p(h − 1)).Moreover, construction of the sets K and N has complexity O(ph) in the worst case.Finally, the construction of the constant offset c and the generator matrix G has complexity O(nh) in the worst case.The overall complexity is therefore O(ph log(h)) + O(p(h − 1)) + O(ph) + O(nh) = O(ph log(h) + nh).

Fig. 5
Fig. 5 Visualization of order reduction using Prop.12 for the CPZ from Example 3. The original set CPZ (right) and the corresponding constraint (left) are depicted in red, and the reduced order CPZ reduce(CPZ, 6) (right) is depicted in blue.

p
k=1 α R (k,s)k with the solved constraint in (35), we lose the ability to enforce the bounds p k=1

p
k=1 α R (k,s) k with the solved constraint in (35), the dependency between the factors α k in p k=1 α R (k,s) k and the factors α k in other terms gets lost.Let us demonstrate this by an example: Example 5 We consider the CPZ

Fig. 7
Fig. 7 Visualization of constraint reduction using Prop.13 for CPZ from Example 5 (red, right), where the corresponding constraint is visualized on the left.Due to the loss of dependency, constraint reduction reduceCon(CPZ, 1, 2, 2) (blue) results in an over-approximation.

Fig. 8
Fig. 8 Visualization of the polytope P (left, red), the boundary of the region 0.5 x 2 1 ≤ x 2 in (38) (left, blue), and the computed image of P under the nonlinear function f (x) in (38) (right).
The outline of the proof is as follows: We first show in Sec.A.1 that the constraints of the resulting CPZ c, G, E, A, b, R CP Z from Theorem 1 restrict the values for the factors α k to the domain D = D 1 ∪ D 2 corresponding to the union of two domains D 1 and D 2 .Afterward, in Sec.A.2, we apply Lemma 2 to express the resulting CPZ from Theorem 1 as c, G, E, A, b, R CP Z = S 1 ∪ S 2 , where S 1 , S 2 represent the sets corresponding to the domains D 1 , D 2 , respectively.Finally, we show in Sec.A.3 that S 1 = CPZ 1 and S 2 = CPZ 2 holds, which concludes the proof.A.1 Domain Defined by the Constraints First, we show that the constraints of the resulting CPZ from Theorem 1 restrict the values for the factors α k to the domain D = D 1 ∪ D 2 .The matrices A, R and the vector b in Theorem 1 define the constraint α 1 α 2 = 1.(40) The only solutions for (40) within the domain α 1 ∈ [ 1, 1], α 2 ∈ [ 1, 1] are the two points α 1 = 1, α 2 = 1 and α 1 = 1, α 2 = 1.The constraint (40) therefore restricts the values for the factors α k to the domain