COMPUTATION OF LOCAL ISS LYAPUNOV FUNCTIONS WITH LOW GAINS VIA LINEAR PROGRAMMING

. In this paper, we present a numerical algorithm for computing ISS Lyapunov functions for continuous-time systems which are input-to-state stable (ISS) on compact subsets of the state space. The algorithm relies on a linear programming problem and computes a continuous piecewise aﬃne ISS Lyapunov function on a simplicial grid covering the given compact set excluding a small neighborhood of the origin. The objective of the linear programming problem is to minimize the gain. We show that for every ISS system with a locally Lipschitz right-hand side our algorithm is in principle able to deliver an ISS Lyapunov function. For C 2 right-hand sides a more eﬃcient algorithm is proposed.


1.
Introduction. The concept of input-to-state stability (ISS) was first introduced by Sontag [20] in the late 1980s and has soon turned out to be one of the most influential concepts for characterizing robust stability of nonlinear systems. Several results on ISS can be found in [20,21,22]. In [24], different equivalent formulations of ISS are given. In particular, it is shown that the ISS property is equivalent to the existence of an ISS Lyapunov function. The ISS notion is also useful in stability analysis of large scale systems. Stability of large interconnected networks of systems can be analyzed by means of ISS small gain theorems if all the subsystems are ISS [6,7,8,9,16]. Motivated by these results, in this paper we focus on the computation of ISS Lyapunov functions for small systems, as the knowledge of ISS Lyapunov functions immediately leads to the knowledge of ISS gains which may be used in a small gain based stability analysis.
Based on [23, Lemma 2.10-2.14], it was shown in [4] that ISS Lyapunov functions in implication form may be calculated using a Zubov approach. An alternative Zubov type approach was developed in [17]. In these two papers the problem of computing ISS Lyapunov functions was transformed into the problem of computing robust Lyapunov functions for suitably designed auxiliary systems such as the auxiliary system described in Remark 3.2. This robust Lyapunov function can be characterized by the Zubov equation, a Hamilton-Jacobi-Bellman type partial differential equation, which can be solved numerically [3]. However, this approach only yields a numerical approximation of a Lyapunov function but not a true Lyapunov function. This means that the computed function is close to a Lyapunov function in some appropriate norm, but due to numerical errors it may fail to satisfy some of the properties a Lyapunov function is supposed to have, for instance it may not be decreasing along solution trajectories everywhere in its domain. For discrete time systems, following the same auxiliary system approach, true Lyapunov functions can be computed by a set oriented approach [11]. This numerical approach, however, does not carry over to the continuous time setting. Moreover, the detour via the auxiliary system introduces conservatism, since the resulting Lyapunov function and ISS gains strongly depend on the way the auxiliary system is constructed.
We thus propose a linear programming based algorithm for computing true ISS Lyapunov functions without introducing auxiliary systems. The approach to use linear programming for the computation of continuous, piecewise affine Lyapunov functions was first presented in [18]. In [12], it was proved that for exponentially stable equilibria the approach proposed in [18] always yields a solution. This result was extended to asymptotically stable systems [13], to asymptotically stable, arbitrarily switched, non-autonomous systems [14], and to asymptotically stable differential inclusions [2]. The approaches proposed in these papers yield true Lyapunov functions on compact subsets of the state space except possibly on arbitrarily small neighborhood of the asymptotically stable equilibrium. Mainly inspired by [2], in this paper we will propose an analogous linear programming based algorithm for computing true ISS Lyapunov functions for locally ISS systems.
The paper is organized as follows. In the ensuing Section 2, we introduce the notation and preliminaries. In Section 3 we present our algorithm along with a couple of auxiliary results needed in order to formulate the constraints in the resulting linear program. The algorithm is formulated in two variants for Lipschitz continuous and C 2 right hand sides. Section 4 contains the main results of the paper: we prove that upon successful termination the algorithm yields an ISS Lyapunov function outside a small neighborhood of the equilibrium, and that successful termination is guaranteed if the system admits a C 2 ISS Lyapunov function and the simplicial grid is chosen appropriately. In Section 5, we illustrate our algorithm and results by two numerical examples.
2. Notations and preliminaries. Let R + := [0, +∞) and for a set C ⊂ R n denote by cl C and int C its closure and interior respectively. For a vector x ∈ R n we denote its transpose by x . The standard inner product of x, y ∈ R n is denoted by x, y . We use the standard norms x p := ( n i=1 |x i | p ) 1/p for p ≥ 1 and x ∞ := max i∈{1,2,...,n} |x i | and let B p (z, r) denote the open ball, closed ball of radius r around z in the norm · p , respectively. The induced matrix norm is defined by A p := max x p =1 Ax p . By u ∞,p = ess sup t≥0 u(t) p we denote the essential supremum norm of a measurable function u : The closed convex hull of vectors x 0 , x 1 , . . . , x m ∈ R n is given by This definition is independent of the numbering of the x i , that is, of the choice of the reference point x 0 .
A continuous function α : R + → R + is said to be positive definite if it satisfies α(0) = 0 and α(s) > 0 for all s > 0. A positive definite function is of class K if it is strictly increasing and of class K ∞ if it is of class K and unbounded. A continuous function γ : R + → R + is of class L if γ(r) is strictly decreasing to 0 as r → ∞ and we call a continuous function β : R + × R + → R + of class KL if it is of class K ∞ in the first argument and of class L in the second argument.
In this paper we consider nonlinear perturbed systems described by the ordinary differential equationẋ with vector field f : R n × R m → R n , f (0, 0) = 0, state x(t) ∈ R n , and perturbation input u(t) ∈ R m , t ≥ 0. The admissible input values are given by U R := cl B 1 (0, R) ⊂ R m for a constant R > 0 and the admissible input functions by u ∈ U R := {u | u : R + → U R measurable}. The solution corresponding to an initial condition x(0) = x 0 and an input u ∈ U R is denoted by x(·, x 0 , u). For our algorithmic construction of Lyapunov functions, we need certain regularity properties of f which also determine certain inequalities imposed in the algorithm. To this end, we require one of the following two hypotheses. (H1) The map f : R n × R m → R n is locally Lipschitz continuous. (H2) The map f is twice continuously differentiable. Given a compact set G ⊂ R n , as regards (H1) we fix the following notation: For each u ∈ U R , L x (u) is a Lipschitz constant of the map x → f (x, u), and for each x ∈ R m , L u (x) is a Lipschitz constant for the function u → f (x, u). Moreover, by (H1) L x (u), L u (x) may be chosen so that there exist constants L x and L u such that for all x ∈ G, u ∈ U R . The following definition specifies the stability property we are considering in this paper.
Observe that ISS implies that the origin is an equilibrium of (1) which is locally asymptotically stable for u ≡ 0. The function γ ∈ K ∞ is called an ISS gain.
It is known that the ISS property of (1) is equivalent to the existence of a smooth, i.e. C ∞ , ISS Lyapunov function for (1), see [23]. While this result guarantees the existence of smooth ISS Lyapunov functions our numerical techniques will not generate a smooth function. In the following we will numerically construct piecewise affine and thus nonsmooth Lyapunov functions. In order to define these functions, we need a generalized notion of the gradient and for our purpose Clarke's subdifferential turns out to be useful. Since we are exclusively dealing with Lipschitz functions, we can use the following definition, cf.
Definition 2.2. For a Lipschitz continuous function V : R n → R, Clarke's subdifferential is given by Now we can state the definition of a nonsmooth ISS Lyapunov function.
hold for all x ∈ int G, u ∈ U R and ξ ∈ ∂ Cl V (x). If G = R n and R = ∞ then V is called a global nonsmooth ISS Lyapunov function. The function β ∈ K ∞ is called Lyapunov ISS gain or briefly gain in what follows. The gain is called linear if β is linear.
Remark 2.4. The particular norms chosen in the formulation of the ISS property in (3) or of an ISS Lyapunov function in (6) do not play a role from the conceptual point of view: as all norms in R n are equivalent, different norms will only lead to different numerical values of the gains. The particular formulations we have chosen will turn out to be useful in deriving easy estimates, see the last paragraph of the proof of Theorem 4.1.
In order to simplify the algorithm to be proposed in this paper, we will restrict ourselves to Lyapunov functions which satisfy (6) with linear functions α(s) = s and β(s) = rs for some fixed r > 0. The following proposition shows that on compact subsets of the state space excluding a ball around the origin this can be done without loss of generality. Proposition 2.5. Let G ⊂ R n with G ⊂ cl int G and 0 ∈ int G be compact. If there exists a Lipschitz continuous ISS Lyapunov function W for system (1) on G, then for any > 0 and σ > 0 there exist positive constants C, r > 0 such that and with U R from Definition 2.3.
Proof. From the assumption, there exist α, β ∈ K ∞ such that W satisfies For the construction of C and r we now distinguish two cases.
Case 1. lim sup s→0 β(s)/s is bounded. In this case we define Then there exists a constant r > 0 satisfying Case 2. lim sup s→0 β(s)/s is unbounded. In this case we choose C as Then it is possible to find a constant r > 0 such that for all u ∈ U R satisfying Cβ( u 1 ) ≥ .
In both cases, a straightforward calculation shows that V (x) = CW (x) satisfies the desired inequalities.
Remark 2.6. It may not always be possible to choose = 0 in Proposition 2.5. However, in general the linear programming approach to computing Lyapunov functions only works outside a neighborhood of the origin, anyway, cf. Remark 3.3, such that the need to remove B 2 (0, ) does not introduce additional limitations into our approach.
3. The algorithm. In this section we introduce an algorithm to compute a local ISS Lyapunov function defined on a compact set G ⊂ R n with 0 ∈ int G and valid for perturbation inputs from U R ⊂ R m . The algorithms uses linear programming and the representation of the function on a simplicial grid in order to obtain a numerical representation as a continuous, piecewise affine function. By taking into account interpolation errors, the algorithm yields a true Lyapunov function, not only an approximative one.
3.1. Definitions. We recall the following basic definitions: A simplex in R n is a set of the form Γ = co{x 0 , x 1 , . . . , x j }, where x 0 , x 1 , . . . , x j are affinely independent. The faces of Γ are given by co{x i0 , x i1 , . . . , An n-simplex is generated by a set of n + 1 affine independent vertices. A collection S of simplices in R n is called a simplicial complex, if (i) for every Γ ∈ S, all faces of Γ are in S, (ii) for all Γ 1 , Γ 2 ∈ S the intersection Γ 1 ∩ Γ 2 is a face of both Γ 1 and Γ 2 (or empty). Some authors consider the empty simplex to be a face of Γ, so that the last statement in (ii) is superfluous, but this will have no relevance in the present paper. The diameter of a simplex Γ is defined as diam(Γ ) := max x,y∈Γ x − y 2 .
We now return to our problem. We assume that G ⊂ R n may be partitioned into finitely many n-simplices T = {Γ ν | ν = 1, 2, . . . , N }, so that T defines a simplicial complex. By assumption, we may also partition U R into m-simplices T u = {Γ u κ | κ = 1, 2, . . . , N u } defining a simplicial complex. We briefly write h x,ν = diam(Γ ν ), h u,κ = diam(Γ u κ ) and h x = max ν=1,2,...,N h x,ν , h u = max κ=1,2,...,Nu h u,κ . For each x ∈ G we define the active index set For the simplices T u , we assume additionally that for each simplex Γ u κ ∈ T u , the vertices of Γ u κ are in the same closed orthant. (12) Let P L(T ) denote the space of continuous functions V : G → R which are affine on each simplex, i.e., there are a ν ∈ R, w ν ∈ R n , ν = 1, 2, . . . , N , such that Let ∇V ν,k (k = 1, 2, . . . , n) denote the k-th component of the vector ∇V ν for every Γ ν ∈ T . Similarly, we define P L u (T u ). Observe that (12) implies that the map u → u 1 is contained in P L u (T u ).
Remark 3.1. The algorithm will construct an ISS Lyapunov function V ∈ P L(T ). In particular, this means that the inequality (8) has to be satisfied. To this end, observe that from the definition it follows that for any function V ∈ P L(T ) Clarke's subdifferential is given by Hence, for fixed x ∈ int G and u ∈ U R inequality (8) becomes Linearity of the scalar product in its first argument implies that this is equivalent to An inequality of this type will be used for ensuring (8) in the algorithm.
Then, using arguments similar to [17], it can be shown that a robust Lyapunov function for (17) is an ISS Lyapunov function for (1). However, it turns out that for computational purposes this detour via the auxiliary system is not efficient, as it leads to an algorithm in which two linear programs have to be solved and furthermore introduces conservatism into the estimates. Therefore we will not explicitly use the auxiliary system. Our approach is based on the algorithm from [2] for the computation of robust Lyapunov functions. The way we adapt it to the ISS case is, however, inspired by the structure of (17).
3.2. Interpolation errors. As in [14,2], the key idea for the numerical computation of a true Lyapunov function lies in incorporating estimates for the interpolation errors on T -and in this paper also on T u -into the constraints of a linear program. In this section we analyze the error terms we need for this purpose. Let The basic idea of the algorithm is to impose conditions on V ∈ PL(T ) in the vertices x i of the simplices Γ ν ∈ T which ensure that the function V satisfies the inequalities (5) and (8) for σ = 1 on the whole set G \B 2 (0, ). Note that V ∈ PL(T ) is completely determined by its values at the vertices of the simplices in T .
In order to ensure the properness condition (5), we impose the condition for every vertex x i ∈ Γ ν , V (0) = 0 and V ∈ PL(T ). According to (18), for x ∈ Γ ν we have Note that (19) and V (0) = 0 can only be true if the origin is a vertex of the grid. This we always ensure in our computations.
In order to make sure that V (x) satisfies (16) for all x ∈ Γ ν ⊂ G, u ∈ Γ u κ ⊂ U R via imposing inequalities in the node values V (x i ), we need to incorporate an estimate of the interpolation error into the inequalities. To this end, we demand that for all i = 0, 1, . . . , n, j = 0, 1, . . . , m. Here A ν,κ ≥ 0 is a bound for the interpolation error of f in the points ( Close to the origin the positive term ∇V ν 1 A ν,κ may become predominant on the left hand side of (20), thus rendering (20) infeasible. This is the reason for excluding a small ball B 2 (0, ) in the construction of V . Under certain conditions on f this problem can be circumvented by choosing suitably shaped simplices near the origin. In order to keep the presentation in this paper concise we do not go into details here and refer to [10], instead.
Since (20) will be incorporated as an inequality constraint in the linear optimization problem, we need to derive an estimate for A ν,κ before we can formulate the algorithm. For this purpose we introduce the following Proposition 3.4. Here, for a function g : R n × R m → R which is twice continuously differentiable with respect to its first argument, we denote the Hessian of g(x, u) with respect to x at z by For the first argument x ∈ Γ ν , let and let K x : U R → R + , K x , respectively, denote a bounded function and a positive constant satisfying max z∈Γν r,s=1,2,...,n In the next proposition which is proved in a similar way to [2, Proposition 4.1, Lemma 4.2 and Corollary 4.3], we state properties of scalar functions g : G ×U R → R or vector functions g : G ×U R → R p with respect to their first argument. Analogous properties hold with respect to the second argument.
(a) If g(x, u) is Lipschitz continuous in x with the bounds L x (u), L x from (2), then holds for all x ∈ G, u ∈ U R . (b) If g j (x, u) is twice continuously differentiable with respect to x with the bound H x (u) from (21) on its second derivative for some j = 1, 2, . . . , p, then Under the same differentiability assumption for all j = 1, 2, . . . , p, the estimate holds for all u ∈ U R by assuming the bounds from (22).
(c) The analogous estimates hold for interpolation in u when the Hessian and the second derivatives are computed w.r.t. u. The corresponding constants will be denoted by H u (x), L u (x), K u (x) and K u .
Proof. We only prove the estimate (25) which is an immediate consequence of (24) and the estimate The proof of (26) follows from the following observation. Let M ∈ R n×n , |M | the matrix obtained by taking the absolute value component-wise, r an upper bound for the absolute values of the entries in M and E the matrix with all entries equal to 1. Then we have M 2 ≤ |M | 2 ≤ r E 2 = nr.
3.3. The algorithm. Now we have collected all the preliminaries to formulate the linear programming algorithm for computing an ISS Lyapunov function V for (1). In this algorithm, we introduce the values V xi = V (x i ) as optimization variables. It is desirable to obtain an ISS Lyapunov function in which the influence of the perturbation represented by the value r in the gain β(s) = rs in (8) is as small as possible. Since in general there is a tradeoff between σ and r in (8), see also [15], here we fix σ = 1. The objective of the linear program will then be to minimize the number r in (8) for σ = 1. As explained in Remark 3.3, we only consider x satisfying x ∈ G \ B(0, ) for a small > 0. To this end we define the subsets Thus T is a triangulation of a compact set G away from the origin, with a distance of at least . We note that since 0 ∈ int G we may choose > 0 small enough so that the origin is separated from R n \ G, loosely speaking G \ G is inside of G. In the remainder of the paper we will assume that > 0 is chosen in this manner. In the following algorithm, we will only impose the conditions (18) and (20) in those nodes x i which belong to simplices Γ ∈ T . Moreover, we use the estimates of the interpolation errors A ν,κ obtained from Proposition 3.4. Here we distinguish between vector fields f satisfying (H1) and (H2), respectively. While the stronger assumption (H2) allows for improved estimates, the weaker assumption (H1) applies to a larger class of vector fields f .

Algorithm
We solve the following linear optimization problem.
Optimization problem: For all vertices x i of each simplex Γ ν ∈ T , all vertices u j of each simplex Γ u κ ∈ T u , one of the following conditions is required: where A ν,κ = L x h x,ν + L u h u,κ if f satisfies (H1) and A ν,κ = nK x h 2 x,ν + mK u h 2 u,κ if f satisfies (H2). Remark 3.5.
(i) By (19), the condition (A1) yields V (x) ≥ x 2 for x ∈ G ε . (ii) The condition (A2) ensures that V (x) does not explode and defines linear constraints on the optimization variables V xi , C ν,k . (iii) For (A3) note that ∂(G \ G ) is the inner boundary of G and ∂G is the outer boundary of G . Define M to be the maximum value of V (x) at the inner boundary. The linear constraint (A3) implies that the sublevel set {x ∈ G | V (x) ≤ M } ⊂ int G can be assumed to include the set B 2 (0, ) as we may assume V (x) < M in the interior of G \ G without violating any constraints. If system (1) is locally ISS, then the condition (A3) is not necessary.
Remark 3.6. If the above linear optimization problem has a feasible solution, then the values V xi = V (x i ) from this feasible solution at all vertices x i of all simplices Γ ν ∈ T and the condition V (x) ∈ P L(T ) uniquely define a continuous, piecewise affine function Remark 3.7. It follows from Proposition 3.4 that instead of the term nK x h 2 x,ν + mK u h 2 u,κ in (A4) when f fulfills (H2), one may use the sharper bound . The latter was used in our numerical experiments.
4. Main results. In this section we formulate and prove our two main results. We show that any feasible solution of our algorithm defines an ISS Lyapunov function on G and give conditions under which our algorithm will yield such a feasible solution.
We start with the former.  (5) and (6) for all x ∈ G and all u ∈ U R .
In order to prove inequality (8) for σ = 1 we compute According to Proposition 3.4, the constraints from (A4) ensure that V satisfies In the last step we have used the equality m i=0 µ j u j 1 = u 1 , which is true by assumption (12) and because we used the 1-norm for u. Indeed, this assumption ensures that the signs of the entries of u j coincide in each entry location. Thus we have shown (8) with σ = 1 for all x ∈ G and all u ∈ U R . Now we turn to the second objective of this section. We derive conditions under which the linear programming problem has a feasible solution. To this end, we need a certain regularity property of the simplices in our grids. In order to formalize these, we need the following notation.
For each Γ ν = co{x 0 , x 1 , . . . , x n } ∈ T , let y = x 0 , and define the n × n matrix X ν,y by writing the components of the vectors x 1 − y, x 2 − y, . . ., x n − y as row vectors consecutively, i.e., Let X * ν,y := X −1 ν,y 2 . In Part (ii) of the proof of [2,Theorem 4.6] it is proved that X * ν,y is independent of the order of x 1 , x 2 , . . . , x n . Moreover, X * ν,y = λ −1 min holds, where λ min is the smallest singular value of X ν,y . We define The regularity property now demands that we need to avoid grids with arbitrarily flat simplices. Formally, this means that there exists a positive constant R 1 > 0 such that all simplices Γ ν ∈ T in the considered grids satisfy the inequality for X * ν and X * from (33), cf. [2,Remark 4.7], and h x = max ν=1,2,...,N h x,ν from Section 3.1.
Theorem 4.2. Consider a system (1) which satisfies (H1) or (H2) and which is ISS. Let > 0 and R 1 > 0. Then, for all grids T and T u such that T satisfies (34) and for which h x and h u are sufficiently small, the linear programming problem from our algorithm has a feasible solution and delivers an ISS Lyapunov function V ∈ P L(T ) on G . Proof. Since system (1) is ISS, there exists a C 2 ISS Lyapunov function W : G → R, see [23]. Applying Proposition 2.5 we may without loss of generality assume that W satisfies (7) and (8) with σ = 2 and some r > 0. Now consider an arbitrary but fixed Γ ν = co{x 0 , x 1 , . . . , x n } ∈ T . Let y = x 0 , and define . . .
where H W (z) denotes the Hessian of W (x) at z and z i = y + ξ i (x i − y) for some ξ i ∈ [0, 1]. Thus, using Proposition 3.4 (by ignoring the dependence on u), we get where A := max z∈G i,j=1,2,...,n ∂ 2 W ∂xi∂xj (z) . Applying Proposition 3.4 again, we obtain After these preliminary considerations, we now assign values to the variables V xi and C ν,k of the linear programming problem from the algorithm and show that they fulfill the constraints.
For each vertex It thus remains to show (A4) for some r > 0.
Thus C ν,k ≥ |∇V ν,k | for each Γ ν ∈ T . Since ∇W (x) is bounded on G and (34) holds, there exists a positive constant C such that holds for all ν and k. From this analysis and the fact that W satisfies (8) with σ = 2, we obtain that where D := max x∈G,u∈U R f (x, u) 2 < ∞. Now let h = max{h x , h u }. If (H1) holds, i.e., in the Lipschitz case, the linear constraint from (A4) is fulfilled whenever h > 0 is so small that for all vertices x i of simplices in T ε nAh( In case (H2), i.e., f is C 2 , the linear constraint from (A4) is satisfied if h > 0 is so small that where K x , K u are the constants satisfying inequality (22) with respect to x, u respectively. Thus, the theorem is proved.

5.
Examples. In this section we illustrate the algorithm by two examples. In order to highlight the fact that our algorithm minimizes the gain r in (8) for σ = 1, in our first example we compare the result of our algorithm with two piecewise affine Lyapunov functions, for which a closed-form expression is derived following the construction of the proof of Theorem 4.2. Our second example shows the result of our algorithm for an example for which no closed-form ISS Lyapunov function is known.
Example 5.1. We consider the following system which is adapted from [19] x For this example, we obtain two ISS Lyapunov functions V 1 and V 2 on G based on two different functions W 1 (x), W 2 (x) following the construction in the proof of Theorem 4.2 and compare them with the numerical ISS Lyapunov function V delivered by the algorithm.
(1) For constructing a theoretical ISS Lyapunov function V 1 we start with the quadratic function candidate W 1 (x) = x 2 1 + x 2 2 . It is obvious that W 1 is twice differentiable. For this function we obtain for the dynamics from (42) is an increasing function whenever x 2 ∈ [0, √ 2 2 ] and can thus be extended as a K ∞ function for x 2 > √ 2 2 . Hence, for β(|u|) = 0.05|u| 2 ∈ K ∞ , W 1 is an ISS Lyapunov function for system (42) on G.
Now we follow the proof of Theorem 4.2 in order to construct a piecewise affine Lyapunov function V 1 satisfying the constraints in our algorithm. To this end, let = 0.048. Then the appropriate rescaling constant C in Proposition 2.5 is given by for x ∈ G , and |u| ≤ 4.41. For u satisfying |u| ≤ 4.41, we now need to find r 1 > 0 with Since in the algorithm the objective is to minimize r, we select the minimal r satisfying this inequality which is given by r 1 = 4.41 20 (1− 2 ) = 4.6044. Now, linear interpolation of this W 1 (x) on a sufficiently fine grid T yields the desired function V 1 (x) = C 1 W 1 (x), which is plotted in Figure 4.
(2) Since in the construction in the proof of Theorem 4.2 we rescale the function via Proposition 2.5 to satisfy W (x) ≥ x 2 , it appears reasonable to start with W 2 (x) = x 2 as a Lyapunov function candidate. Following the same steps as in (1), we can show that W 2 is also an ISS Lyapunov function. A rescaling V 2 (x) := C 2 W 2 (x) along Proposition 2.5 yields C 2 = 2 1− 2 and r 2 = 4.41 10(1− 2 ) = 0.4410. The resulting interpolated V 2 is shown in Figure 5.
(3) From the algorithm, we get the numerical ISS Lyapunov function V shown in Figure 6 with r = 0.420909. The simplicial complex is obtained in the same way as in [10,Section 2]: First, we pick points x = (x 1 , x 2 ) ∈ Ω ∩ Z 2 as vertices. The resulting triangulation of Ω is shown in Figure 1.
Second, let M 1 , M 2 , M 3 , M 4 be positive integers, typically much smaller than N 1 , N 2 , N 3 , N 4 , and ⊂ Ω. On the small neighborhood Ω 1 around the equilibrium, the previous triangulation is replaced by triangles with one vertex in the origin and two vertices on the boundary of Ω 1 . The resulting second triangulation of Ω for M 1 = . . . = M 4 = 2 is shown in Figure 2.
Third, by the map F : we transfer the vertices from the second triangulation into new vertices from which we construct the simplices used in computation, cf. V (x)}. The resulting set G is shown in Figure 3. In the map F , the parameter ρ > 0 controls the size of the resulting vertices. For our computations we used N i = 7, M i = 2, i = 1, 2, 3, 4, and ρ = 0.012.   The triangulation of U R is obtained by the same method in 1d with N 1 = N 2 = 21, M 1 = M 2 = 0 and using the map G : R → R with instead of F and γ = 0.01. As expected, the optimization based algorithmic approach yields the smallest possible gain parameter r. For finer grids, even smaller values of r can be obtained and it appears that r converges to a lower bound r > 0.4. However, we are not aware of an analytical method to confirm this numerically computed lower bound.
In Figures 8-9 we include a comparison of level curves of the calculated ISS Lyapunov function V with these of the theoretical obtained functions V 1 and V 2 . Note that the ISS Lyapunov function is not unique, but the calculated one is more similar to V 2 .   Example 5.2 (Synchronous generator with varying damping). We consider the following model adapted from [1] which is described bẏ The triangulation of U R is obtained with N 1 = N 2 = 5, M 1 = M 2 = 0, and the map from (47) with γ = 0.012. By the algorithm, we get the numerical ISS Lyapunov function V shown in Figure 10 for system (48). Some level curves of the ISS Lyapunov function V are shown in Figure 11. Note that for this example an analytical ISS Lyapunov function is not known and that our numerical analysis yields a numerical value for the (in our approach ) linear ISS gain.  Remark 5.3. In our algorithm, we construct grids on U R and G, respectively. If G is a two-dimensional set and the number of vertices for gridding U R increases by 1, then the number of constraints in the linear program increases at least by (48 + 6(N v +2)(N v −4)+12(N v −2)), where N v = min{N v,1 , N v,2 } and N v,i is the number of vertices intersecting the x i -axis. Similarly, the number of constraints increases in higher space dimensions. Thus, the gridding of U R renders the number of constraints much larger than the number of optimization variables. It is hence much faster to solve the corresponding dual optimization problem than to solve the primal problem. For numerical computations we used the GNU Linear Programming Kit (GLPK) 1 , Gurobi 2 and CPLEX 3 respectively. We experienced that Gurobi and CPLEX carry out a significantly better preprocessing of the constraints which eliminates much more redundant ones and thus both methods can solve the optimization problem much faster than GLPK.
Since f (x, u) and the interpolation errors on T and T u are incorporated in constraints (A4), the grids of x, u have an influence in computing the gain parameter. From numerical experience, the grid of x plays a more important role in obtaining a good estimate of r via the solution of the linear optimization problem (30). 6. Conclusions. In this paper, we proposed a new method of computing a continuous piecewise affine ISS Lyapunov function for a dynamic system with input perturbation (1). For suitable triangulations of the state space and the input perturbation space, the algorithm delivers a true ISS Lyapunov function with a gain function for system (1) on a compact subset of the state space minus a small neighborhood of the origin (Theorem 4.1). Such an ISS Lyapunov function satisfies a linear inequality. We think the new computational approach will be helpful in analyzing stability of interconnected systems which are locally ISS, which is one topic of our future research. If system (1) has a C 2 ISS Lyapunov function, then there exist triangulations such that the algorithm has a feasible solution (Theorem 4.2). It is known that if system (1) is ISS, there exists a smooth ISS Lyapunov function [23]. Therefore, our proposed algorithm always has a feasible solution.