Steplength Thresholds for Invariance Preserving of Discretization Methods of Dynamical Systems on a Polyhedron

Steplength thresholds for invariance preserving of three types of discretization methods on a polyhedron are considered. For Taylor approximation type discretization methods we prove that a valid steplength threshold can be obtained by finding the first positive zeros of a finite number of polynomial functions. Further, a simple and efficient algorithm is proposed to numerically compute the steplength threshold. For rational function type discretization methods we derive a valid steplength threshold for invariance preserving, which can be computed by using an analogous algorithm as in the first case. The relationship between the previous two types of discretization methods and the forward Euler method is studied. Finally, we show that, for the forward Euler method, the largest steplength threshold for invariance preserving can be computed by solving a finite number of linear optimization problems.

1. Introduction. Invariant set is an important concept in the theory of dynamical systems and it has a wide range of applications in control. One of the reasons of the interest is due to the fact that invariant sets enable us to estimate the attraction region of a dynamical system. We consider linear continuous dynamical systems in the formẋ and discrete dynamical systems in the form where x k , x(t) ∈ R n are the state variables, A c , A d ∈ R n×n are the coefficient matrices, and t ∈ R and k ∈ N indicate continuous and discrete time steps, respectively. Note that equations (1) and (2) can be treated as autonomous systems or as controlled systems. In the latter case, the coefficient matrix A c (or A d ) can be represented in the form of A + BF , where A is the open-loop state matrix, B is the control matrix, and F is the gain matrix. For simplicity, we use the term system to indicate dynamical system. Intuitively, a set S is called an invariant set for a system, if all the trajectories of the system, which are starting in S, remain in S. Numerous surveys on the theory and applications of invariant sets are published in the recent decades, see e.g., Blanchini [3]. Recently, several sufficient and necessary conditions, which are simply refereed to as invariance conditions, are derived to verify if a set is an invariant set for a continuous or discrete system. Various convex sets with different characteristics are considered as candidates for invariant sets. Invariance conditions for polyhedra are given in [5,6]. Ellipsoidal sets as invariant sets are analyzed in [4]. Cones as invariant sets are studied in [15,18,20]. A novel unified approach to derive invariance conditions for polyhedra, ellipsoids, and cones is presented in [13].
Although many mathematical techniques are developed to directly solve continuous systems, in practice, one usually solves a continuous system by applying certain discretization methods. Assume that a set is an invariant set for a continuous system, then it should be also an invariant set for the discrete system, which is obtained by the discretization method, i.e., discretization should preserve the invariance. However, this is not always true for every steplength used in the discretization method, thus it will be convenient if there exists a predictable threshold for valid invariance preserving steplength. The existence of such steplength thresholds of invariance preserving on various sets is thoroughly studied in [12]. In this paper, we consider three types of discretization methods on polyhedra and we aim to derive valid thresholds of the steplength in terms of explicit form or obtained by using efficiently computable algorithms. The popularity of polyhedra as invariant sets is due to the fact that the state and control variables are usually represented in terms of linear inequalities. For Taylor approximation type discretization methods, i.e., the coefficient matrix of the discrete system is derived from the Taylor expansion of e Ac∆t , we present an algorithm to derive a valid steplength threshold for invariance preserving. In particular, the algorithm aims to find the first positive zeros of some polynomial functions related to the system and the polyhedron. For general rational function type discretization methods, i.e., the coefficient matrix of the discrete system is a rational function with respect to A c and ∆t, we derive a valid steplength threshold for invariance preserving that can be computed by using analogous methods as for the case of Taylor approximation type methods. This steplength threshold is related to the steplength threshold for the forward Euler method and the radius of absolute monotonicity of the discretization method. We note that this result is similar to the one presented in [10,11], where Runge-Kutta methods are considered. Finally, we propose an optimization model to find the largest steplength threshold for the forward Euler method.
Notation: For the sake of simplicity, the following notational conventions are introduced. A nonnegative matrix, denoted by H ≥ 0, means that all entries of H are nonnegative. An off-diagonal nonnegative matrix, denoted by H ≥ o 0, means that all entries, except the diagonal entries, of H are nonnegative.
The paper is organized as follows. In Section 2, some fundamental concepts, theorems, and the key problems in this paper are introduced. In Section 3, we present our main results, i.e., deriving valid steplength thresholds for invariance preserving, for the three types of discretization methods. Finally, conclusions are provided in Section 4.

2.
Background. We now introduce the definitions of invariant sets for continuous and discrete systems.
According to the definitions of invariant sets, we have that an invariant set means that the continuous (or discrete) trajectory of the system remains in the same set. In fact, there is an alternative perspective, see e.g., [13]. In that interpretation S is an invariant set for (1) if and only if e Act S ⊆ S for any t ≥ 0, and S is an invariant set for (2) if and only if A d S ⊆ S.
In this paper, candidate invariant sets are restricted to convex polyhedron in R n . A polyhedron P in R n can be characterized as the intersection of a finite number of half spaces.
Two classical subsets of polyhedra are extensively studied in many applications. One is called polytope, which is a bounded polyhedron. The other one is called polyhedral cone, a polyhedron with b = 0 in (3), and the origin is its only vertex.
Given a system and a polyhedron, the invariance condition indicates sufficient and necessary condition such that the polyhedron is an invariant set for the system. There are many such equivalent invariance conditions, e.g., [2,5]. The most common ones are presented in Theorem 2.3. A novel and unified approach to derive these invariance conditions is proposed in [13]. The invariance conditions in Theorem 2.3 provide powerful and practical tools to verify whether a polyhedron is an invariant set for a given system.
• the discrete system (2) if and only if there exists anH ∈ R m×m , such that From the theoretical perspective, when a discretization method is applied to a continuous system, the invariant polyhedron for the continuous system should also be an invariant set for the discrete system. This means that conditions (4) and (5) are satisfied simultaneously, when the system, polyhedron, and discretization method are given. However, this is not always true. Intuitively, the smaller steplength used in the discretization method has larger possibility to yield that the polyhedron is also an invariant set for the discrete system. For the sake of self-contained presentation, the formal definitions of invariance preserving and steplength threshold are introduced as follows.
Definition 2.4. Assume a polyhedron P is an invariant set for the continuous system (1), and a discretization method is applied to the continuous system to yield a discrete system. If there exists a τ > 0, such that P is also an invariant set for the discrete system for any steplength ∆t ∈ [0, τ ], then the discretization method is invariance preserving for ∆t ∈ [0, τ ] on P, and τ is a steplength threshold for invariance preserving of this discretization method on P.
The steplength threshold in Definition 2.4 implies that any value smaller than this threshold is also a valid steplength threshold 1 . This is an important property. In certain cases, a discretization method may be invariance preserving on a set in the form of [0, Here we are only interested in finding τ 1 . We also note that the steplength threshold in Definition 2.4 is uniform 2 on P, i.e., τ needs to be a valid steplenth threshold for every initial point in P.
Since a continuous system is usually solved by using various discretization methods in practice, invariance preserving property of the chosen discretization method plays an important role. Further, a larger steplength threshold has many advantages in practice. For example, for larger steplength, the size of the discretized system is smaller, which yields that the computation is less expensive. Thus, we introduce the key problem in the paper: Find a valid (if possible the largest) steplength threshold τ > 0, such that a discretization method is invariance preserving for every ∆t ∈ [0, τ ] on P.
3. Main Results. In this section, we present the approaches for computing a valid (or largest) steplength threshold such that three classes of discretization methods are invariance preserving on a polyhedron. These three classes of discretization methods are considered in the following order: Taylor approximation type discretization methods, rational function type discretizatin methods, and the forward Euler method. The Taylor approximation type represents a family of explicit methods. The rational function type is an extended family of the Taylor approximation type, which also includes some implicit methods. The relationship between these discretization methods and the forward Euler method is also studied. Finally, for the forward Euler method, we derive the largest steplength threshold for invariance preserving.
3.1. Taylor Approximation Type Discretization Methods. We first consider the Taylor approximation type discretization methods. Note that the solution of the continuous system (1) is explicitly represented as x(t) = e Act x 0 , thus one can use the Taylor approximation to numerically solve the continuous system. The p-order Taylor approximation of e Ac∆t is given as follows: The discrete system obtained by applying the Taylor approximation type discretization methods is given as (6). In fact, the Taylor approximation type methods form a family of discretization methods. For example, p = 1 corresponds to the forward Euler method, p = 2 corresponds to the general Runge-Kutta 2nd order methods.
3.1.1. Existence of Steplength Threshold. Our approach to derive steplength threshold is based on the invariance conditions presented in Theorem 2.3. The basic ideas is that we build the relationship between these two invariance conditions of the continuous and discrete systems. In fact, conditions (4) and (5) are essentially linear feasibility problems [16]. The unknowns in the two invariance conditions are the matrix H andH given by (4) and (5), respectively. Thus, the key is to find relationship between those matrices.
We note that γ in Lemma 3.1 is not unique, e.g., any value greater than a valid γ is also valid. We will show more about the effect of γ to the steplength threshold in Section 3.2, and the way to derive a larger steplength threshold based on γ is also presented.
We note that p i=0 γ i f i (∆t) is equivalent to replacing I andĤ in (8) by 1 and γ, respectively. Then, according to (11), we have For implying that the right hand side of (12) equals to 1, thus (10) follows immediately. b). We note that for every i the first term of f i (∆t) given as in (9) is 1 i! ∆t i . Then we can write Thus, we have that there exists a τ i > 0, i.e., the first positive zero of thus we have f i (∆t) ≥ 0 for all ∆t ∈ [0, τ ] and i = 0, 1, ..., p. According to (8), and by noting thatĤ i ≥ 0 for any i = 1, 2, ..., p, we have thatH(∆t) ≥ 0 for all ∆t ∈ [0, τ ], where τ is defined by (14). Thus, we have proved that the first condition in (5) is satisfied. By recursively using HG = GA c , for any i, we have (15) Then, according to (15), and substituting (7) and (6), we havẽ Thus, we have proved that the second condition in (5) is satisfied.
Since H satisfies (4), we have Hb ≤ 0. Also, note that H =Ĥ − γI, thus we Then, according to (16) and (10), we havẽ Thus, we have proved that the third condition in (5) is satisfied. The proof is complete.
Lemma 3.2 presents an important relationship between the two matrices H and H corresponding to the continuous and discrete systems invariance conditions. This relationship is explicitly represented in (7), which is derived from the Taylor approximation (6). According to Lemma 3.2 and Theorem 2.3, we have the following theorem. Theorem 3.3. Assume a polyhedron P be given as in (3) is an invariant set for the continuous system (1), and a Taylor approximation type discretization method (6) is applied to the continuous system (1). Then, the steplength threshold τ > 0 as given in Lemma 3.2 is a valid steplength threshold for invariance preserving for the given Taylor approximation type discretization method (6) on P.
According to the proof of Lemma 3.2, we have that a valid τ requires f i (∆t) ≥ 0 for all ∆t ∈ [0, τ ] and all i = 0, 1, ..., p, where f i (∆t) given as (9). Since each f i (∆t) can be represented in the form of (13), the following corollary is immediate. Corollary 1. The value of τ given in Theorem 3.3 (or Lemma 3.2) is a valid steplength threshold for invariance preserving on P for the Taylor approximation type discretization methods (6). To compute τ , one needs to find the first positive zeros of finitely many polynomial functions in the form where α 1 , α 2 , ..., α q ∈ R and q ∈ N.
In fact, Lemma 3.2 can be extended to a more general case for polynomial approximation rather than Taylor type discretization methods.
The proof for Part b) is the same as the one presented for Part b) in Lemma 3.2, thus we are not presenting here.
3.1.2. Compute Steplength Threshold. We now consider the value of τ , i.e., the steplength threshold. In this section, we present an algorithm to numerically compute τ . In particular, this algorithm aims to find the first positive zero of a polynomial function in the form of (17).
The proof is complete.
In fact, the value t * given in Lemma 3.6 can be considered as one of the termination criteria of the algorithm to find the first positive zero of f (∆t), where f (∆t) is defined as (17).
The Sturm sequence {s i (t)} of f (t) and the Sturm Theorem presented in the following definition play a key role in our algorithm. The Sturm Theorem aims to give the number of real zeros of a univariate polynomial function in an interval by using the property of Sturm sequence on the end points of the interval.
where s i (t) is the negative of the remainder on division of s i−2 (t) by s i−1 (t).
For the sake of simplicity, we introduce the following definition and notation, which are used in the statement of the Sturm Theorem.
According to Lemma 3.6 and Theorem 3.9, we now propose our algorithm to numerically find the first positive zero of f (∆t) where f (∆t) is defined as (17). Let us denote #f [δ] the number of positive zeros of f (∆t) at interval [0, δ]. The value of #f [δ] can be computed by Sturm Theorem 3.9. The basic idea in our algorithm is by using the bisection method to shrink the interval, which contains the first positive zero of f (t), by 2 in each iteration. Our algorithm is presented as follows.
Step 1: [Initial Setting]: Set t l = t • , t r = t * , and ǫ be the precision. The correctness of the termination condition in Step 2 is ensured by Lemma 3.6. If neither of the termination conditions in Step 2 and 3 are satisfied, then it means that the first positive zero of f (t) exists and is located in the interval (t l , t r ). The second case in Step 4 means that the first positive zero of f (t) is located in the interval (t l , t m ). Analogously, the third case in Step 4 means that the first positive zero of f (t) is located in the interval (t m , t r ). In Step 5, we conclude that the first positive zero of f (t) is located in the interval (t l , t r ). Recall that we are interested to find a value τ , such that f (t) ≥ 0 for all [0, τ ], thus we return t l , i.e., the left end of the interval. Remark 1. If all coefficients σ i ≥ 0 for i = 1, 2, ..., p in (18), then the algorithm is also applicable to compute a valid steplength threshold for invariance preserving for the polynomial approximation (18).

Rational Function Type Discretization
Methods. The previous discussion is mainly about a steplength threshold for invariance preserving for a Taylor approximation type discretization methods as specified in (6). In this section, we consider more general discretization methods, which are refereed to as the rational function type discretization methods. To be specific, these discretization methods applying to the continuous system yield the discrete system where r(t) : R → R is a rational function defined as where λ 0 , λ 1 , ..., λ p ∈ R, µ 0 , µ 1 , ..., µ q ∈ R, and p, q ∈ N. It is clear that Taylor approximation type discretization methods belong to this type. Some implicit methods are also in this type, e.g., the backward Euler method, Lobatto methods [7], etc. The radius of absolute monotonicity of a function is extensively used in the analysis of positivity, monotonicity, and contractivity of discretization methods for ordinary differential equations, see e.g., [10,14,17].
Theorem 3.12. Assume r(t) is a rational function with r(0) = 1. Let ρ be the radius of absolute monotonicity of r(t). Assume a polyhedron P be given as in (3) is an invariant set for the continuous system (1), and the rational function type discretization method given as in (24) is applied to the continuous system (1). Then τ = ρ γ , where γ is given in Lemma 3.1, is a valid steplength threshold for invariance preserving of the rational function type discretization method given as in (24) on P.
Proof. The framework of this proof is similar to the one presented for Lemma 3.2. Since P is an invariant set for the continuous system, according to Theorem 2.3 and Lemma 3.1, there exists an H, and γ > 0, such that H + γI ≥ 0, HG = GA c , and Hb ≤ 0.
Then, according to Theorem 2.3, to ensure P is also an invariant set for the discrete system, we need to prove that there exists anH(∆t) ∈ R m×m , such that LetH(∆t) = r(H∆t). Now we prove thatH(∆t) satisfies (28).
For the first condition in (28), we use the Taylor expansion of r(t) at the value −ρ as By substituting t = H∆t into (29) we havẽ Since ρ is the radius of absolute monotonicity of r(t), we have r (i) (−ρ) i! ≥ 0 for all i. Also, according to (27), and ∆t ≤ ρ γ , i.e., ρ ∆t ≥ γ, so we have H + ρ ∆t I ≥ H +γI ≥ 0. Then we have (H + ρ ∆t I) i ≥ 0 for all i. According to (30), we haveH(∆t) ≥ 0 for ∆t ≤ ρ γ , thus the first condition in (28) is satisfied. For the second condition in (28), according to Definition 3.10, the second condition in (28) can be rewritten as (h(H∆t)) −1 g(H∆t)G = Gg(A c ∆t)(h(A c ∆t)) −1 , i.e., According to (25), we have By recursively using HG = GA c , for any i, j, we have According to (32) and (33), we have that (31) is true, i.e., the second condition (28) is satisfied. For the third condition in (28) we havẽ Thus, the third condition in (28) is also satisfied. The proof is complete.
The assumption r(0) = 1 in Theorem 3.12 is a fundamental condition for most discretization methods. This is since the steplength ∆t = 0, yielding that the coefficient matrix of the discrete system is the identity matrix.

Compute
Steplength Threshold. The steplength threshold given in Theorem 3.12 is related to ρ and γ. Recall that γ is given in Lemma 3.1, thus we only consider the computation of ρ.
Since r(t) is a rational function, all of its derivatives r (i) (t) have the same format, i.e., they are represented as quotients of two polynomial functions. Now recall that the radius of absolute monotonicity ρ is defined as r (i) (t) ≥ 0 for t ∈ [−ρ, 0]. This requires that the polynomial function in the numerator of r (i) (t) is nonnegative for t ∈ [−ρ, 0]. Thus, a valid ρ is the negative of the first negative real zero of this polynomial function. Then an algorithm similar to the one presented in Section 3.1.2 can be proposed to numerically compute ρ. We are not repressing the algorithm here due to the space consideration.

3.3.
Parameter of Steplength Threshold. According to Theorem 3.3 and Theorem 3.12, we have that the parameter γ plays an important role to derive a large valid steplength threshold. In this section, we consider the effect of γ to the steplength threshold.
3.3.1. Best Parameter. Let us first consider the case for Taylor approximation type discretization methods. By simple modification, we have that f i (∆t) defined in (9) can be written as which means that smaller γ will yield larger steplength threshold for Taylor type discretization method given as in (6). Similarly, according to Theorem 3.12, we also have that smaller γ will yield larger steplength threshold for the rational function type discretization methods (24). Thus we prefer the smallest possible γ, which in fact can be computed by solving the following optimization problem min{γ | H + γI ≥ 0, HG = GA c , and Hb ≤ 0}.
In optimization problem (35), the variables are H and γ, while G, A c and b are known, thus problem (35) is a linear optimization problem, which can be easily solved by existing optimization algorithms, e.g., simplex methods [1] or interior point methods [16]. In particular, if there exists an H ≥ 0 such that HG = GA c and Hb ≤ 0, then the optimal solution, denoted by γ * , of (35) is nonpositive. In this case, according to (34), we have f i (∆t) ≥ 0 for all ∆t ≥ 0. Then according to the proof of Lemma 3.2, we have that the steplength threshold for invariance preserving for Taylor approximation type discretization methods (6) on polyhedron P is infinity. Similarly, if γ * ≤ 0, according to Theorem 3.12, we have that the steplength threshold for invariance preserving for rational function type discretization methods (24) on polyhedron P is also infinity. Thus, we have the following theorem.
Theorem 3.13. If the optimal solution of (35) is nonpositive, then the steplength threshold for invariance preserving on the polyhedron P is infinity for Taylor approximation type discretization methods (6) and rational function type discretization methods (24).
One should note that the steplength thresholds given in Theorem 3.3 and Theorem 3.12 may not be the largest steplength thresholds. For example, for the Taylor approximation type discretization methods, we aim to find the first positive zeros of finitely many polynomial functions. In fact, the first positive zeros may not be the best in some cases. For example, if the function is given as f (∆t) = (∆t − 1) 2 (∆t − 2) 2 , then its first positive zero is 1. Then, by our methods, we have τ = 1. However, it is clear that f (∆t) ≥ 0 for any ∆t ≥ 0. Thus, in this case, we have τ = ∞.
If the first zero, ∆t * , of a function is a local minimum of this function, i.e., f ′ (∆t * ) = 0, then the first zero should not be used for computing the steplength threshold. This is since the function is tangent to the x axis at the first zero. To verify if a zero is a local minimum, one can check the first order and second order directives f ′ (∆t * ) and f ′′ (∆t * ). If f (∆t * ) = 0 and f ′ (∆t * ) < 0, then we can say that ∆t * is not a local minimum, and thus it is a valid positive zero. If f (∆t * ) = 0, f ′ (∆t * ) = 0, and f ′′ (∆t * ) > 0, we can say that ∆t * is a local minimum. Then we have to make ∆t to be larger, and use an algorithm similar to the one presented in Section 3.1.2 to find the next zero of f (∆t).

3.3.2.
Relation to the Forward Euler Method. The following lemma presents the relationship between γ that satisfies the constraints in (35) and the operator I + γ −1 A c on P. Recall that I + ∆tA c is the coefficient matrix of the discrete system by using the forward Euler method. Proof. " ⇒ " For x ∈ P, i.e., Gx ≤ b, we have Thus we have (I + γ −1 A c )x ∈ P, i.e., (I + γ −1 A c )P ⊆ P.
" ⇐ " We note that (I + γ −1 A c )P ⊆ P means that P is an invariant set for the following discrete system: Then according to Theorem 2.3, we have that there exists anH ∈ R m×m , such that H ≥ 0,HG = G(I + γ −1 A c ), andHb ≤ b. LetĤ = γH, and then we havê i.e., Thus replacingĤ − γI by H, the proof is complete.
3.4. Forward Euler Method. As illustration, we consider the simplest discretization method, the forward Euler method, in this section. For simplicity, a polytope, i.e., a bounded polyhedron, is chosen as the invariant set for the forward Euler method. A polytope can be defined in terms of convex combination of its vetices, i.e., where {x i } are the vertices of P. A sufficient and necessary condition under which a polytope is an invariant set for the continuous system is presented below.
Lemma 3.15. [13] The polytope P defined as in (36) is an invariant set for the continuous system (1) if and only if A c x i ∈ T P (x i ), for i = 1, 2, ..., ℓ, where T P (x i ) is the tangent cone 3 at x i , which can be given An alternative is presented by the following discussion. Equation (37) implies that Ax i is a feasible direction, i.e., x i + τ i A c x i ∈ P, for sufficiently small τ i > 0. Then we can formulate the following linear optimization problem: Optimization problems (41) and (42) are equivalent problems, i.e., we claim that τ i is equal to ǫ i . Observing that n j=1 β (i) j = 1 for the first constraint in (42), we have i.e., A c x i = j =i According to the argument for ǫ i above, we have the following theorem.
Theorem 3.16. Assume that the polytope P defined as in (36) is an invariant set for the continuous system (1), and the forward Euler method is applied to (1). Then, τ = min i=1,2,...,ℓ {ǫ i }, where ǫ i is defined as in Corollary 2, is the largest steplength threshold τ > 0 for invariance preserving of the forward Euler method on P.
Proof. For any x ∈ P, and ∆t ∈ [0, τ ], we have x + ∆tA c x = ℓ i=1 λ i (x i + ∆tA c x i ). According to Corollary 2 and 0 ≤ ∆t ≤ τ ≤ ǫ i , we have x i + ∆tA c x i ∈ P. Thus we have x + ∆tA c x ∈ P. The proof is complete.

4.
Conclusions. Many real world problems are studied by developing dynamical system models. In practice, continuous systems are usually solved by using discretization methods. In this paper, we consider invariance preserving steplength thresholds on polyhedron, when the discrete system is obtained by using special classes of discretization methods. We particularly study three classes of discretization methods, which are: Taylor approximation type, rational function type, and the forward Euler method.
For the first class of discretization methods, we show that a valid steplength threshold can be obtained by finding the first positive zeros of a finite number of polynomial functions. We also present a simple and efficient algorithm to numerically compute these positive zeros. For the second class of discretization methods, a valid steplength threshold for invariance preserving is presented. This steplength threshold depends on the radius of absolute monotonicity, and can be computed by analogous method as in the first case. For the forward Euler method we prove that the largest steplength threshold can be obtained by solving a finite number of linear optimization problems.