Computing Convex Hulls of Trajectories

We study the convex hulls of trajectories of polynomial dynamical systems. Such trajectories include real algebraic curves. The boundaries of the resulting convex bodies are stratified into families of faces. We present numerical algorithms for identifying these patches. An implementation based on the software Bensolve Tools is given. This furnishes a key step in computing attainable regions of chemical reaction networks.


Introduction
Dynamics and convexity are ubiquitous in the mathematical sciences, and they come together in applied questions in numerous ways.We explore an interaction of dynamics and convexity that is motivated by reaction systems in chemical engineering [10,18].Consider an autonomous system of ordinary differential equations where x : R → R n is an unknown function, and φ : R n → R n is a given polynomial map.By the Picard-Lindelöf Theorem, each initial value problem for (1) has a unique solution on a local interval.Although φ is assumed to be polynomial, for most of the techniques we develop it suffices for φ to be locally Lipschitz continuous.Any starting point x(0) in R n gives rise to a unique trajectory C := {x(t) | t ∈ [0, a) for a > 0}.This curve may or may not converge to a stationary point, and its dynamics can be chaotic.We consider the case where the trajectory C is bounded.If it is not bounded, we restrict time t to a finite interval.
Our object of study is the convex hull in R n of the trajectory starting from t = 0: convtraj x(0) := conv(C). ( We call this set the convex trajectory of the point x(0).By definition, it is the smallest convex set containing the trajectory.The convex trajectory need not be closed.In that case, we usually replace convtraj x(0) by its topological closure, so that the convex trajectory becomes a convex body, that is, a compact convex set with nonempty interior.We ensure full-dimensionality of the convex trajectory by restricting to the space in which this curve lies.This article introduces numerical methods that solve the following problems for the dynamical system (1) with respect to a given starting point y = x(0) in R n .
(i) Compute convtraj(y).This is a convex body in R n .The output should be in a format that represents the boundary as accurately as possible.
(ii) Decide whether the convex trajectory convtraj(y) is forward closed.This happens if and only if the vector field φ(z) at every boundary point z points inwards, for every supporting hyperplane.
A subset S ⊂ R n is forward closed if every trajectory of the dynamical system (1) that starts in S remains in S. In the literature, forward closed sets are also referred to as forward invariant sets.If convtraj(y) is forward closed then it equals the attainable region of the point y for (1).The attainable region is the smallest subset of R n that contains y and is both convex and forward closed.We refer to [10,15,18] for motivation and many details.
This article enlarges the repertoire of convex algebraic geometry, a field that studies convex semialgebraic sets, and their role in polynomial optimization [3].Indeed, every algebraic curve is locally the trajectory of a polynomial dynamical system.Hence, our results apply to convex hulls of algebraic curves [19].While Problem (ii) is specific to dynamical systems, Problem (i) makes sense for any smooth curve that can be approximated sufficiently well.In particular, this is the case for parametric curves t → x(t).See Figures 1, 2 and 4.
Approximations by polygonal curves are crucial for our approach.The boundary of the convex hull of a generic curve has a distinct structure that is characterized by families of polyhedral faces spanned by curve points.Our algorithms account for this by constructing polyhedral approximations of convex trajectories, using points that are sampled from the curve of interest.We will discuss conditions on C and the sampling that assures the significance of our methods.More specifically, our results establish a connection between the facial structure of the convex bodies we study and the facets of the approximating polytopes.
Our presentation is organized as follows.Section 2 treats the planar case (n = 2).Example 2.3 illustrates our solution for the Hamiltonian system associated with the Trott curve.Section 3 develops a general geometric theory for taking limits of convex polytopes.This will be applied to identify the desired convex body (2) from a sequence of inner approximations by polytopes.We compute these polytopes using the Matlab/Octave package Bensolve Tools [5].The optimization methodology that underlies this approach is explained in Section 4. Section 5 describes a stratification of a piecewise smooth convex body of dimension n.The strata are (n − k − 1)-dimensional families of k-dimensional faces, called patches, for various k.Section 6 concerns dynamical systems (1) whose trajectories are algebraic or trigonometric curves.This includes linear dynamical systems.Table 1 offers data on their patches for n = 3. Section 7 presents our solution to Problem (ii).Algorithm 7.1 decomposes each patch of the convex trajectory (2) into two subsets, depending whether the vector field φ(z) points inward or outward.If the outward set is always empty then (2) is forward closed.In Section 8 we focus on chemical reaction networks with mass action kinetics [12,14].These motivated our study.We characterize planar algebraic curves that are trajectories of chemical reaction networks, we study the van de Vusse system [10,18], and we exhibit a toric dynamical system [6] with convtraj(y) not forward closed.This resolves [15, Conjecture 4.1].

Planar Scenarios
In this section we study the convex trajectories of dynamical systems in the plane (n = 2).Our input is a pair of polynomials φ = (φ 1 , φ 2 ) along with a starting point y = (y 1 , y 2 ) in R 2 .The resulting trajectory of ( 1) is a smooth plane curve C that is parametrized by time t.We write C = conv(C) for the associated convex trajectory.This is a convex region in R 2 .
The boundary ∂C of the convex trajectory C consists of arcs on the curve C and of edges that connect them.Each edge of C is a line segment between two points on C.These are either points of tangency or endpoints of the curve.This partitions ∂C into patches, described in general in Section 5. Here, k-patches are edges (for k = 1) and arcs (for k = 0).
If we had an exact algebraic representation of the curve C then we could use symbolic methods to compute its bitangents and derive from this a description of ∂C.For instance, if C is an algebraic curve of degree four, as in Example 2.3, then it has 28 bitangent lines (over C) which can be computed using Gröbner bases.In [2] and [9], methods for reducing the number of bitangent line computations are introduced for convex hull computations in the planar case.But, such algebraic representations are not available when we study dynamical systems.Each trajectory is an analytic curve t → x(t).This parametrization is given indirectly, namely by the differential equation (1) it satisfies.Furthermore, our aim is not only to obtain a representation for the convex hull of a curve but also to characterize the structure of its boundary.We are not aware of any reports on such studies.Output a list of curve points that are endpoints of those edges of A, that belong to G i .This represents the ith arc of ∂C.6 end 7 Edges H j of A that correspond to isolated nodes of G represent edges of C.
We now assume that a polygonal approximation is given for the curve C. Our input is a finite list A of points x(t i ) on C. In our computations we solve the differential equation (1) numerically using the versatile solver ode45 in Matlab.This generates the set A of sample points which we assume to be reliably accurate.Using ode45 also allows us to control the quality of the approximation.To ensure a certain precision of the approximation one could employ better methods for integrating dynamical systems like the one presented in [13].
As a first step we address Problem (i) in the Introduction.Algorithm 2.1 computes a representation of the boundary ∂C of the convex trajectory C from a polygonal approximation A of the trajectory C. A key idea is the identification of long edges in A := conv (A).In Section 5 we generalize to curves in R n .Algorithm 2.1 is a special case of Algorithm 5.4.
We next solve Problem (ii) from the Introduction.For each point z on an arc of C, the vector φ(z) is tangent to the curve, so there is nothing to be checked at the arcs.To decide whether C is forward closed with respect to the dynamics (1), we must examine the edges of C. Suppose that conv x(t i ), x(t j ) is an edge of C, and consider its relative interior points where 0 < λ < 1. ( The following 2 × 2 determinant is a polynomial in the parameter λ: We compute all real zeros of the polynomial f (λ) in the open interval (0, 1).The zeros partition the edge of C into segments where φ(z) points either inward or outward, relative to the convex region C.If there are no zeros then the entire edge of C is either inward pointing or outward pointing.In this manner we partition ∂C, and thereby solve (ii).In order to compute the attainable region, we start new trajectories from a sample of outward points.We tested our method on trajectories that are algebraic curves.Let C be the curve in R 2 defined by a polynomial equation h(x, y) = 0.The associated Hamiltonian system equals ẋ = ∂h ∂y (x, y) and ẏ = − ∂h ∂x (x, y).
Henceforth, we only consider Hamiltonian systems associated with polynomials.Thus h is a polynomial in x and y.At any point that is not a critical point of h, the right hand side of ( 5) is orthogonal to the gradient vector of h.This means that the vector field is tangent to the level curves h(x, y) = c, where c ranges over R. From this we infer the following well-known result.
We close this section by illustrating Hamiltonian systems and our solutions to Problems (i) and (ii) for n = 2, for a curve that is familiar in computational algebraic geometry.
The vector field for the Hamiltonian system ( 5) is shown in Figure 1, along with two of its trajectories.Fix any point (u, v) in R 2 and set c = h(u, v).Then the trajectory of (5) that starts at (u, v) travels on the quartic curve defined by h(x, y) = c.Consider the starting point (0, −1), which lies on the original Trott curve h(x, y) = 0. Its trajectory is one of the four ovals, namely the oval at the bottom that is red in Figure 2 (left) and blue in Figure 1.This convex trajectory is not forward closed.This can be seen by taking any starting point (x, b) where 0 < x < 0.239173943.The resulting trajectory is a convex curve that lies above the original oval.It is shown in red in Figure 1.The set of all trajectories as x ranges from 0 to a sweeps out the attainable region for (0, −1).This region is a semialgebraic set.Its boundary consists of parts of two ellipses.Their union is the zero set of h(x, y) − 1053 638 .The containment of the Trott curve in the ellipse is shown on the left in Figure 2. The attainable region of the lower oval in the Trott curve is the convex set of the right in Figure 2.

Limiting Faces in Polyhedral Approximations
We now study the limiting behavior of the facets of a sequence of polytopes.Each polytope is the convex hull of points sampled on a smooth, compact curve C in R n whose convex hull C = conv(C) has dimension n.Our aim is to derive information about faces of C from the facets of its polyhedral approximations.Even though there are lots of profound results on polytopal approximation of convex bodies, see e.g.[4] for a survey, to the best of our knowledge there is no study that gives results similar to the ones presented here.In addition to the approximation of C in terms of Hausdorff distance, in the limit we also approach the special facial structure of the convex hull C of the curve.Our approximation method, which is based on points sampled on that curve, accounts for this.We will elaborate on the structure of ∂C in Section 5.Here we establish a theoretical basis for our algorithms.
We wish to compute the boundary of C using a sequence of inner approximations by convex polytopes, each obtained as the convex hull of a path that approximates C. The Hausdorff distance of two compact sets B 1 and B 2 in R n is defined as A point x of a compact convex set B is extremal if it is not a proper convex combination of elements of B, that is, {x} is a zero-dimensional face of B. Extremal points of a polytope are called vertices.Even if {B ν } ν∈N is a sequence of polytopes that Hausdorff converges to a polytope B, the limit of a convergent sequence of vertices x ν of B ν is not necessarily a vertex of B. For instance, consider B ν = conv({(0, 0), (1, 1 ν ), (2, 0)}).However, a converse holds.Lemma 3.1.Let {B ν } ν∈N → B be a Hausdorff convergent sequence of compact convex sets in R n .For every extremal point x of B there exist extremal points x ν of B ν converging to x.
Proof.Let x be an extremal point of the limit body B. By Hausdorff convergence, there exists a sequence {x ν } ν∈N with x ν ∈ B ν that converges to x.By Carathéodory's Theorem, each x ν is a convex combination of at most n + 1 extremal points v ν 0 , v ν 1 , . . ., v ν n of B ν .Fix some ε > 0. Assume there is an infinite subset N of N such that By compactness, there is a subsequence N of N such that, for each i, the sequence {v ν i } ν∈N converges to some v j .Hence x is a proper convex combination of v 0 , v 1 , . . ., v n ∈ B. This contradicts x being extremal in B. Therefore, for every ε > 0 there exists ν 0 ∈ N such that We consider a sequence {A ε } ε 0 of ε-approximations, where ε 0 stands for a decreasing sequence {ε ν } ν∈N of positive real numbers ε ν .The polytopes A ε = conv(A ε ) can be described by their facets.Our goal is to study convergent sequences of facets F ε of A ε in order to get information about the facial structure of the convex hull C of the curve C. Proposition 3.2.Let {F ε } ε 0 be a Hausdorff convergent sequence of proper faces F ε of the polytopes A ε .Then its limit F is contained in an exposed face of C.
Proof.We write the face F ε of the polytope A ε in the form Proof.A proper face F of C belongs to some hyperplane.By (H3), the set C ∩ F is finite.Since F is a face of C, an extremal point of F is also an extremal point of C. All extremal points of C belong to C, since they cannot be expressed as a proper convex combination of curve points.Thus, F is a polytope as it has only finitely many extremal points.Let F be a proper face of C. We seek a sequence F ε of facets of A ε Hausdorff converging to F .In general, such a sequence does not exist, even under the assumptions (H1), (H2), (H3).We need to additionally require the face F to be uniquely exposed, that is, there is a unique halfspace H + with C ⊂ H + and F = C ∩ −H + .For an example see Figure 3 (right).Theorem 3.5.Assume (H1) and let F be a simplex which is a uniquely exposed face of C. Then F is the Hausdorff limit of a sequence {F ε } ε 0 of facets of A ε .
Proof.Let v 0 , . . ., v k be the vertices of F .Without loss of generality, , where γ = 0 and h ∈ R n with h = 1 is unique.We claim that for every δ > 0 there exists γ > 0 such that We prove this by contradiction.Suppose there exists δ > 0 such that for all γ > 0 there exists y ∈ C\H + γ with y − v i ≥ δ for all vertices v i of F .Thus we can construct a sequence of curve points approaching −H + = −H + 0 but maintain a distance of at least δ from each v i .By compactness of C, this sequence can be assumed to converge to some z ∈ −H + ∩ C ⊂ F .Since the curve point z belongs to the boundary of C, assumption (H1) ensures that z is an extremal point of C. Hence, z is a vertex of F different from v 0 , . . ., v k , a contradiction.
For ε > 0, we consider the linear program We claim that the following holds for sufficiently small ε > 0: Assume the contrary.Then, for all ε > 0, span h and A ε can be separated weakly by a hyperplane γ} weakly separates span h and C. Since 0 is contained in both C and span h, we have γ = 0. Since 0 is a relative interior point of F , we must have F ⊂ H. Hence, F is exposed with respect to a halfspace H+ := {x ∈ R n | hT x ≥ γ} corresponding to H. Since H + = H+ , this contradicts the assumption that F is uniquely exposed with respect to H + .
Let µ ε be the optimal value of ( 7).We have {µ ε } ε 0 = 0. We will use linear programming duality to show that, for ε > 0 sufficiently small, µ ε h belongs to a facet of the form Let (y ε , η ε ) denote an optimal solution of (9).By duality, µ ε = η ε .We conclude that the set F ε is a face of the polytope A ε .
To see that F ε is a facet of A ε , we replace ( 7) and ( 9) by the pair of dual problems Using complementary slackness, we conclude from ( 8) that (y, η) = (0, 0) is the unique optimal solution of (11).Hence the set of optimal solutions of ( 9) is bounded, and we can choose (y ε , η ε ) to be a vertex.At least n linearly independent inequalities in (9) hold with equality at (y, η) = (y ε , η ε ).These correspond to n affinely independent points in A ε , all belonging to the hyperplane 6), we conclude that, for sufficiently small ε > 0, the point µ ε h ∈ F ε (which approaches the mean of the vertices of F ) can be represented only by elements y ∈ A ε with y − v < δ for some vertex v of F .Since F is a simplex, each vertex v of F is used in this representation.Hence, for each vertex v of F there exists a vertex y of F ε with y − v < δ.
We claim that, if ε > 0 is sufficiently small then for every vertex y of F ε there exists a vertex v of F with y − v < δ.Assume the contrary.Then, by (6), for any small ε > 0, there is a vertex y ε of F ε with y ε ∈ C ∩ H + γ .By compactness of C, we may assume that {y ε } ε 0 converges to some ȳ ∈ C ∩ H + γ .By Proposition 3.2, conv (F ∪ {ȳ}) belongs to an exposed face of C. Since ȳ ∈ −H + , this contradicts the assumption that F is uniquely exposed.We conclude that for every δ > 0 we find ε 0 > 0 such that d(F ε , F ) < δ for all 0 < ε ≤ ε 0 .Hence, the simplex face F of C is the Hausdorff limit of the sequence {F ε } ε 0 , as desired.
Remark 3.6.In the proof of Theorem 3.5, we construct a sequence {F ε } ε 0 of facets of A ε whose vertices converge to the vertices of F .This is stronger than Hausdorff convergence.Remark 3.7.In our proofs we had assumed, for simplicity, that the curve C is smooth and that the sampled points lie exactly on C. Our results should extend to nonsmooth curves and to points that are sampled around C with a given accuracy.The investigation of polyhedral approximations of C = conv(C) under such weaker assumptions is left for future work.

Convex Hulls in Bensolve
The key step in our solution to Problem (i) is the computation of the convex hull of an ε-approximation of a curve C.There are many methods and implementations for convex hulls.For this paper, the software Bensolve Tools [5] was used.It is based on Benson's algorithm; see e.g.[8].One reason for that choice is the output sensitivity of the underlying method.This means that the runtime is mainly dependent on the number of facets and vertices of the polytope that is the convex hull and its dimension.In particular, the number of sampled points in the interior of C A only marginally influences the computation time.Another advantage of Benson's algorithm is the possibility to set the parameter ε in Algorithm 4.4.This feature enables the approximative representation of highly complex convex hulls in a reasonable amount of time.In addition, the process can be aborted at any point while still providing an outer approximation.For small values of ε, we obtain exact solutions, up to numerical inaccuracy of the vertex enumeration routine.
We next discuss this software, its underlying methodology, and how we apply it.Bensolve [17] is a solver for multiple objective linear programs (MOLP).In Bensolve Tools it is utilized to perform many polyhedral calculus operations, among them convex hull.The key insight behind this is that multiple objective linear programming is equivalent to polyhedral projection [16].Convex hull computation is a special case of polyhedral projection.This follows from standard facts in Ziegler's textbook [22,Chapter 1].We state it as follows: where V ∈ R n×k is the matrix with set of columns V and e = (1, 1, . . ., 1) T .Hence, the convex hull of V is a projection into R n of the polytope To understand the computation of conv(V) from Q, let us turn to an arbitrary polyhedral projection problem.By Fourier-Motzkin Elimination, every linear projection of a polyhedron is a polyhedron.This leads to the concept of a P-representation of a polyhedron.Let M ∈ R n×k , B ∈ R m×k , a ∈ R m be given.The triple (M, B, a) represents the polyhedron In what follows, we restrict to polytopes (bounded polyhedra).Given a Prepresentation (12) of a polytope, the polyhedral projection problem is to compute an irredundant V-representation, i.e. a representation as convex hull of finitely many points, and an irredundant H-representation, i.e. a representation by finitely many linear inequalities (cf.[22]).
Given a triple (M, B, a) as above, the associated multiple objective linear program is min M x s.t.Bx ≥ a. (MOLP) The upper image of the program (MOLP) is the polyhedron A solution of (MOLP) consists of irredundant V-and H-representations of P.This concept of solution can be used to address the polyhedral projection problem: yields an irredundant V-and H-representation of the P-represented polyhedron (12).
The upper image of the MOLP in ( 14) is the polyhedron Corollary 4.3.The polytope P = {M x | Bx ≥ a} is obtained from P by setting An irredundant V-representation of P derives from the set of vertices of P by deleting their last components.An H-representation of P gives an H-representation of P by adding the equation z = −e T y.
Bensolve computes V-and H-representations of the upper image (15) using Algorithm 4.4.This is a version of Benson's algorithm.It applies to upper images satisfying P ⊆ y + R n ≥0 for some y ∈ R n .This version suffices for handling projections of polytopes including the representation of the convex hull of finitely many points.Since the algorithm is numerical, we work with a prescribed tolerance ε > 0. The output is an ε-approximation to the upper hull P, i.e. it is a polyhedron O that is ε-close to P in the sense that εe + O ⊆ P ⊆ O.One starts with an initial outer polyhedral approximation of P. Both an Hrepresentation and a V-representation are stored.Until the tolerance ε is reached, each iteration adds a linear inequality to refine the outer approximation of P.An iteration step starts by choosing a vertex v of the current polyhedron.The Vrepresentation is updated after adding an inequality.From v one moves in direction e = (1, . . ., 1) T to the boundary point y = v+t * e of P. To this end, a linear program has to be solved.The solution of the dual linear program yields the desired linear inequality which cuts off v and holds with equality in y.Algorithm 4.4 terminates and computes both V-and H-representation of an ε-approximation of P. For computations in this paper we use the dual Benson algorithm [8].It is dual to Benson's algorithm and provides an inner approximation for the upper hull P.
We employ Bensolve for obtaining a polyhedral approximation of the convex hull of a smooth curve C in R n .This is done by computing the convex hull of a sufficiently large finite subset A of C. The output gives both an irredundant Hand V-representation of an inner ε-approximation C A of conv(A), and this is our approximation to conv(C).All facets and all vertices of C A are known after such a computation.The output also contains the incidence matrix I A for facets and vertices of C A and the adjacency matrix A A for vertices of C A .
Example 4.5.Let C be the trigonometric space curve parametrically given by θ → (cos(θ), sin(2θ), cos(3θ)) . ( Its convex hull C = conv(C) is shown in [19, Figure 1].We select the sample points Using Bensolve, as described above, we can compute irredundant V-and Hrepresentations of the inner approximation C A of the polytope conv (A) for various values of N and with specified accuracy ε.For instance, let N = 100 and ε = 10

Detection of Patches
Let C be a full-dimensional compact convex body in R n .The boundary of ∂C is an (n − 1)-dimensional set whose subset ∂C sm of smooth points is dense.We shall stratify ∂C sm into finitely many manifolds we call patches.Each patch is an (n − k − 1)-dimensional family of k-faces of C. For a typical convex body of dimension n = 3, the boundary is comprised of surfaces of extreme points (k = 0), curves of edges (k = 1), and finitely many facets (k = 2).For the general definition, we use the concept of the normal cycle of a convex body.Let S n−1 denote the unit (n − 1)-sphere.Following [11, eqn (10)], the normal cycle of C equals If ∂C is smooth then N (C) is a Legendrian submanifold of dimension n − 1.If C is not smooth then we can approximate C by nearby smooth convex bodies C ε , for ε > 0. By [11, Theorem 3.1], the normal cycle N (C) is the Hausdorff limit of the manifolds N (C ε ) for ε → 0. The normal cycle N (C) is pure (n − 1)-dimensional, and its smooth points are dense.
There are several other ways of defining the normal cycle.The one we like best uses the dual convex body C ∨ .Assuming that the origin is in the interior of C, The normal cycle comes naturally with two surjective maps Let E ⊆ ∂C ∨ be the set of exposed points of C ∨ .We have v ∈ E if and only if there the fibers of π 2 vary continuously in the Hausdorff metric, and ψ is maximal with these properties.We say that ψ is a k-patch if dim(π 2 (ψ)) = n − k − 1.This means that π 2 (ψ) is an (n − k − 1)-dimensional manifold of exposed points of C ∨ , and these exposed points support continuously varying k-faces of C. Remark 5.1.If the trajectory C is algebraic then its convex hull C is semialgebraic.Also the normal cycle N (C) and all its patches ψ are semialgebraic.This follows from Tarski's theorem on quantifier elimination, and we find that the number of patches of C is finite.We believe that finiteness holds more generally for compact trajectories.But we do not yet know the precise statement.Real analytic geometry is much more delicate than real algebraic geometry.For instance, the family of semianalytic sets is not closed under projection.We refer to [1] for a recent account.The concept of C-semianalytic sets, introduced in [1] and named after Cartan, might be appropriate for our setting.One can hope that the convex trajectories and their patches are C-semianalytic when φ in (1) is polynomial.
We now study the following computational problem.The input is a smooth curve C in R n , typically arising as trajectory of a dynamical system (1).We seek the boundary of the convex trajectory C=conv(C).The output is the list of all patches.
Example 5.2 (n = 3).The convex body in Examples 4.5 and 6.2 has six patches.It has two 2-patches, namely the two triangles.It has two irreducible edge surfaces (cf.[19]).These have degrees 3 and 16.Each contributes two 1-patches to ∂C.All six patches are visible from the edges and triangles of the polytope in Figure 4.
Example 5.3 (n = 4).If C is the convex hull of a curve C in R 4 then its 1-patches are surfaces of edges, its 2-patches are curves of 2-faces, and its 3-patches are the facets of C.
Our goal is to identify all patches of C = conv(C) from ε-approximations, using the results in Section 3. We assume that C is a simplicial curve, by which we mean that it satisfies the hypotheses (H1), (H2) and (H3), and the number of patches of C is finite (cf.Remark 5.1).
The special case n = 2 is solved by Algorithm 2.1.Algorithm 5.4 computes the patches for n ≥ 3. We implemented this algorithm for n = 3 and n = 4.A detailed theoretical analysis of this algorithm is left for future work.The task is to identify the precise conditions under which the output detects the true patches when applied to ε-approximations with ε → 0.
In the next section we report on some experiments with our methods.The code for dimensions 2, 3 and 4 is made available at http : //tools.bensolve.org/trajectories.
We discuss the steps in Algorithm 5.4.Step 1 is executed using Bensolve as discussed in Section 4. Each facet H comes with its unit normal vector v(H).
Step 2 reflects the conditions in the definition of patches.For instance, the criterion d(H 1 , H 2 ) ≤ δ, stating that H 1 and H 2 are close in Hausdorff distance, reflects the continuous variation of k-faces.The exposed points in π 2 (ψ) ⊂ E are represented by the vectors v(H i ).The requirement that they are δ-close along the edges of G is our discrete version of the smoothness of ψ.In step 3 we identify the connected components of G, and these represent the patches of C.
The inner loop in steps 5-7 reflects our results in Section 3. By Theorem 3.5, every k-face of C is approximated by a facet H ∈ G. Here, for each vertex of the k-face, the algorithm chooses a nearby vertex u i of H.In the loop between steps 4 and 8, it can happen that a δ-proximity cluster corresponds to more than one vertex of the k-face.This happens for k-faces with an edge shorter than δ.For that reason, we take the maximum in step 9.
Step 10 is very important and requires some explanation.In a connected component of G, some facets will be misclassified: they represent k-faces of C, but step 5 identifies less than k + 1 proximity clusters at the fixed tolerance level δ.For such facets, step 10 adds additional points u i from an existing cluster to get up to the correct value of k for that patch.Associate the tuple (u 0 , . . ., u k ; v) with that node of G i .
Adjust all tuples with smaller indices found in step 7 to that common value of k. end Output (# 1 , . . ., # n−1 ), where # k is the number of graphs G i representing k-patches.
Example 5.5 (n = 4).We illustrate the output of Algorithm 5.4 when C is a random trigonometric curve of degree six in R 4 , as in Section 6, and A ⊂ C is a finite approximation.Figure 5 shows the graph G.Each node in G is a face of the polytope conv(A).We find # 3 = 0.There are # 1 = 3 patches for k = 1, represented by the three connected components of G on the right in Figure 5.These three connected graphs encode surfaces worth of edges.The number of patches for k = 2 is # 2 = 2.These two components of G are shown on the left in Figure 5.
Each node represents a triangle face of conv(C).So, the picture on the left shows two curves of triangle faces in the boundary of our 4-dimensional convex body.

Algebraic and Trigonometric Curves
In what follows we note that every algebraic curve can be realized locally as the trajectory of an autonomous polynomial dynamical system.This generalizes the Hamiltonian systems (5) seen in Section 2. Hence, the computation of the convex hull of a real algebraic curve in R n , discussed in e.g.[19], is a special case of the problem we addressed in Sections 3-5.
Let C be an algebraic curve in R n and let z be a regular point on C. We construct an appropriate vector field φ(x) on R n as follows.Let f 1 , f 2 , . . ., f n−1 be polynomials in x = (x 1 , x 2 , . . ., x n ) that cut out the curve C locally near its point z.Let J denote their Jacobian matrix.Thus J is the (n − 1) × n matrix whose entry in row i and column j is the partial derivative ∂f i /∂x j .Let J i be (−1) i+1 times the determinant of the submatrix of J obtained by deleting the ith column.Fix the vector of polynomials φ = (J 1 , J 2 , . . ., J n ) T .Locally at z, the kernel of J is the line spanned by the vector φ.This follows from Cramer's rule, and it implies that φ(z) is a tangent vector to the curve C at its point z.We are interested in the dynamics of the system ẋ = φ(x) when the starting point is z ∈ C. Proposition 6.1.The trajectory of the dynamical system ẋ = φ(x) that starts at a point z on the algebraic curve C remains on the curve C. It either cycles around one nonsingular oval of C, or it diverges towards infinity in R n , or it converges to a singular point of C. Proof.The proof is analogous to Corollary 2.2 which dealt with the case n = 2. Example 6.2 (n = 3).Let C be the trigonometric curve (16) in Example 4.5.By [19, §1], this is an algebraic curve, namely it is the zero set of the two polynomials With these two polynomials we associate the dynamical system ẋ = −2y and ẏ = 12x 3 − 5x + z and ż = −24x 2 y + 6y.
Suppose we start this at a point on the curve C, such as (1, 0, 1).The trajectory travels on C and it stops at the singular point (0, 0, 0).To get conv(C), we compute the convex hull of two trajectories obtained by using two different starting points on the curve given by (19).The convex body has six patches, as shown in [19, Figure 1] and in our Figures 4 and 6.
Using the methods in Section 7, we analyzed the vector field on these facets and patches, and we found points with both inward and outward pointing directions on each of them.Figure 6 shows a point in a triangle facet with outward pointing direction.The resulting trajectory is also depicted.This solves Problem (ii) from the Introduction for this example.
Trigonometric curves also arise from linear dynamical systems.Here (1) takes the form ẋ = Ax, where A is a real n × n-matrix.We tested our convex hull algorithms on linear systems for n = 3, 4. We sampled matrices A with no real eigenvalues.This ensures that the trajectories are bounded in R n .They can be written in terms of trigonometric functions.It was shown in [15] that every convex trajectory of a linear dynamical system is forward-closed.Thus, computing the convex trajectory is equivalent to computing the attainable region.
Consider the generalized moment curve, whose convex hull was studied in [21, Theorem 1].Let z = (1, 0, 1, 0) and consider the linear dynamical system given by where p and q are relatively prime positive integers.The trajectory is the curve The curve is closed, and we can restrict to 0 ≤ t < 1.The convex hull of the curve is a 4-dimensional convex body.By [21, Theorem 1], there are no 3-dimensional faces.Assuming p, q ≥ 3, there are two 1-patches and two 2-patches.The explicit description in [21] makes this a useful test case.
Example 6.3 (p = 3, q = 4).The segment conv{x(s), x(t)} is an edge if and only if In addition to this surface of edges, there are two curves of 2-faces, namely the triangles conv x t , x t and the squares conv x t , x t + 1 4 , x t + 1 2 , x t + 3 4 for 0 ≤ t < 1 4 .These are all the exposed faces of the convex trajectory.Even though the curve is not simplicial, Algorithm 5.4 works well, and we verified Smilansky's findings using our software.
We experimented with our Bensolve-based code for random trigonometric curves x : [0, 1] → R n .The coordinates of x are trigonometric polynomials of the form We write the coefficients as a pair of n×d matrices A and B together with a column vector C, all filled with real numbers.For general matrices, the resulting curve is an algebraic curve of degree 2d in R n .We computed many examples and recorded the features seen in the boundary.We were most interested in the maximal numbers of facets that were observed.We sampled random data (A, B, C) and computed the convex hull of the resulting curves.For n = 3 we recorded the number of triangles (= 2-patches).The second row in Table 1 shows the maximal number of triangles that was observed for given degree 2d.Each triangle spans a real tritangent plane of the curve.The third row lists the number of complex tritangent planes for this curve, which is a generic space curve of genus 0 of degree 2d.The edge surface of the curve is an irreducible ruled surface that defines the nonlinear part of the boundary of the convex hull.Its degree is listed in the fifth row.This surface is the Zariski closure of any of the  We illustrate our computational results in Table 1 for a curve of degree 2d = 14.
Example 6.4 (2d = 14).We consider the curve defined by the 3 × 7 matrices   We seek to partition the boundary ∂C into two regions.In one region, the vector φ(z) points inward and in the other it points outward, as in Figure 6.For n = 2 this partition is determined by the formulas (3) and (4).In what follows we generalize this method to n ≥ 3. Let ψ be a k-patch of C. The output of Algorithm 5.4 represents ψ by a connected graph G i .Each node of G i is a face F = conv{u 0 , . . ., u k } along with a normal vector v at C. We are interested in the restriction of the boundary above to the patch ψ of interest: We obtain an approximate representation of ∂C employing Algorithm 5.4.Algorithm 7.1 computes a partition of this approximation into inward and outward pointing regions.The simplex ∆ k is the convex hull of the unit vectors in R k .Consider a fixed graph G i in the loop started in step 1.The value of k is constant in step 2. This is ensured by step 10 in Algorithm 5.4.The graph G i represents a k-patch of C. We algorithmically realize the restriction of the hypersurface (21) to the k-faces in that patch in step 4. If k = 1 then this results in a finite partition of a line segment.For k = 2 we obtain a curve in a triangle, and for k = 3 we obtain a surface in a tetrahedron.The latter case happens only for n ≥ 4.
The paradigm for our computations is the algebraic case.Suppose that C is an algebraic curve, for instance obtained from a dynamical system as in Proposition 6.1.In that case, the equation φ(u)•v = 0 is a polynomial in k unknowns λ 1 , . . ., λ k , after setting λ 0 = 1 − k j=1 λ j .To be precise, let v = (v 1 , . . ., v n ) be the normal vector of the k-face in question.In the situation of Proposition 6.1, the polynomial equation we are solving on ∆ k takes the form In some situations, we know the equation f = 0 of the hypersurface π 1 (ψ) in R n .Here f is analytic or polynomial, depending on the instance.With this, we write This formula allows us to solve the equation φ(u)•v = 0 simultaneously on the entire and exact k-patch, and not just on each approximated k-face of ψ individually, as it is done in step 4 of Algorithm 7.1.(19).The degree 16 polynomial g is displayed in [19, §1].On the cubic patches, the equation ∇f 2 •φ = 0 holds identically, so these patches are not partitioned.Hence, every trajectory that starts on a cubic patch remains in that patch.The two degree 16 patches are partitioned by a curve of degree 262, obtained by intersecting the patches with the surface defined by ∇g•φ = 0.The two triangle facets lie in the planes z = ±1.They are partitioned by the lines y = 0 and x = ±1/2.Figure 6 shows the trajectory that starts at a red point in the outward pointing region of the top triangle.
In the next section we apply our methods to partition the boundary of convex trajectories of dynamical systems that arise from chemical reaction networks.Figure 9 shows our partition for a convex trajectory arising in an application.

Chemical Reaction Networks
Our interest in convex trajectories and attainable regions is motivated by dynamical systems for chemical reaction networks.The aim of attainable region theory [18] is to design chemical reactors that are optimal for chemical reactions of interest.This research topic was pioneered by Feinberg and Hildebrand in [10].They argued that optimal reactors are often found at the extreme points of the attainable region and showed that these extreme points are realizable by parallel operations of elementary reactor types.We refer to [15] for a recent study in the setting of convex algebraic geometry [3,19].The notion of protrusions in [10, §2.6] is dual to our notion of patches in Section 5, in the sense that a k-patch on C corresponds to an (n − k − 1)-dimensional protrusion on C ∨ under the self-duality of the normal cycle N (C) in (18).Patches and protrusions are interesting for further research.
To explain the connection to chemical reactions, we work in the setting of mass action kinetics.Let G be a directed graph with m vertices, each labeled by a monomial x ai in n unknowns x = (x 1 , . . ., x n ).These unknowns are the concentrations of n chemical species.The m monomials are the chemical complexes.Each x j = x j (t) is a function of time t.With each edge (i, j) of G, connecting two monomials x ai and x aj , we associate a parameter κ ij , which is the rate constant for that reaction.The associated dynamical system is given by φ This is a row vector of length n, written as a product of three matrices, of formats 1 × m, m × m, and m × n.The middle matrix Λ G (κ) is the Laplacian of the graph G, with entry κ ij for each edge (i, j), all other off-diagonal entries set to zero, and diagonal entries inferred so that the row sums of Λ G (κ) are zero.An important special case of ( 22) are the toric dynamical systems [6].For details see the forthcoming book by Dickenstein and Feliu [7] and Shiu's dissertation [20, §1.3].
A feature of many chemical reaction systems (22) is the existence of conservation relations.These arise when the entries of the polynomial vector φ(x) are linearly dependent over R. If this happens then all trajectories lie in certain lower-dimensional subspaces of R n .In such cases, the ambient dimension n can be reduced.Namely, we always transform our dynamical systems so that each trajectory affinely spans R n .We examined some important classes with n ≤ 4. Of particular interest are chemical reaction networks that admit multistationarity.The smallest such networks were characterized by Joshi and Shiu [14].
Given an arbitrary polynomial dynamical system (1), it is natural to ask whether it arises from some chemical reaction network G.The solution to this inverse problem was given by Hárs and Tóth [12].They showed that φ = (φ 1 , . . ., φ n ) is realized by a graph G as above if and only if each monomial with negative coefficient in φ i is divisible by x i , for i = 1, 2, . . ., n.
We discussed Hamiltonian systems in Section 2. The following result characterizes chemical reaction dynamics in R 2 that is Hamiltonian.It would be interesting to study such reaction networks, along with the higher-dimensional versions arising from Proposition 6.1.Proposition 8.1.Let h(x, y) be a polynomial.The Hamiltonian system (5) can be realized as a mass action system (22) if and only if the coefficients of all powers of y in h(x, y) are nonnegative and those of all pure powers of x are nonpositive, i.e.Proof.This is immediate from Theorem 3.2 in [12].
The convex trajectory we compute as an answer to problem (i) in the Introduction is a first approximation to the attainable region and its representation in [10].In Section 7 we presented an algorithm for partitioning the approximated boundary of the convex trajectory.We next apply that algorithm to two interesting chemical reaction networks.
Example 8.2 (n = 4, m = 5).We revisit the Van de Vusse reaction.This is studied extensively in the chemistry literature (cf.[18,Chapter 6]).The network equals The rate constants κ ij are written over the edges.The mass action system (22) where x i is the concentration of species X i .Explicitly, this is the system (1) with φ(x 1 , x 2 , x 3 , x 4 ) = −x 1 − 20x 2 1 , x 1 − x 2 , x 2 , 10x 2 1 .Computations of the critical reactors of this system are found in [18,Section 5.3].
Fix the starting point y = (1, 0, 0, 0).The trajectory starting at y converges to the steady state y * = (0, 0, 0.1522, 0.4238).The dynamics takes place in R 4 , but the stoichiometry space has dimension 3.In our analysis we use the projection onto the first three coordinates.With this, the trajectory is an arc in a 3-dimensional space, shown in blue in Figure 9.
We computed the convex trajectory starting at y for (23) using Algorithm 5.4, and we then partitioned its boundary using Algorithm 7.1.The result is shown in Figure 9.The convex body has two 1-patches, obtained by joining each of the two endpoints with each point on the curve.One of the patches is entirely green.This means that the vector field is pointing inward on that patch.The other patch is partitioned into a green region and a red region, as shown on the right in Figure 9. Red color indicates that the vector field points outward.In particular, the convex trajectory is strictly contained in the attainable region.The mass action system (22) is called weakly reversible if every connected component of the underlying directed graph is strongly connected in G, i.e. there is The proof is by computation using our algorithms.Here is the counterexample: Example 8.4 (Weakly Reversible System).Consider the following weakly reversible network The three coordinates for (1) are explicitly given by φ 1 = −10x 2 1 x 2 + 10x 2 2 x 3 − 4x 1 x 2 + 2x 1 x 3 , φ 2 = 2x 2 1 x 2 − 6x 2 2 x 3 + 4x 1 x 2 + 2x 1 x 3 , φ 3 = 6x 2 1 x 2 − 6x 2 2 x 3 + 4x 1 x 2 − 2x 1 x 3 .This system has deficiency zero, and it is a toric dynamical system [6].There are no conservation relations.The trajectories are curves that span the ambient space R 3 .Figure 10.Convex trajectory of a weakly reversible system that is not forward closed.
Let y = (4, 4, 2).The convex body C = convtraj(y) was computed using Algorithm 5.4 and is shown in Figure 10.The triangle shown in gray is a 2-patch of C. The vector field given by (φ 1 , φ 2 , φ 3 ) points inward at all points on that triangle facet.We also show the partition of the 1-patches of C, as computed by Algorithm 7.1.One of the patches is partitioned into a green region and a red region.As before, the vector field points outward at each red point.We conclude that the convex trajectory C of y is not forward closed.

Algorithm 2 . 1 .
(Detection of edges and arcs for n = 2) input : A list A of points on a curve C in R 2 ; a threshold value δ > 0 output: The numbers # 0 and # 1 of arcs and edges of C = conv (C)For each i: list of curve points that represent the ith arc of ∂C.List of line segments that represent the edges of C. 1 Compute the vertices V and edges H of A = conv (A).

Figure 1 .
Figure 1.The Hamiltonian vector field defined by the Trott curve and two of its trajectories.

Figure 2 .
Figure 2. A pair of ellipses encloses the Trott curve and bounds the attainable region.
Proposition 3.2 may not be a face of C.This is shown in Figure 3 on the left.The following genericity assumptions on C will ensure that a Hausdorff convergent sequence of proper faces F ε of the polytopes A ε converges to a proper face F of the body C: (H1) Every point on the curve C that is in the boundary of C is an extremal point of C. (H2) Every polytope face of C is a simplex.(H3) Intersecting the curve C with a hyperplane always results in a finite set.We now give a sufficient condition that proper faces of C are polytopes.Proposition 3.3.If C satisfies (H3), then every proper face of C is a polytope.

Figure 3 .
Figure 3.A Hausdorff convergent sequence of facets F ε of A ε need not converge to a face F of C. The face F on the left contains a curve point y ∈ C which is not extremal in C. The endpoint F of the curve C on the right is an exposed face of C but it is not uniquely exposed.There is no sequence of facets F ε of A ε that Hausdorff converges to that face F .

Algorithm 4 . 4 . 5
(Benson's algorithm) input : (MOLP) given by the matrices M , B and vector a; a tolerance ε ≥ 0. output: ε-close V-representation V and H-representation H of P in (13). 1 T ← ∅ 2 Compute the H-representation H of an outer approximation of P having the same recession cone as P. Compute the corresponding V-representation V. 3 while (V \ T ) = ∅ do 4 Choose a vertex v ∈ V \ T .Compute the solution t * of the linear program min {t | v + te ∈ P}. 6 Compute the solution (u * , w * ) of the dual linear program max a T u − v T w | B T u = M T w, e T w = 1, w ≥ 0, u ≥ 0 .7 if t * ≥ ε then 8 Refine H by adding y | (w * ) T y ≥ a T u * to the description.9 Update V by performing vertex enumeration on H. else T ← T ∪ {v} end end

− 9 .
The sample(17) is shown on the left in Figure4.Its convex hull is the polytope C A on the right in Figure 4.It has 70 vertices and 102 facets, so, by Euler's relation, it has 170 edges.Thus the incidence matrix I A is of size 102 × 70 and has 340 nonzero entries.The 70 × 70 adjacency matrix A A also has 340 nonzero entries.The polytope C A in Figure 4 already looks like [19, Figure 1] and Figure 6.The picture of C A reveals the edge surfaces of C and the two triangles in ∂C.The identification of such patches from the Bensolve output is our theme in Section 5.

Figure 4 .
Figure 4.A sample of points (left) from a space curve and its convex hull (right).

Algorithm 5 . 4 .
(Detection of patches for n ≥ 3) input : Finite list A of points on a curve C in R n ; a threshold value δ > 0 output: For each k ≥ 1: the expected number # k of k-patches of C = conv (C) For each i: list of k-polytopes that represent the k-patch G i 1 Compute vertices V, facets H, incidence list I A and adjacency list A A of conv (A). 2 Build a graph G with node set H as follows: two facets H 1 , H 2 form an edge if their unit normals

Figure 5 .
Figure 5. Two 2-patches (left) and three 1-patches (right) in the boundary of a 4-dimensional convex body.It is the convex hull of a trigonometric curve of degree six.The picture shows the graph G, with five connected components G i , found by Algorithm 5.4.

Figure 6 .
Figure 6.The red boundary point shows that the convex trajectory is not forward closed.

5 Figure 7 .
Figure 7.The convex hull of a trigonometric curve of degree 14 in 3-space.The boundary of this convex body consists of triangles and of 1-patches in a ruled surface of degree 286.

, 7 .
along with the vector C = 0.39768 0.42346 0.23797T .The convex hull of this curve has 20 triangle facets.It is shown in Figure7.The planes that define the triangles are tritangent planes.The curve is generic and has 1320 tritangent planes over C. The nonlinear part of the boundary is the edge surface[19].This is an irreducible ruled surface of degree 286.Partitioning the Boundary Problem (i) from the Introduction was addressed in the previous sections.In this section we propose a solution to Problem (ii).Our input now is the output of Algorithm 2.1 or Algorithm 5.4.If n = 3 then we are given all 2-patches (triangles) and all 1-patches (in the edge surface) of the convex hull C = conv(C) of a trajectory C of the dynamical system (1).For instance, the output of Algorithm 5.4 might be Figure8.This is a variant of Figure7, derived from a trigonometric curve C of degree 14.Here C is simplicial, and the boundary ∂C consists of 20 triangles (= 2-patches) and 30 1-patches, shown in different colors in Figure8.

Figure 8 .
Figure 8.The convex hull of a trigonometric curve of degree 14.
For a pair (u, v) in the normal cycle N (C), v is an outward pointing unit normal vector at C. The right hand side of (1) points inward at u ∈ ∂C if φ(u) • v ≤ 0 for all v ∈ π 2 (π −1 1 (u)).It points outward otherwise.Taking into account that the set {u ∈ R n | π 2 (π −1 1 (u)) = {v}} is dense in ∂C, we see that the boundary between inward and outward pointing vectors φ(z) is the image under the projection π 1 : N (C) → ∂C of the set (u, v) ∈ N (C) : φ(u) • v = 0 .This image has dimension n − 2 inside the (n − 1)-dimensional normal cycle N (C).

Algorithm 7 . 1 . 5 end 6 end
(Partitioning the boundary of a convex trajectory) input : The graphs G i representing the patches of a convex trajectory of (1) output: Partition of the boundary into inward and outward pointing regions 1 foreach connected graph G i in the output of Algorithm 5.4 do 2 foreach node ({u 0 , . . ., u k }, v) of the graph G i do 3 Set u = k i=0 λ j u j where λ j are nonnegative unknowns satisfying k j=0 λ j = 1. 4 Compute the (k − 1)-dimensional hypersurface in ∆ k defined by φ(u) • v = 0 and identify inward and outward pointing regions.Example 7.2 (n = 3).We partition the boundary of the convex body in Figure 6.Its edge surface has two irreducible components, of degrees 3 and 16.Each contributes two patches.The cubic is f h(x, y) = xy • a(x, y) − b(x) + c(y), where b and c have nonnegative coefficients.

Figure 9 .
Figure 9. Convex trajectory of the Van de Vusse reaction and the partition of its boundary.
Output the number # 1 of isolated nodes of G and the number # 0 of remaining connected components G i .4 foreach nonsingleton connected component G i do 2 Build a graph G with node set H such that two distinct edges H 1 , H 2 of A form an edge of G if H 1 ∩ H 2 = ∅ and both H 1 and H 2 have length ≤ δ. 3

Table 1 .
Census of random trigonometric curves in 3-space 1-patches.The fourth row shows the maximal number of observed 1-patches.The numbers in the third and fifth row are taken from [19, Corollary 3.1].