Polyhedral Lyapunov Functions with Fixed Complexity

Polyhedral Lyapunov Functions with Fixed Complexity arXiv:2103.03613v1 [math.OC] 5 Mar 2021 Dimitris Kousoulidis and Fulvio Forni Abstract— Polyhedral Lyapunov functions can approximate any norm arbitrarily well. Because of this, they are used to study the stability of linear time varying and linear parameter varying systems without being conservative. However, the computational cost associated with using them grows unbounded as the size of their representation increases. Finding them is also a hard computational problem. Here we present an algorithm that attempts to find polyhedral functions while keeping the size of the representation fixed, to limit computational costs. We do this by measuring the gap from contraction for a given polyhedral set. The solution is then used to find perturbations on the polyhedral set that reduce the contraction gap. The process is repeated until a valid polyhedral Lyapunov function is obtained. The approach is rooted in linear programming. This leads to a flexible method capable of handling additional linear constraints and objectives, and enables the use of the algorithm for control synthesis. set. This is one of the features that distinguish our approach from others available in the literature [2], [3], [10], [14], [17]. We provide a comparison in Section VI. Another relevant feature is that we can easily extend the framework to include additional constraints and new variables, which is of practical significance for control synthesis. This allows us to develop most of our discussion for autonomous linear systems, simplifying the exposition. We then show how the algorithm is generalized to control design and nonlinear systems in Section V. We also show how our algorithm can be used to assess contraction/incremental stability of a system. The effectiveness of the algorithm is illustrated on a DC motor example, in Section VII. I. I NTRODUCTION We use C to denote a convex set and int C to denote its interior, which we define as the largest open subset of Rn in C. Inequalities on vectors and matrices are used in the element-wise sense. We represent polyhedra in two ways: V-representation given n × m matrix V , In this paper we provide a tool for the analysis and design of Linear Time Varying (LTV) and Linear Parameter Varying (LPV) systems. We first focus on verifying the stability of these systems. Then, we focus on the design of linear feedback for the stabilization problem. For the stability problem we adopt the framework of polyhedral Lyapunov functions [3]. Those are asymmetric norms (finitely valued, satisfy the triangle inequality, positively homogeneous, and positive definite) whose unit balls are polyhedra. We provide the necessary background on polyhedra and polyhedral Lyapunov functions in Section II, and show how to use them for verifying stability in Section III. The polyhedral Lyapunov framework is less conservative than standard quadratic Lyapunov functions. But this comes at a price: whereas finding a quadratic Lyapunov function to assess stability is a convex problem, efficiently solved using linear matrix inequalities [4], the same problem is hard to solve for polyhedral Lyapunov functions. In this paper we tackle the problem by developing an approach based on a two step iterated procedure, along the lines of [11]. Each step is a linear program (LP) and the iterations build a sequence of polyhedral sets, each reducing the gap to a contractive set, i.e. to a decaying Lyapunov function. The full algorithm is presented in Section IV. A relevant feature of our algorithm is that we fix the complexity of the polyhedral Lyapunov function, that is, the number of vertices or constraints that define the polyhedral D. Kousoulidis is supported by the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom. D. Kousoulidis and F. Forni are with the Department of Engineering, University of Cambridge, CB2 1PZ, UK [email protected] & [email protected] II. M ATHEMATICAL P RELIMINARIES A. Polyhedra CV (V ) = {x : x = V p, p ≥ 0, 1T p = 1} (1) H-representation given m × n matrix H, CH (H) = {x : Hx ≤ 1} (2) For this paper we exclusively use V-representation polyhedra. All results can be readily extended to H-representation polyhedra by duality, [3, Proposition 4.35]. Geometrically, the columns of the V correspond to the vertices of CV (V ), which is formed by taking their convex hull. We make the standing assumptions that our polyhedra • are Bounded (A1) • contain 0 in their interior (A2) Assumption (A1) is always true for V-representation polyhedra. Eq. (A2) requires that there exists a p+ > 0 such that 0 = V p+ , which we can test using linear programming (LP). B. Polyhedral Lyapunov Functions Any convex set C defines a non-negative function through its Minkowski functional. The Minkowski functional of a set C is given by ΨC (x) = inf{r ∈ R : r > 0 and x ∈ rC} (3) For V-representation polyhedra this means ΨV (x) = min{1T p : x = V p, p ≥ 0} (4) Using the strong duality of LP problems [5, Section 5.2.1], this is equivalent to ΨV (x) = max{hT x : hT V ≤ 1T } (5) Using LP, (4) and (5) can be evaluated numerically at any given point. For V-representation polyhedra satisfying (A2), ΨV is always a real valued function, strictly positive for all x 6= 0, and radially unbounded. As such they are valid Lyapunov function candidates. However, for us to fully establish the link with Lyapunov theory we also require the notion of a derivative of the candidate function. Unfortunately, the Minkowski functionals of polyhedra are not everywhere differentiable and we need to resort to the more general setting of subdifferentials: ∂Ψ(x) = {h : ∀y ∈ Rn , Ψ(y) − Ψ(x) ≥ hT (y − x)} (6) For differentiable functions the subdifferential at x corresponds to the gradient. For Minkowski functionals of V-representation polyhedra, we can obtain an explicit representations of ∂ΨV (x). Proposition 1 (V-representation Subdifferential): For any given matrix V , consider the Minkowski functional ΨV . Then, ∂ΨV (x) = {h : hT x = ΨV (x), hT V ≤ 1T } (7) Proof: (7) ⊆ ∂ΨV (x): We need to show that ΨV (y) ≥ hT y. From (5), this must be true. ∂ΨV (x) ⊆ (7): We prove this in two steps. We first show that any h ∈ ∂ΨV (x) must satisfy hT V ≤ 1T by contradiction. If h doesn’t satisfy hT V ≤ 1T , there must exist a column of V , [V ]i , such that hT [V ]i = c > 1 and ΨV ([V ]i ) = 1. If we set y = k[V ]i , for h ∈ ∂ΨV (x) we must then have k − ΨV (x) ≥ kc − hT x for all k > 0. However, if we set k > max(0, (hT x − ΨV (x))/(c − 1)) the inequality is invalidated, leading to a contradiction. As such, if h ∈ ∂ΨV (x) then hT V ≤ 1T . If we set y = 0 we get hT x ≥ ΨV (x). But from before hT V ≤ 1T , and so h must be included in the domain of (5) and therefore hT x = ΨV (x). Since h satisfies both hT V ≤ 1T and hT x = ΨV (x), h ∈ (7). We use subdifferentials to quantify how our candidate Lyapunov function ΨV changes along system trajectories. Proposition 2: For any trajectory x(·) of the autonomous system ẋ = Ax, Z t sup hT Ax(τ )dτ (8) ΨV (x(t)) − ΨV (x(0)) = 0 h∈∂ΨV (x(τ )) Proof: See [6] and [3, Sections 2.2]. The integrand is equivalent to max hT Ax subject to hT x = ΨV (x), hT V ≤ 1T , which is an LP that we can directly numerically compute. III. C ONTRACTION , L INEAR P ROGRAMMING , AND LYAPUNOV S TABILITY Convex sets can be used to show stability through the notion of contraction (also frequently referred to as strict positive invariance). For continuous time systems, this requires us to first define tangent cones. For a general convex set C, the interior of the tangent cone at a point x ∈ C is given by int KC (x) = {w : λ > 0, x + λw ∈ int C} (9) We then say that C is a contracting convex set under the dynamics ẋ = Ax if ∀x ∈ C, Ax ∈ int KC (x) (10) This can be interpreted in words as saying that the flow induced by the dynamics maps all points in C to points in int C for all t > 0. For polyhedra, we can express this geometric condition in a more convenient form. Proposition 3 (Contraction in V-representation): For a polyhedron CV (V ), (10) is equivalent to ∃M ∈ Rm×m and scalar η > 0 such that −η1T = 1T M, AV = V M, off-diag(M ) ≥ 0 [≡ M is Metzler] (11a) (11b) (11c) Where we use off-diag(M ) ≥ 0 to denote that all elements of M outside of the main diagonal are non-negative. We call η the rate of contraction. Proof: For CV (V ), we use (1) to express (10) as ∀(p ≥ 0, 1T p = 1) : ∃(λ > 0, p̂ ≥ 0, 1T p̂ < 1) : AV p = λ−1 V (p̂ − p) This is equivalent to ∃(Λ > 0, P̂ ≥ 0, 1T P̂ < 1T ) : AV = V Λ−1 (P̂ − I) where Λ is a diagonal matrix. Setting M̂ = Λ−1 (P̂ − I), we rewrite this as: ∃(η > 0, M̂ ) : 1T M̂ ≤ −η1T , AV = V M̂ , off-diag(M̂ ) ≥ 0 From (A2), ∃p+ > 0 such that 0 = V p+ . We can then add this p+ to each column of M̂ until it satisfies (11a) without affecting (11b) or (11c). Testing contraction is thus an LP optimization problem. max η M,η (12) s.t.: (11) The conditions of Proposition 3 are verified if the maximization in (12) is greater than 0. It is fairly straightforward to show that set contraction implies stability through the homogeneity of the linear dynamics. However, using (7), we can build a direct connection to Lyapunov theory. Proposition 4 (Lyapunov stability): CV (V ) is a contracting set for A with rate η if and only if ∀x ∈ Rn sup hT Ax ≤ −ηΨV (x) (13) h∈∂ΨV (x) Proof: (11) =⇒ (13): Starting with (7) and (4) we have that there exists some p such that hT x = 1T p = ΨV (x), x = V p, p ≥ 0 for the h ∈ ∂ΨV (x) that maximizes hT Ax. Also note that we can rewrite a matrix M satisfying (11c) as M = P + δI, where P is an positive matrix. Therefore hT Ax = hT AV p = hT V M p (11b) T = h V (P + δI)p = hT V P p + δhT V p = hT V P p + δ1T p [x = V p, 1T p = hT x] ≤ 1T P p + δ1T p = 1T M p [h ∈ ∂ΨV (x), P p ≥ 0] ≤ −η1T p (11a) ≤ −ηΨV (x) (13) =⇒ (11): We show this by contradiction. The optimization problem in (11) can be broken down to multiple multiple sub-problems, one for each column of V , [V ]i : ∃p ∈ Rm : η = −1T p, A[V ]i = V p, pj6=i ≥ 0 (14) where {j 6= i} is [1, m] \ {i}. If (11) is infeasible for a given η then (14) must be infeasible for at least one [V ]i . This is because otherwise a feasible solution for (11) could be generated by aggregating {p} into a matrix P . As such we begin by assuming that there exists some [V ]i for which (14) is infeasible. Then, by Farkas’ Lemma, there exists an h such that hT A[V ]i > −η, hT [V ]i = 1, hT V ≤ 1T (15) This means that h ∈ ∂ΨV ([V ]i ). But, since, ΨV ([V ]i ) = 1, hT A[V ]i > −ηΨV ([V ]i ), contradicting (13). From Lyapunov theory, (13) with η > 0 combined to (8) imply global exponential stability of the system with convergence rate η. The connection between contraction and Lyapunov stability is a well-known result in system theory [3, Thm 4.24, 4.33]. However, we include this proof to show how we can use (7) to directly proof results about polyhedral Lyapunov functions. We hope to use this approach to extend results that have no direct interprettation as set contraction, such as dissipativity, to polyhedral Lyapunov functions. This can be particularly useful for control design and for the analysis of open systems. IV. A LGORITHM FOR F INDING P OLYHEDRA A. Overview Testing if a set is contracting under the linear dynamics ẋ = Ax is a much easier problem than finding a contractive set. In fact, treating V as a variable makes (11b) nonconvex and nonlinear (since M and η are also variables). However, we can instead try to find such a V iteratively, by decomposing (12) into a pair of LP problems. One considers V fixed and estimates the contraction gap of V . The other, introduces variations on V to reduce the contraction gap. To initialize our algorithm we need to build an initial candidate polyhedron which satisfies our standing assumptions. A possible construction is suggested in Remark 1. The number of vertices m of our polyhedron will remain the same throughout the whole procedure. This means that we can control the complexity of our Lyapunov function and treat m as an input parameter. (i) Estimation of the contraction gap: for any given candidate polyhedron V , the contraction gap is estimated via the LP optimization (12) The following Proposition justifies the use of (12). Proposition 5: Given a CV (V ) that satisfies (A2), there exists some finite η (not necessarily > 0) that satisfies (11). Additionally, if CV (V ) satisfies (11) for some η ∗ , it will also satisfy (11) for any η < η ∗ . Proof: Note that (A2) implies that V is full row rank since otherwise the strict interior of CV (V ) would be empty. We rewrite p = pb +δi p+ , where pb is a solution to A[V ]i = V pb (implied to exist because V is full row rank) and p+ > 0 is a positive solution to 0 = V p+ (implied to exist because of (A2)). There then exists a δi∗ such that A[V ]i = V p and pi6=j ≥ 0 for all δi ≥ δi∗ . Since 1T p+ > 0, increasing δi decreases ηi = −1T pb − δi 1T p+ . As such, there exists {δi } such that for all i ηi∗ = − mini {1T pb + δi∗ p+ }. Aggregating {p} into a matrix P produces a solution to (11) with η = ηi∗ . By increasing {δi } further, (11) can be satisfied for any η < η∗ . If η > 0 we conclude that our polyhedron is contracting and we are done. (ii) Reduction of the contraction gap: if η ≤ 0, the following LP program searches for variations on V that are compatible with a reduction of the contraction gap. For a small ε > 0, solve: max δη (16a) −δη1 = 1T δM, (16b) δη,δM,δV T AδV = δV M + V δM, 0 ≤ off-diag(M + δM ) kδV k1 ≤ ε, kδM k1 ≤ ε (16c) (16d) (16e) The idea is to look for variations δV in a small neighborhood of V that would improve the gap η < η + δη. The search space is limited to the tangent space of the constraints manifold, at the point (η, V, M ). Therefore, for small enough variations, any solution to the LP above improves the contraction rate while preserving the integrity of the constraints. We thus build the new polyhedron as Vnew := V + δV . A small ε ensures that the Vnew remains in a small neighborhood of V to be compliant with the linear approximation above. It also guarantees that CV (Vnew ) satisfies our standing assumption (by continuity argument). The procedure iterates between these two steps until the estimation of the contraction gap returns a positive η. Remark 1: We compute an initial candidate polyhedron by randomly sampling m−1 points in Rn and rescaling them to have length 1. We then add an additional vertex equal to minus their sum and again normalised to have length 1. This way we obtain a polyhedron CV (V ) with m irreducible vertices that, as long as m ≥ n + 1, must also satisfy (A2). B. Fast reduction of the contraction gap We can speed up (16) by not having δη and δM as variables in the second optimization problem. We do this by considering directions in the tangent space that are compatible with the optimality conditions of (12) around a fixed previous solution. With these extra restrictions in place, we can obtain a directional derivative of η with respect to V. To estimate this derivative we need to rewrite (12) into standard LP form max ĉT x̂ (17a) x̂ s.t.: Âx̂ = b̂ x̂ ≥ 0 (17b) (17c) where  is a k × l matrix. A detailed transformation procedure is provided in Section IV-C. Using the standard form, we can then make use of the notion of Basic Feasible Solutions (BFS): from [15, Section 13.2], a BFS of (17) is a set B ⊂ {1, .., l} and a vector x̂ such that: • Âx̂ = b̂ and x̂ ≥ 0 • B contains k elements • i 6∈ B =⇒ x̂i = 0. • The k × k matrix B̂ defined as B̂ = [Âi ]i∈B is non-singular, where Âi is the ith column of Â. When solving (17) it is sufficient to only consider BFSs [15, Theorem 13.2]. We can obtain a BFS that is also an optimal solution of (17) via the Simplex algorithm [7]. Many commercial interior point solvers also offer the option of deriving a BFS as a cheap post-processing stage. With the optimal BFS we can now find the derivative of η with respect to changes in V . In general, we need to implicitly differentiate equation (17b) to obtain a derivative of x̂ (which represents our variables η and M in standard form) with respect to  and b̂ (which are functions of V ). However, an explicit solution can only be obtained if  is invertible [8]. Instead, we restrict ourselves to the BFS and consider only variables that aren’t constrained to be 0 by our choice of B. Indeed, given ŷ = {x̂i : i ∈ B}, we can always obtain the derivative of ŷ with respect to B̂ and b̂. Because B̂ is invertible (from the properties of the BFS), this derivative is always well defined. Using the transformation in Section IV-C, we can map the LP problem back to its original coordinates and build an explicit representation of DV η and DV M . We can then replace (16) with the following LP optimization problem: max (DV η)(δV ) subject to: kδV k1 ≤ ε δV and : (18) off-diag(M +(DV M )(δV )) ≥ 0 By constraining δV to be small enough we can also ensure that the new polyhedron CV (V +δV ) satisfies (A1) and (A2). This can also be enforced directly, as in [11], but keeping δV small seems to work sufficiently well in practice. We summarize the full algorithm below: Algorithm 1 Stability Verification Algorithm Data: The matrix A, the maximum number of iterations, and the number of vertices in the matrix m Result: V satisfying (11) if found, else F alse Procedure: V (0) = Initialize(m) k=0 while k ≤ max iter: η, M, B = Optimal BFS of (12) with V = V (k) if η > 0: return V (k) DV η, DV M = Derivatives of (12) using B and V = V (k) δV = (18) using M, DV η, and DV M V (k+1) = V (k) + δV k := k + 1 return F alse C. The Derivative Computation in Detail To go from (12) to (17) we use the vectorization operator vec(·), which converts a n × m matrix into a column vector with nm entries by stacking up its columns and the Kronecker product (⊗):     Im ⊗ V − vec(V ) 0nm 0nm vec(AV ) b̂ =  = 0m Im ⊗ 1Tm −1m 1m −1m     0mm vec(P )  0   −δ    ĉ =  x̂ =   1   η+  −1 η− (19) where η+ − η− = η and P + δI = M . We can also rewrite (17b) as f = Âx̂ − b̂ = B̂ ŷ − b̂ = Ĝ vec(V ) = 0 where Ĝ =   M T ⊗ In − Im ⊗ A 0m×(nm) (20) (21) From (20), we use the implicit function theorem to compute the Jacobian of ŷ with respect to vec(V ) ∇vec(V ) ŷ = −(∇ŷ f )−1 ∇vec(V ) f = −B̂ −1 Ĝ (22) such that (DV ŷ)(δV ) = (∇vec(V ) ŷ) vec(δV ). For our derivatives of interest (DV η)(δV ) = (DV η+ )(δV ) − (DV η− )(δV ) and, ∀i 6= j, (DV Mij )(δV ) = (DV Pij )(δV ) Note that all derivatives on the right hand side are of elements of x̂. To compute each of them we first obtain the corresponding index i in x̂ from (19). Then, if i 6∈ B, the derivative is zero. Else, we obtain the corresponding index j in ŷ so that x̂i = ŷj . To compute the derivative, we then multiply the j th row of (22) with vec(δV ). Our derivative calculations can be interpreted as only considering the effect of changes in V to our variables under a fixed basis. As long as the solution with the current basis is feasible for V + δV , this will still be a BFS (although it might not be the optimal one). V. G ENERALIZATIONS A. Feedback Design Algorithm 1 can be readily extended for output-feedback design for systems of the form ẋ = Ax + Bu y = Cx. This includes state feedback as a special case when C = I. We replace (12) with max η M,K,η (23) (A + BKC)V = V M, off-diag(M ) ≥ 0 [≡ M is Metzler] The two step procedure outlined in Section IV reduces the nonlinear (for non-fixed V ) program (12) into an iteration involving linear programs. A similar reduction can be applied to (23). Just like (12), (23) is a linear program for fixed V , which means that the estimation of the contraction gap is a straightforward operation. Likewise, for (23), the reduction of contraction gap step is performed on a larger manifold, which now also contains K. However, the search procedure remains essentially the same. This means that the approach in (IV-B) can also be adapted to control synthesis (it just requires different transformations to get the constraints into standard form). We can also consider scenarios beyond simple stabilization. New formulations may be considered as long as they are encoded into linear constraints. Natural directions to explore are testing/enforcing input/output gains, the design of dynamic output feedback, and the design of nonlinear state-feedback controllers taking advantage of the Lyapunov framework. This will be the goal of future research. B. From Linear Stability to LDIs The discussion in this paper has been developed for simple linear systems but, since (12) and (23) are convex, the same approach can be used for the analysis LTV and LPV systems, and even for nonlinear systems. The idea is that systems of the form ẋ = A(w(t))x + B(w(t))u (24) satisfy, for all t ≥ 0, A(w(t)) ∈ ch{Ai }, B(w(t)) ∈ ch{Bi } where ch{Ai } = {A : A = X i pi Ai , pi ≥ 0, −η1T = 1T Mi and Ai V = V Mi . Indeed, the extension of the constraint set does not change the structure of the solution. Combined with (23), this allows us to tackle the design of a stabilizing linear state-feedback law u = Kx. A similar approach can be explored to verify contraction [12], [18], [16], [9] of nonlinear systems of the form ẋ = f (x). The problem reduces to finding a contracting polyhedral set V for the linearization ˙ = ∂f (x)δx δx ∀x . The solution is obtained as above, by relaxing the problem to a suitable convex hull ∂f (x) ∈ ch{Ai } for all x. VI. C OMPARISON WITH R ELATED A PPROACHES −η1T = 1T M, s.t.: and likewise for ch{Bi }. Stability is thus verified as in (12) where constraints (11a) and (11b) are replaced by several constraints of the form X i pi = 1}; In systems theory polyhedral Lyapunov functions are usually constructed using algorithms that progressively add new vertices or constraints until a contractive set is reached. Other algorithms exploit optimization based approaches to directly construct contractive sets. [3], [10] and [14] provide examples of the iterative approach. The new vertices are generated by iterating the update equation of discrete-time linear systems. The main advantage of these approaches is that in many cases the sequence of polyhedra produced can be guaranteed to converge to a polyhedral Lyapunov function, if one exists. The main disadvantage is that not all linear constraints can be incorporated easily. Additionally, it is difficult to bound the complexity of the resulting polyhedron. Examples of optimization-based algorithms can be found in [2] and [17]. [17] considers positive rescalings of the vertices of the initial candidate polyhedron. This allows for fast computation but may lead to high complexity polyhedra even when a simple contractive polyhedron exists. The algorithm in [2] makes clever use of both representations of a polyhedron. However, this is also a limiting factor because of the possible difference in the number of vertices and constraints of the two representations. In higher dimensions, a few vertices/constraints in one representation may lead to a very large number of constraints/vertices in the other representation [13], which hinders the use of the algorithm in high dimensional settings. The approach proposed in this paper is also optimizationbased. It takes into account explicitly the problem of complexity of the polyhedron by taking the number of vertices of the polyhedron as a parameter. This comes at the price of weaker convergence guarantees (if the predefined number of vertices is too small the problem is unfeasible). Our algorithm generalizes the approach of [17] by allowing vertices to move freely, and, in constrast to [2], only uses a single representation. Our approach also tackles the problem linear control synthesis which is not considered by any of the approaches above. Using our algorithm we find that we simultaneously stabilise these under a common Lyapunov function using output feedback law   u = −19.7 −18.7 y (30) 0.2 i 0.1 0.0 The identified polyhedron is also shown in Fig. 1 and has only 9 vertices. −0.1 −0.2 −2.0−1.5−1.0−0.5 0.0 0.5 1.0 1.5 2.0 VIII. C ONCLUSIONS AND F UTURE W ORK ω Fig. 1: Left: Invariant polyhedron for (25) and (26), for γs = 10 (η ≥ 0.07). Right: Invariant polyhedron for (27) and (29), with γp = 4 and u as in (30) (η ≥ 0.004, 0.03 for nominal) R EFERENCES VII. E XAMPLE : DC M OTOR We test our optimization-based approach on the DC motor speed and position examples from [1]. For the DC motor speed example, we carry out a stability analysis under various possible parameter variations. We then solve the feedback design problem for the DC motor position example using both static and dynamic state feedback. We use the parameters given in the website for the speed example as our nominal values throughout (inertia J0 = .01, viscous friction constant b0 = .1, emf constant K0 = .01, resistance R0 = 1, and inductance L0 = .5). The DC motor speed example is governed by the following dynamics:      d ω −b/J K/J ω = (25) −K/L −R/L i dt i In this case, we consider 2 orders of magnitude of parameter perturbations of the form J ∈ [J0 /γs , γs J0 ], b ∈ [b0 /γs , γs b0 ], K ∈ [K0 /γs , γs K0 ] (26) for γs = 10. Using Algorithm 1, we are able to find a polyhedron satisfying (13) and (26) (with η = .07). It is shown in Fig. 1 and has only 6 vertices. For comparison, we are only able to find a common quadratic function using LMIs up to around γs = 8.7. We next consider the DC motor position example. The dynamics are now:        θ 0 1 0 θ 0 d    ω = 0 −b/J K/J  ω  +  0  u (27) dt i 0 −K/L −R/L i 1/L And we assume that we can only measure θ and i:     θ 1 0 0   ω y= 0 0 1 i (28) We also have parameter perturbations of the form J ∈ [J0 /γp , γp J0 ], K ∈ [K0 /γp , γp K0 ], with γp = 4. We presented an algorithm for finding polyhedral Lyapunov functions of fixed complexity for feedback design. We believe that this can be used to tackle interesting research directions. We plan to extend the theory to open systems and interconnections, and apply it to the nonlinear contraction setting. (29) [1] University of Michigan control Tutorials for MATLAB and Simulink. https://ctms.engin.umich.edu/CTMS/index. php?aux=Home. [2] R. Ambrosino, M. Ariola, and F. Amato. A Convex Condition for Robust Stability Analysis via Polyhedral Lyapunov Functions. 50(1):490–506. [3] F. Blanchini and S. Miani. Set-Theoretic Methods in Control. Systems & Control: Foundations & Applications. Birkhäuser, 2 edition. [4] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix Inequalities in System and Control Theory. Studies in Applied and Numerical Mathematics. Society for Industrial and Applied Mathematics. [5] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press. [6] F. Clarke. Lyapunov Functions and Feedback in Nonlinear Control. In M. S. Queiroz, M. Malisoff, and P. Wolenski, editors, Optimal Control, Stabilization and Nonsmooth Analysis, volume 301 of Lecture Notes in Control and Information Sciences, pages 267–282. Springer Berlin Heidelberg. [7] G. B. Dantzig. Origins of the simplex method. In S. G. Nash, editor, A History of Scientific Computing, pages 141–151. ACM. [8] A. L. Dontchev and R. T. Rockafellar. Implicit Functions and Solution Mappings: A View from Variational Analysis. Springer Monographs in Mathematics. Springer New York. [9] F. Forni and R. Sepulchre. A Differential Lyapunov Framework for Contraction Analysis. 59(3):614–628. [10] N. Guglielmi, L. Laglia, and V. Protasov. Polytope Lyapunov Functions for Stable and for Stabilizable LSS. 17(2):567–623. [11] D. Kousoulidis and F. Forni. An Optimization Approach to Verifying and Synthesizing K-cooperative Systems. In 21st IFAC World Congress (IFAC-V 2020). [12] W. Lohmiller and J. Slotine. On Contraction Analysis for Non-linear Systems. 34(6):683–696. [13] P. McMullen. The maximum numbers of faces of a convex polytope. 17(2):179–184. [14] S. Miani and C. Savorgnan. MAXIS-G: A software package for computing polyhedral invariant sets for constrained LPV systems. In Proceedings of the 44th IEEE Conference on Decision and Control, pages 7609–7614. [15] J. Nocedal and S. J. Wright. Numerical Optimization. Springer Series in Operations Research. Springer, 2nd ed edition. [16] A. Pavlov, N. van de Wouw, and H. Nijmeijer. Uniform Output Regulation of Nonlinear Systems: A Convergent Dynamics Approach. Systems & Control: Foundations & Applications. Birkhäuser, 2005. [17] A. Polański. On absolute stability analysis by polyhedral Lyapunov functions. 36(4):573–578. [18] G. Russo, M. Di Bernardo, and E. Sontag. Global entrainment of transcriptional systems to periodic inputs. PLoS Computational Biology, 6(4):e1000739, 04 2010.