A Lyapunov method for stability analysis of piecewise-affine systems over non-invariant domains

ABSTRACT This paper analyses stability of discrete-time piecewise-affine systems, defined on possibly non-invariant domains, taking into account the possible presence of multiple dynamics in each of the polytopic regions of the system. An algorithm based on linear programming is proposed, in order to prove exponential stability of the origin and to find a positively invariant estimate of its region of attraction. The results are based on the definition of a piecewise-affine Lyapunov function, which is in general discontinuous on the boundaries of the regions. The proposed method is proven to lead to feasible solutions in a broader range of cases as compared to a previously proposed approach. Two numerical examples are shown, among which a case where the proposed method is applied to a closed-loop system, to which model predictive control was applied without a-priori guarantee of stability.


Introduction
The research area of piecewise-affine (PWA) systems (Sontag, 1981), has received a lot of attention during the last 15 years. Different techniques have been recently studied for their stability analysis, mostly based on the computation, through semi-definite programming, of common quadratic or piecewise-quadratic Lyapunov functions (Eghbal, Pariz, & Karimpour, 2013;Feng, 2002;Ferrari Trecate, Cuzzola, Mignone, & Morari, 2002;Mason, Sigalotti, & Daafouz, 2012;Sun, 2010). Other methods are based on piecewise-polynomial Lyapunov functions (Prajna & Papachristodoulou, 2003), and on PWA Lyapunov functions (Grieder, Kvasnica, Baotić, & Morari, 2005). In spite of their simple formulation, PWA Lyapunov functions can in some cases find a solution where higher order piecewise-smooth functions fail to succeed. This is due to the non-conservative nature of the computational scheme employed for determining PWA Lyapunov functions: if such a function exists on a given partition, then the associated optimisation problem will be able to find it (this property is no longer valid for higher order functions).
In particular, PWA Lyapunov functions are obtained through linear programming (LP), imposing positivedefiniteness and decay conditions at the vertices of the polytopes that compose the (bounded) domain, therefore CONTACT Matteo Rubagotti mr@le.ac.uk; rubagotti@ieee.org † Present address: Université de Toulouse, LAAS, F- Toulouse, France.
enforcing the same properties for all the points of interest. Usually, the set (henceforth referred to as X ) where the PWA dynamics is defined is a known positively invariant set, because the notion of stability has no practical relevance if the state trajectory can exit the domain where the dynamics is defined (Biswas, Grieder, Löfberg, & Morari, 2005). However, there are many cases when the PWA system to be analysed is not defined in a positively invariant set. A typical example is when an explicit model predictive control (MPC) control law (Bemporad, Morari, Dua, & Pistikopoulos, 2002) is synthesised for a linear system without a-priori guarantees of stability and/or positive invariance of some suitable set for the closed-loop system. This can occur, for instance, when approximations of the optimal control law are introduced to obtain lowcomplexity solutions (see, e.g. Bemporad, Oliveri, Poggi, & Storace, 2011;Grieder et al., 2005;Jones & Morari, 2010).
In case of a non-invariant domain, a possible approach is to perform an extensive reachability analysis to find, through a recursive procedure, the maximum positively invariant (MPI) set included in X (see Blanchini & Miani, 2008, Ch. 4-5;Rakovic, Kerrigan, Mayne, & Lygeros, 2006 and references therein). Then, the Lyapunov stability analysis can be carried out in the MPI set. However, it often happens that this latter one is not a domain of attraction for the origin, because the domain of attraction is a proper subset of it. The domain of attraction (or a positively invariant subset of it) must then be determined in order to get a feasible solution applying one of the previously mentioned methods (Feng, 2002;Ferrari Trecate et al., 2002;Grieder et al., 2005;Prajna & Papachristodoulou, 2003). However, this procedure, when applied to PWA systems, can lead to computationally intractable solutions due to the exponential complexity of reachability analysis. Moreover, in many cases, searching for the MPI set is an undecidable problem.
An alternative solution is proposed in Rubagotti, Trimboli, Bernardini, and Bemporad (2011), where an invariant set is determined a posteriori by defining a fictitious dynamics that extends the actual dynamics of the system defined on X . In this way, a larger domain X e is considered, which is positively invariant for the extended system (i.e. the system including the regions where the fictitious dynamics is defined). If a PWA Lyapunov function can be determined for the extended system, a positively invariant (not necessarily maximal) set P e included in X is determined for the actual system. This solution would avoid performing an extensive reachability analysis. However, its main drawback arises from the arbitrariness in the definition of the fictitious dynamics. This can lead, for instance, to the artificial introduction of limit cycles on the extended dynamics, which would make it impossible to find a PWA Lyapunov function. This paper proposes a method based on PWA Lyapunov functions to assess exponential stability of the origin of a PWA system defined on a possibly non-invariant domain X , and determines an estimate P of the region of attraction contained in X . Even though the goal is the same as that of Rubagotti et al. (2011), no fictitious dynamics is required, and the PWA Lyapunov function (directly defined in P) and the set P itself are found simultaneously. Discontinuities at the boundaries of the polytopic sets are allowed for both the system dynamics and the PWA Lyapunov function. In particular, the system dynamics on the regions boundaries can assume any of the neighbouring values, and the presence of a finite number of affine dynamics is allowed within each polytopic region, aspects not considered in Rubagotti et al. (2011).
It is proven that the proposed method leads to a feasible solution whenever there exists a fictitious dynamics leading to a feasible solution by applying a generalisation of the method in Rubagotti et al. (2011). The focus of this paper is on (possibly multiple) nominal dynamics: the reader interested in the direct analysis of systems with parametric uncertainties or additive disturbances is referred to Trimboli, Rubagotti, and Bemporad (2011) and Rubagotti, Trimboli, and Bemporad (2013), respectively, and the references therein. A preliminary version of the theory here developed for the construction of the PWA Lyapunov function was described in Rubagotti, Zaccarian, and Bemporad (2012). As compared to Rubagotti et al. (2012), this paper analyses the presence of overlapping dynamics on the boundaries, uses a different construction of the LP which enlarges the domain of application of the proposed technique, and tests the new method on different examples. It will be shown that the complexity of the proposed LP (i.e. number of variables and constraints) increases only linearly with the order of the system. However, since typically the number of regions increases as well if more state variables come into play, the proposed method (like all the mentioned methods) is suitable to analyse small-sized systems (realistically up to 5-6 state variables) if a standard desktop computer is used. After introducing the main notation, definitions, and background results in Section 2, the stability problem is formulated in Section 3. The proposed LP-based PWA Lyapunov function construction is introduced in Section 4, and the related theoretical results are proven. Section 5 shows the broader applicability of the proposed method with respect to the method based on fictitious dynamics proposed in Rubagotti et al. (2011). Numerical results are discussed in Section 6, and conclusions are finally drawn in Section 7.

Notation, definitions, and background
Let R, R >0 , Z >0 , and Z ≥0 denote the sets of reals, strictly positive reals, strictly positive integers, and non-negative integers, respectively. Given a vector v ∈ R n , let |v| denote any vector norm. Given a discrete-time signal w : Z ≥0 → R p , the sequence of the values of w from the zero instant to the kth instant is denoted by w [k] . The norm of a sequence is defined as w [k] sup i∈{0,...,k} |w(i)|. Given a set D ⊆ R n , its interior is denoted by int(D), its closure byD, its boundary by ∂D, and its convex hull by conv(D). Given a finite number of sets A polyhedron is a set given by the intersection of a finite number of (closed or open) half-spaces. A polytope D is a bounded polyhedron, and the set of the vertices of its closureD is denoted by vert (D). With respect to a given vector norm, we define the norm ball of radius χ ∈ R >0 as B χ {a ∈ R n : |a| ≤ χ}. Given two sets D 1 , D 2 ∈ R n , their Minkowski sum is D 1 ⊕ D 2 {d 1 + d 2 : d 1 ∈ D 1 , d 2 ∈ D 2 }. A function γ : R ≥0 → R ≥0 of class K (γ ∈ K) is continuous, positive definite, and strictly increasing.
Consider a discrete-time nonlinear difference inclusion, where X ⊂ R n , and x ∈ R n is the state vector, X being a compact set that contains the origin in its interior. The notation , and will be always used in the following. Function ϕ( · ) is a locally bounded and outer semi-continuous set-valued mapping whose domain contains X . According to Goebel, Sanfelice, and Teel (2012, Ch. 6), these properties of ϕ ensure well posedness of (1) and suitable existence properties for solutions to (1).
The following ones are standard definitions for difference equations/inclusions. Due to the generality of inclusion (1), we characterise positive invariance as 'strong' , in the first definition below, to emphasise the fact that all solutions starting in D remain in D for all forward times.

Definition 2.1:
Remark 2.1: Note that set X is not assumed to be SPI for (1), so that some solutions can leave X and are therefore defined only on a bounded time domain. Definition 2.3: Consider dynamics (1) and an SPI set D ⊆ X with 0 ∈ int(D). The origin is exponentially stable in D (ES(D)) if there exist c, ρ ∈ R >0 , such that all solutions to (1) satisfy where η, α i ∈ R >0 , i = 1, 2, 3. Then, the origin is ES(D) for (1).

Proof:
The proof follows from Teel, Forni, and Zaccarian (2013, Theorem 1), by only concentrating on the jump dynamics (namely selecting the flow set to be empty). Note that Teel et al. (2013, Theorem 1) assumes continuous differentiability of U but this property is never used in the discrete-time part of the proof. Therefore, the same proof technique can be used to prove our Theorem 2.1.
Note that Theorem 2.1 allows U( · ) (called uniformly strict Lyapunov (USL) function, cf. Lazar, Heemels, & Teel, 2009) to be a discontinuous function. Continuity at the origin is implied by condition (2a), but the continuity on a neighbourhood of the origin is not required. In the remainder of the paper, discontinuous USL functions will be exploited to obtain theoretical results that avoid the conservativity possibly arising from enforcing continuity conditions.

Stability analysis problem and basic definitions
Let X ⊂ R n be a compact polytope that includes the origin in its interior. Consider a partition {X i } of X that consists of a finite number s of polytopes, where H i ∈ R n×q i , and h i ∈ R q i are constant vectors, and the inequality should be understood component-wise.
The number of vertices of X i is denoted by m i . The subset of indices I 0 is defined as I 0 {i ∈ I : 0 ∈ X i }, and it is assumed without loss of generality that, for all i ∈ I 0 , 0 ∈ vert(X i ).
The system dynamics can be defined by different affine functions g j (x) = A j x + a j , with A j ∈ R n×n , a j ∈ R n , j ∈ S {1, . . . , n g }. We denote S i ⊆ S as the subset of indices of functions g j associated to region X i : for example, if S 4 = {1, 5, 8}, it means that, if x ∈ X 4 , the state update can be defined according to g 1 , g 5 , or g 8 . The autonomous discrete-time PWA system is formally defined as Note that, being all X i defined as closed sets, (4) allows x to evolve according to any of the neighbouring dynamics when x is on a boundary shared by more regions. For instance, if x ∈ ∂X 1 and x ∈ ∂X 2 , then x + can be determined according to any g j with j ∈ S 1 ∪ S 2 . The idea of taking into account PWA systems with overlapping dynamics on the boundaries has already been considered by Hovd and Olaru (2013). Given the PWA system (4), for which X is not necessarily an SPI set, our goal is to prove the exponential stability of the origin and provide an SPI set P ⊆ X contained in its region of attraction. In order to introduce the Lyapunov stability analysis framework, suitable sets related to the possible transitions between regions are defined, as follows. We define the polytope, with R(X ) = i∈I R i , and where R i is the set of points reachable in one step from region X i , that can be computed (Blanchini & Miani, 2008) as The closure of the set X out X e \ X , namelyX out , is in general a non-connected set, that can always be expressed asX out = s i=s+1 X i , where X i are closed polytopes with (possibly) overlapping boundaries, but with disjoint interiors. An extended set of indices is also defined asĨ {1, . . . ,s}. In this way, new regions X i , i = s + 1, . . . ,s, are defined outside X , but no dynamics is associated with them. Notice that, in the particular case when X is an SPI set, then X e ≡ X , and X out = ∅. For each pair (i, k) ∈ I ×Ĩ, and each dynamics g j ∈ S i , we define the closed transition sets, of states that may be mapped into the polytope X k in one step from the polytope X i under dynamics g j . The number of vertices of each region X j ik is denoted by m j ik . The sets X j ik can be conveniently expressed as

PWA Lyapunov analysis
To the end of synthesising a USL function for system (4), where in (7a) F i ∈ R 1×n and f i ∈ R are free variables to be determined by the optimiser. Then, define V : X → R as Note that V i (x) and V j (x) may be different on points x ∈ X i ∩ X j . Therefore, for the values of x on the intersections, the required conditions on V i are imposed for all i ∈ N (x), although only the maximum value is taken in (7b), as V(x) must be single valued. The constraints are imposed for all the m i vertices v i,h ∈ vert(X i ), i ∈ I, h = 1, … , m i , where α 1 is a free parameter, while 0 < ϵ 1 is a fixed parameter. The use of ϵ arises from the fact that strict inequalities cannot be imposed in standard LPs. Condition (8a) will lead to V(x) ࣙ α 1 |x| in X , as will be formally shown in the remainder of the paper. The value of V( · ) at the vertices of each region X i is limited by imposing for all the m i vertices v i,h ∈ vert(X i ), i ∈ I, h = 1, … , m i , where M i for i ∈ I are free variables used by the optimiser. In order to obtain V(0) = 0, it is also required that Due to the boundedness of V(x), (8c) will make it possible to prove that there exists With constraints (8) in place, the vector of variables to be determined via a suitable optimisation is composed of α 1 , α 3 , and the terms M i , F i and f i , for all i ∈ I.
A procedure is now proposed so as to determine a choice for such variables by solving the following LP: subj. to (8).
Once (9) has been solved, the function V(x) is defined for all x ∈ X , based on (7). The following main result holds for the set: Theorem 4.1: Consider system (4), whose dynamics is defined on X , and assume that a bounded solution to problem (9) exists. Then, with reference to dynamics (4), P in (10) is an SPI set, and the origin is ES(P ).
Proof: Given x ∈ X i , for any i ∈ N (x) (see (7b)), define γ i, h ࣙ 0, such that m i h=1 γ i,h = 1, as a set of coefficients defining x as a convex combination of the vertices of X i . We obtain from (8a) that Considering that (8c) implies V i (0) = 0 for all i ∈ I 0 , and that V(x) is upper bounded by max i∈I M i < +∞, there exists a scalar α 2 ࣙ ϵ, such that V(x) ࣘ α 2 |x| for all x ∈ X . We conclude that Consider any x ∈ X I (i, j,k)∈I×S i ×I X j ik ⊆ X . Using again a convex combination by expressing Note that, after defining P as in (10), one has P ⊆ X I . The reason why is that, from (8e), one can easily obtain that F i x + g i ࣙ 1 for all x contained in the closure of X \ X I . Therefore, conditions (11) and (13) hold for all x ∈ P ⊆ X I ⊆ X . This fact leads to two conclusions: first, since P is a sublevel set of V(x), then from (13), P is an SPI set; second, from (11)-(13), applying Theorem 2.1, system (4) is ES(P ).

Remark 4.1:
Even if the formulation of the LP (9) will not lead to an explicit maximisation of the domain of attraction, the minimisation of the sum of the M i leads to values of the Lyapunov function as close as possible to zero, preventing them to be assigned large values. Intuitively, if more values of the PWA Lyapunov function on the vertices of the partition are below 1, then the set P is larger.

Remark 4.2:
The preliminary version of the algorithm here proposed (presented in Rubagotti et al., 2012), apart from considering a single affine dynamics for each region X i , would produce an unbounded LP solution in the special case in which X were already SPI. Therefore, in Rubagotti et al. (2012) it had to be assumed that X was not an SPI set. This assumption is not needed in this revised algorithm, since, in that special case, conditions (8e) do not appear in the LP, and the proof of Theorem 4.1 is still valid.

Remark 4.3:
It can be inferred from the description of the LP (9) that the total number of variables is n v = 2 + s(n + 2), while the total number of inequality constraints is

Remark 4.4:
In case (8) is infeasible, a possibility is to increase the number of regions X i , therefore providing more flexibility in synthesising the PWA Lyapunov function (7). A possible way is to consider the sets X j ik as the new sets X i and restart the one-step reachability analysis. As an alternative, one can employ a grid of the set X , which, redefining the partition together with the old sets X i , defines the new sets X i .

Remark 4.5:
In order to make the proposed approach more general, one might want to consider the case in which affine dynamics, additional with to those defined in the full-dimensional regions X i , are present defined only on given facets of dimension 1, . . . , n − 1 of one or more region X i (a thorough definition of the concept of facet can be found, for instance, in Spjøtvold, Kerrigan, Jones, Tøndel, & Johansen , 2006). In this way, specific (possibly multiple) dynamics can be defined on the lowerdimensional intersections between regions X i , including the case of dynamics defined only on one vertex. By considering each facet as a new and independent region, the results proved in Theorem 4.1 still hold.

Comparison with a previous approach
In Rubagotti et al. (2011), an analogous problem has been addressed by proposing a different solution. Since the approach of Rubagotti et al. (2011) considered the case of a single dynamics for each region X i , in order to make a comparison we will consider a particular sub-case of our approach in which S i = {i}, i.e. the i − th dynamics is associated with the i − th region only. Therefore, we can omit the index j from the subsequent formulation of this section. In Rubagotti et al. (2011), a contractive PWA fictitious dynamics (x) is defined a priori in the polytopic regions X i , i = s + 1, . . . ,s, as . . . ,s, (14) such that X e is an SPI set for the so-called extended system, defined as In Rubagotti et al. (2011), due to the difficulty of defining the function (x), it was suggested to define it as a contractive dynamics of type (x) = ρx, ρ ∈ R [0,1) , which was noticed to give good results in practice. In the following we generalise such an approach, by assuming that (x) can be any PWA dynamics. Then, a PWA Lyapunov function V e : X e → R is determined by LP, in order to satisfy Equation (2) with η = 1 for all x ∈ X e . More precisely, we have with F e i ∈ R 1×n and f e i ∈ R, and then we define V e as (16b) If a feasible realisation of function V e (x) can be determined, the set P e is defined as and system (4) is proven to be exponentially stable in P e . The fictitious dynamics provides an additional degree of freedom, but it is in general very hard if not impossible to know (except perhaps for very simple examples) what choice of the fictitious dynamics would lead to a larger set P e , or even if there exists a realisation of (x), such that a set P e = ∅ can be determined. Therefore, in some cases, a wrong choice of the fictitious dynamics can prevent the extended system from converging to the origin. The approach of this paper overcomes this problem, since no dynamics is defined out of the set X , and the Lyapunov function is defined only for x ∈ X . The next result on the comparison between the two methods is now introduced. Theorem 5.1: Given a PWA system in form (4), whose dynamics is defined on X , define the extended set X e as in (5), such that X e is an SPI set for the extended system (15). Assume that there exists a fictitious dynamics (x), x ∈X out , such that a PWA USL function V e (x) is determined for all x ∈ X e . Then, for system (4) defined on X , there exists a feasible solution of problem (9) and a scalar β ∈ R >0 , such that V(x) = βV e (x) for all x ∈ X .
Proof: According to the assumptions, for system (15), V e (x) is defined as a PWA function on X e . Moreover, there exist α e 1 , α e 2 , α e 3 ∈ R >0 , such that We will show that, starting from V e (x), we can always find another USL functionṼ e (x) defined for x ∈ X , which is a feasible solution of (9). First of all, define Then, defineṼ e (x) βV e (x), x ∈ X . It is now possible to state thatṼ e (x) = max i∈N (x)Ṽ e i (x), withṼ e i (x) F e i x +f e i , beingF e i βF e i andf e i β f e i . Note that V e (x) for x ∈ X satisfies, by construction, If any of the conditions (23) is satisfied for all x in a given compact set, it is automatically satisfied also on its vertices. Therefore, conditions (23) imply conditions (8), which means thatṼ e is a feasible solution of problem (9). Remark 5.1: Notice that the existence of a feasible solution of problem (9) does not necessarily imply that, once (x) is fixed, there exists a PWA USL function V e (x) in X e . For example, consider the following PWA sys- The fictitious dynamics in X 4 would be defined a priori as a contractive one. We assume that x + = 0.5x, leading to a positively invariant set X e = [−1, 4]. For any initial condition x(0) ࢠ (1, 4], the state evolution will exhibit a limit cycle, therefore preventing to prove the (in this simple case, apparent) fact that the origin is an exponentially stable equilibrium point with domain of attraction X 1 ∪ X 2 . This result is instead obtained with the approach presented in this paper: both X 1 and X 2 are SPI sets, while x + / ∈ X , ∀x ∈ X 3 . The transition sets would be X 1 11 = X 1 , X 2 22 = X 2 , and X 34 = X 3 . An optimal solution for (9) would be given by the PWA Lyapunov function, In conclusion, it is always possible to obtain P solving (9) when P e can be found with the method based on fictitious dynamics. On the other hand, there are systems for which this latter gives no solution (mainly because of the difficulty, in an n-dimensional case, of understanding what a reasonable choice for (x) can be), while it is possible to find a set P by using the method proposed in this paper.

Example 1
As a first example, a simple first-order piecewise-linear system is considered, which allows one to intuitively understand the possible advantages of the application of the proposed method. The system is defined as follows: x ∈ X 4 [5, 6].
One can immediately see that, because of the dynamics in X 4 , the set X = 4 i=1 X i is not SPI, and therefore classical methods for numerically determining a Lyapunov function cannot be directly applied. The method proposed in this paper, instead, would determine the SPI set P = [−2, 5) and the Lyapunov function shown in Figure 1. In this latter, notice that V(x) = 1 for x ∈ X 4 , while values smaller than 1 are achieved or the other 3 regions. The computation of the transition sets X  example, with j = k) and the other preliminary calculations took 0.36 s on a 2.4 GHz processor with 8Gb RAM. The LP for determining V(x) consists of 32 constraints and 14 variables, and is solved using CVX (Grant & Boyd, 2013) in 0.62 s on the previously-mentioned machine.
In this simple case, the set P coincides with the maximal positively invariant set, which can be determined intuitively, and for which an efficient computational algorithm is in general provided in the MPT toolbox (Herceg, Kvasnica, Jones, & Morari, 2013). After determining the MPI set, classical methods (based on PWA or piecewisequadratic Lyapunov functions implemented for instance in the MPT toolbox) can be successfully employed to prove the exponential stability of the origin. Notice that a common quadratic Lyapunov function does not exist in P for the considered example.

Example 2
In order to show a more complex example of the proposed approach, we consider the application of the PWA control law of Bemporad et al. (2002) to a discretised double integrator, where and where the goal is to optimally stabilise the origin satisfying the input constraints u(k) ∈ U for all k ∈ Z ≥0 , with U {u ∈ R : |u| ≤ 1}. In particular, following Bemporad et al. (2002), we seek for an optimal LQ stabiliser over a prediction horizon N = 5, with weight matrices on the state and input variables, respectively. No terminal weight or terminal constraints are defined, and no state constraints are imposed (as a consequence, no setinvariance or closed-loop stability properties are guaranteed a priori). The arising explicit MPC control law, determined by using the MPT toolbox, is defined as a PWA function for all x ∈ R 2 . If we consider the closedloop dynamics in X x ∈ R 2 : |x| ∞ ≤ 10 , we obtain a partitioning of X into 31 regions, to which 31 different affine dynamics are associated. In order to use a finer partition, we increase the number of regions X i by adding a rectangular partitioning, obtaining a total number of s = 148 regions. The closed-loop dynamics (4) is Since state constraints were not imposed a priori, X is not an SPI set. Therefore, we are interested in finding a set P ⊂ X which is SPI and contained in the region of attraction. As a first step, we find the set X e defined in (5), and the transition sets X j ik . Notice that, in this case, the closed-loop dynamics defined by the MPC controller is continuous, and j = k, since the dynamics is uniquely defined at each x ∈ X . We are then ready to formulate the LP (9), setting ϵ = 10 −5 . The resulting LP is infeasible. Therefore, we decide, as suggested in Remark 4.4, to grid the set X as it can be clearly seen in Figure 2 (a coarser grid would again lead to an infeasible LP). The LP is successfully solved, which proves that the origin of the closed-loop system is ES(P ) for the set P shown in Figure 2. Also, the PWA Lyapunov function V(x) is shown in Figure 3. The LP comprises 3108 constraints and 594 variables, and is solved using CVX in 1.15 s on a 2.4 GHz processor with 8GB RAM. The computation of the sets X j ik , together with all the other required operations apart from the LP, took about 75 s on the same computer.
The same stability analysis problem was considered also using the method proposed in Rubagotti et al. (2011), but no feasible solution was found, meaning that this example falls into the cases considered in Section 5. In this case, P does not necessarily coincide with the MPI set for the closed-loop system. The MPI set could not be numerically determined using the MPT Toolbox, but an overestimation of it can be determined using the same toolbox by computing the maximum control-invariant set for the open-loop system (24) with respect to the given state and input constraints. The result is shown in Figure 1, where the maximum control-invariant set is depicted in black.
This shows that the area of set P is close to that of the MPI set.

Conclusions
This paper proposes a convex optimisation-based method for the analysis of stability and invariance of PWA systems, which is based on the construction of a suitable PWA Lyapunov function defined through LP. Multiple affine dynamics are allowed in each region of the PWA system, and discontinuities are taken into account for both the system dynamics and the Lyapunov function at the regions boundaries. The method is of particular interest when the set of definition of the system dynamics is not positively invariant, but the case of positively invariant sets can be analysed with the same procedure, as a particular case. The method can find a feasible solution in a broader range of cases with respect to a previously-proposed approach, being independent from the definition of fictitious dynamics. Future work comprises extending our approach from mere analysis to controller synthesis. The main difficulty of this resides in preserving the convexity of the proposed optimisation.

Disclosure statement
No potential conflict of interest was reported by the authors.