Learning Lyapunov Functions for Hybrid Systems
Shaoru Chen, Mahyar Fazlyab, Manfred Morari, George J. Pappas, Victor M. Preciado
LLearning Lyapunov Functions for Hybrid Systems
Shaoru Chen, Mahyar Fazlyab, Manfred Morari, George J. Pappas, Victor M. Preciado ∗ Abstract
We propose a sampling-based approach to learn Lyapunov functions for a class of discrete-time autonomous hybrid systems that admit a mixed-integer representation. Such systemsinclude autonomous piecewise affine systems, closed-loop dynamics of linear systems with modelpredictive controllers, piecewise affine/linear complementarity/mixed-logical dynamical systemin feedback with a ReLU neural network controller, etc. The proposed method comprises analternation between a learner and a verifier to find a valid Lyapunov function inside a convexset of Lyapunov function candidates. In each iteration, the learner uses a collection of statesamples to select a Lyapunov function candidate through a convex program in the parameterspace. The verifier then solves a mixed-integer quadratic program in the state space to eithervalidate the proposed Lyapunov function candidate or reject it with a counterexample, i.e., astate where the Lyapunov condition fails. This counterexample is then added to the sample setof the learner to refine the set of Lyapunov function candidates. By designing the learner andthe verifier according to the analytic center cutting-plane method from convex optimization, weshow that when the set of Lyapunov functions is full-dimensional in the parameter space, ourmethod finds a Lyapunov function in a finite number of steps. We demonstrate our stabilityanalysis method on closed-loop MPC dynamical systems and a ReLU neural network controlledPWA system.
Hybrid systems have become widespread within the systems and control community in the lastdecades thanks to their flexibility in modeling the interaction of continuous and discrete dynamicalsystems that frequently arise in cyber-physical systems (CPS) [1, 2]. As today’s cyber-physicalsystems are getting more complex, developing new methods for design, analysis, and control ofhybrid systems is increasingly important.Analysis and control design for general hybrid systems is challenging and therefore, variousmethods have been proposed to tackle special classes of hybrid systems such as linear complemen-tarity (LC) systems [3, 4], mixed logical dynamical (MLD) systems [1], and piecewise affine (PWA)systems [5]. While these classes are mathematically equivalent [6], their representation could havea high impact on their numerical tractability. When it comes to stability analysis, various methodshave been proposed for PWA systems [7, 8, 9] while methods that directly deal with MLD and LCsystems are relatively scarce. Nevertheless, stability analysis tools for PWA systems are applicableto MLD and LC systems as they can be transformed into PWA systems [6].Transforming various types of hybrid systems into a PWA representation for stability analysismay not always be efficient. For example, for a PWA system in feedback with a ReLU neural ∗ Shaoru Chen, Manfred Morari, George J. Pappas, and Victor M. Preciado are with the Department of Electricaland Systems Engineering, University of Pennsylvania. Email: {srchen, morari, pappasg, preciado}@seas.upenn.edu.Mahyar Fazlyab is with the Mathematical Institute for Data Science, Johns Hopkins University. Email: [email protected] a r X i v : . [ m a t h . O C ] D ec etwork controller, although the closed-loop dynamics is PWA, identifying the PWA representationmay be tedious and the stability analysis task may become very challenging since the number ofpartitions generated by the ReLU network can be very large [10].In this paper, we propose a learning-based approach to stability analysis of hybrid systems thatadmit a mixed-integer formulation. These systems include PWA, MLD, LC systems and ReLUnetworks. Our method comprises a learner and a verifier, which iteratively search for a Lyapunovfunction from a target class F of Lyapunov functions (e.g., quadratic or piecewise quadratic). Ineach iteration, the learner uses a set of samples of the hybrid system to localize F by a convexset ˜ F ⊇ F and then solves a semidefinite program (SDP) to select a Lyapunov function candidatefrom ˜ F . The verifier then solves a mixed-integer program in the state space to either validate theLyapunov function candidate or reject it with a counterexample, i.e., a state where the Lyapunovcondition fails. This counterexample is then added to the sample set of the learner to refine theset ˜ F of Lyapunov function candidates. By designing the alternation between the learner and theverifier according to the analytic center cutting-plane method (ACCPM), we show that when theset F of Lyapunov functions is full-dimensional and contains a norm ball with radius (cid:15) > O ( n /(cid:15) ) steps, where n is the ambient dimension. LMI-based stability analysis of PWA systems : Among various stability analysis methods for PWAdynamical systems [8, 9], linear matrix inequality (LMI)-based approaches are relatively prominent.These methods construct an SDP whose solution gives a valid Lyapunov function. For continuous-time PWA systems, LMI-based approaches to synthesize piecewise affine [11], piecewise quadratic(PWQ) [12] and piecewise polynomial Lyapunov functions [13] have been proposed. The adaptationof these Lyapunov function synthesis methods to handle discrete-time PWA systems is summarizedin [7]. For discrete-time PWA systems, a common feature of the LMI-based methods is computingthe transition map between all pairs of modes. This step may become time-consuming when thenumber of modes is large.
Sampling-based Synthesis Methods : The iterative approach of alternating between a learning moduleand a verification module to synthesize a certificate for control systems is known as the Counter-Example Guided Inductive Synthesis (CEGIS) framework proposed by [14, 15] in the verificationcommunity. The application of CEGIS to Lyapunov function synthesis for continuous-time non-linear autonomous systems can be found in [16, 17, 18] using Satisfiability Modulo Theory (SMT)solvers for verification. In general, the termination of the iterative procedures in these works isnot guaranteed. Notably, Ravanbakhsh et al. [19] apply the CEGIS framework to synthesize con-trol Lyapunov functions for nonlinear continuous-time systems and provide finite-step terminationguarantees for the iterative algorithm through careful design of the learner which essentially im-plements the maximum volume ellipsoid cutting-plane method [20]. Our work differs from [19] inthe algorithm design and the application on hybrid systems. Other than the CEGIS framework,a learning-based approach to synthesize control barrier functions for hybrid systems is proposedin [21].
We denote the set of real numbers by R , the set of integers by Z , the n -dimensional real vectorspace by R n , and the set of n × m -dimensional real matrices by R n × m . The standard inner productbetween two matrices A, B ∈ R n × m is given by h A, B i = tr( A > B ) and the Frobenius norm of2 matrix A ∈ R n × m is given by k A k F = (tr( A > A )) / . Denote S n the set of n × n -dimensionalsymmetric matrices, and S n + ( S n ++ ) the set of n × n -dimensional positive semidefinite (definite)matrices. Given a set S ⊆ R n x + n y , Proj x ( S ) = { x ∈ R n x |∃ y ∈ R n y s.t. ( x, y ) ∈ S} denotes theorthogonal projection of S onto the subspace R n x . We denote int( S ) the set of all interior pointsin S . Consider a discrete-time autonomous hybrid system x + = f ( x ) (1)where x ∈ R n x is the state and f : R n x R n x is a continuous function. Without loss of generality,assume system (1) has an equilibrium at the origin, i.e., 0 = f (0). Let R ⊂ R n x be the domain ofthe system and R is a compact set which contains the origin in its interior. Denote x k the state ofsystem (1) at time k and x the initial state. The nonlinear dynamics x + = f ( x ) with domain R can be equivalently described through its graph defined asgr( f ) = { ( x, y ) ∈ R n x × R n x | x ∈ R , y = f ( x ) } . (2)In this paper, we study the Lyapunov stability of the origin of (1) for a class of hybrid systemsthat admit a mixed-integer formulation. Definition 1 (Mixed-integer formulation of a set [22]) . For a set
Q ⊂ R n z , consider the set L Q ⊆ R n z × R n λ × Z n µ in a lifted space given by L Q = { ( z ∈ R n z , λ ∈ R n λ , µ ∈ Z n µ ) | ‘ ( z, λ, µ ) ≤ v } , (3) with a function ‘ : R n z × R n λ × Z n µ → R n ‘ and a vector v ∈ R n ‘ . The set L Q is a mixed-integerformulation of Q if Proj z ( L Q ) = Q . If the function ‘ is linear, we call the related formulationmixed-integer linear (MIL). In the next subsections, we show how to find mixed-integer formulations for PWA, MLD, LCsystems and ReLU networks.
Consider a discrete-time piecewise affine system with control inputs x + = ψ i ( x, u ) = A i x + B i u + c i , ∀ ( x, u ) ∈ R i , (4)where R i = { ( x, u ) ∈ R n x × R n u | F i x + G i u ≤ h i } for i ∈ I = { , , · · · , N mode } are polyhedralpartitions of the state-input space R = S i ∈I R i . We assume that the partitions R i are bounded forall i ∈ I . To make the PWA system (4) well-posed, we assume that int( R i ) ∩ int( R j ) = ∅ , ∀ i = j and f i ( x, u ) = f j ( x, u ) , ∀ ( x, u ) ∈ R i ∩ R j if the intersection is not an empty set.We denote the PWA dynamics (4) collectively as x + = ψ ( x, u ). The graph of ψ ( x, u ) is givenby gr( ψ ) = S i ∈I gr( ψ i ), where each graph gr( ψ i ) is defined asgr( ψ i ) = { ( x, u, x + ) | Q i [ x > u > x > + ] > ≤ q i } , (5)3ith Q i = A i B i − I − A i − B i IF i G i , q i = − c i c i h i . (6)In this paper, we apply a disjunctive programming-based formulation [22] which states that x + = ψ ( x, u ) is equivalent to the following set of constraints [22] F i x i + G i u i ≤ µ i h i , µ i ∈ { , } , ∀ i ∈ I , X i ∈I µ i , x = X i ∈I x i , u = X i ∈I u i x + = X i ∈I ( A i x i + B i u i + µ i c i ) . (7)In the disjunctive programming formulation (7), we have N mode binary variables µ i and 2 N mode auxiliary continuous variables { x i , u i } . The binary variable µ i can be interpreted as the indicatorof the mode where the state-input pair ( x, u ) lives in. Note that µ i = 1 imposes µ j = 0 , ∀ j = i .By the boundedness of the partitions R j , we have F j x j + G j u j ≤ µ j h j = 0 ⇒ x j = 0 , u j = 0,and correspondingly x = P k ∈I x k = x i , u = P k ∈I u k = u i , x + = P k ∈I ( A k x k + B k u k + µ k c k ) = A i x i + B i u i + c i . We can also obtain a mixed-integer formulation of the PWA dynamics (7) throughthe big-M method [23].When the PWA system (4) is interconnected with a controller which also has a mixed-integerformulation, we can describe the closed-loop dynamics through mixed-integer constraints. Whenan autonomous PWA system is considered, we obtain its mixed-integer formulation by removingthe control input related variables in (7). Consider a discrete-time linear complementarity system [6, 3] x + = Ax + B u + B w (8a) v = E x + E u + E w + g (8b)0 ≤ v ⊥ w ≥ v, w ∈ R s and ⊥ denotes that v > w = 0. The complementarity constraint (8c) can beequivalently formulated as a set of mixed-integer linear constraints through the big-M method [23]as v i w i = 0 ⇔ ≤ v i ≤ µ i M ,i , ≤ w i ≤ (1 − µ i ) M ,i , (9)where µ i ∈ { , } and the subscript i denotes the i -th entry of the variable v and w . The binaryvariable µ i either forces v i to be zero ( µ i = 0) or forces w i to be zero ( µ i = 1). In general, selectingthe big-M values M ,i and M ,i is no simple task [24]. In practice, when the big-M values are hardto obtain, we can alternatively use the special-ordered set constraint in Gurobi [25] which forcesconstraint (8c) through a branching rule instead of specifying the big-M’s explicitly.One interesting application of the formulations in (8) and (9) is in the description of the closed-loop MPC systems. We refer the readers to [26] for the details of this formulation.4 .3 Mixed-logical dynamical systems The mixed-logical dynamical system [1] can be written as x + = Ax + B u + B δ + B z (10) E x + E u + E δ + E z ≤ g (11)where x = [ x > r x > b ] > with x r ∈ R n r and x b ∈ { , } n b ( u has a similar structure), z ∈ R n z and δ ∈ { , } n δ are auxiliary variables. The MLD system is explicitly constructed through mixed-integer linear constraints. Consider a PWA system in feedback with an L -layer ReLU neural network controller u = π ( x ),where π : R n x → R n u is given by z = xz ‘ +1 = max( W ‘ z ‘ + b ‘ , , ‘ = 0 , · · · , L − π ( x ) = W L z L + b L . (12)Here z = x ∈ R n ( n = n x ) is the input to the neural network, z ‘ +1 ∈ R n ‘ +1 is the vectorrepresenting the output of the ( ‘ + 1)-th hidden layer with n ‘ +1 neurons, π ( x ) ∈ R n L +1 ( n L +1 = n u )is the output of the neural network, and W ‘ ∈ R n ‘ +1 × n ‘ , b ‘ ∈ R n ‘ +1 are the weight matrix and thebias vector of the ( ‘ + 1)-th hidden layer, respectively.Consider a scalar ReLU function y = max(0 , x ) where x ≤ x ≤ ¯ x . Then it can be shown thatthe ReLU function admits the following mixed-integer representation [27], y = max(0 , x ) , x ≤ x ≤ ¯ x ⇐⇒ y ≥ , y ≥ x, y ≤ x − x (1 − t ) , y ≤ ¯ xt, t ∈ { , } , (13)where the binary variable t ∈ { , } is an indicator of the activation function being active ( y = x )or inactive ( y = 0). For a ReLU network described by the equations in (12), let m ‘ and ¯ m ‘ be the element-wise lower and upper bounds on the input to the ( ‘ + 1)-th activation layer, i.e., m ‘ ≤ W ‘ z ‘ + b ‘ ≤ ¯ m ‘ . Then the neural network equations are equivalent to a set of mixed-integerconstraints: z ‘ +1 = max( W ‘ z ‘ + b ‘ , ⇐⇒ z ‘ +1 ≥ W ‘ z ‘ + b ‘ z ‘ +1 ≤ W ‘ z ‘ + b ‘ − diag( m ‘ )( − t ‘ ) z ‘ +1 ≥ z ‘ +1 ≤ diag( ¯ m ‘ ) t ‘ , (14)where t ‘ ∈ { , } n ‘ +1 is a vector of binary variables for the ( ‘ + 1)-th activation layer. We notethat the element-wise pre-activation bounds { m ‘ , ¯ m ‘ } can be found by, for example, interval boundpropagation or linear programming assuming known bounds on the input of the neural network[28, 29, 30]. Combined with the MIL formulation of PWA systems shown in Section 2.1, theclosed-loop dynamics of a ReLU neural network controlled PWA system admits a mixed-integerformulation. 5 Stability Analysis via Lyapunov functions
The convergence behavior of system (1) around its equilibrium points can be studied by Lyapunovstability.
Definition 2. (Lyapunov stability [31, Chapter 13]) The equilibrium point x = 0 of the autonomoussystem (1) is• Lyapunov stable if for each (cid:15) >
0, there exists δ = δ ( (cid:15) ) such that if k x k < δ , then k x k k < (cid:15), ∀ k ≥ asymptotically stable if it is Lyapunov stable and there exists δ > k x k < δ ,then lim k →∞ k x k k = 0.Since the hybrid dynamics (1) is nonlinear, Lyapunov stability is often a local property and itis of interest to estimate its region of attraction defined as Definition 3 (Region of attraction) . The region of attraction O of the nonlinear system (1) isthe set of states from which the trajectory of system (1) converges to the origin, i.e., O = { x ∈R| lim t →∞ x t = 0 } . In this work, we do not assume the domain R to be positive invariant. Instead, we introducethe region of interest (ROI) X ⊆ R which is a polytopic set given by X := { x ∈ R n x | F X x ≤ h X } , (15)to guide our search for an inner approximation ˜ O of the ROA O . We can certify the asymptoticstability of system (1) and find an ˜ O by constructing Lyapunov functions defined in the followingtheorem: Theorem 1. [31, Chapter 13] Consider the discrete-time nonlinear hybrid system (1). If there isa continuous function V ( x ) : R 7→ R with domain R such that V (0) = 0 and V ( x ) > , ∀ x ∈ X \ { } (16a) V ( f ( x )) − V ( x ) ≤ , ∀ x ∈ X , (16b)where the set X is the region of interest (ROI) satisfying X ⊆ R and 0 ∈ int( X ), then the origin isLyapunov stable. If, in addition, V ( f ( x )) − V ( x ) < , ∀ x ∈ X \ { } , (17)then the origin is asymptotically stable.We call any V ( · ) satisfying (16a) a Lyapunov function candidate . If additionally, V ( · ) satisfiesthe condition (16b) or (17), then V ( · ) is called a valid Lyapunov function candidate, or simply,a Lyapunov function. Since asymptotic stability is our primary focus in this paper, Lyapunovfunctions refer to any V ( · ) satisfying constraints (16a) and (17) unless specified otherwise. Once aLyapunov function V ( x ) is obtained, an inner estimate of the ROA is given by ˜ O = { x | V ( x ) ≤ τ } ,where τ = inf x ∈R\X V ( x ). In other words, ˜ O ⊂ X is the largest sublevel set of V ( x ) that iscontained in X . Remark 1.
The ROI X can be set according to prior knowledge or from simulation of the nonlinearhybrid dynamics (1) . See Section 5 for examples of choosing the ROI. .1 Lyapunov function parameterization Searching for a Lyapunov function in the function space is intractable since the problem is infinite-dimensional in this space. Instead, we reduce our search space to the class of Lyapunov functionsdefined by V ( k ) ( x ; P ) = z ( k ) > P z ( k ) , (18)where P ∈ S ( k +1) n x ++ , and z ( k ) is given by z ( k ) = [ x > f ( x ) (1) , > f (2) , > ( x ) · · · f ( k ) , > ( x )] > (19)Here we use the notation f (0) ( x ) = x, f ( k +1) = f ( f ( k ) ( x )) , ∀ k ≥
0. We call V ( k ) ( x ; P ) the Lya-punov function candidate of order k , which is a quadratic function in composition with the systemdynamics f ( x ) evolved for k steps. Indeed, we can increase the complexity of the function classmonotonically by increasing the order k .Since P is positive definite, searching for a valid Lyapunov function of the form (18) reduces tofind a matrix P that satisfies the Lyapunov difference condition (17). Explicitly, we can characterizethe space of matrices P that admit a valid Lyapunov function candidate as F = { P ∈ S ( k +1) n x | αI (cid:22) P (cid:22) βI, ∆ V ( k ) ( x, P ) < , ∀ x ∈ X \ { }} . (20)where 0 ≤ α < β , and ∆ V ( k ) ( x, P ) is the Lyapunov difference given by∆ V ( k ) ( x, P ) := V ( k ) ( f ( x ); P ) − V ( k ) ( x ; P ) . (21)The constraint αI (cid:22) P guarantees the condition (16a), while the constraint P (cid:22) βI is imposed tomake F bounded without loss of generality since we can always scale P while satisfying (17).We call F the target set. It follows that F is convex since it is defined by semidefinite as wellas linear constraints on P . As αI (cid:22) P (cid:22) βI imposes an additional constraint on the conditionnumber of P , in practice we choose β/α large or simply set α = 0 . Then finding a Lyapunovfunction in a given function class V ( k ) ( x ; P ) is stated as the following problem: Problem 1.
For each parameterized function class V ( k ) ( x ; P ) and the target set F defined in (20) ,find a feasible point in F or certify that F is empty. Although the target set F is convex, the fact that it is characterized by infinitely many linearconstraints (i.e., the constraints ∆ V ( k ) ( x, P ) < , ∀ x ∈ X \ { } ) poses computational challenges tosolving Problem 1. In this work, we propose a learning-based approach to address this challenge byiteratively drawing state samples to refine our over-approximation of the target set F . By designingthe learning strategy based on ACCPM from convex optimization, we show that when the targetset F is full-dimensional in the parameter space of P , our method is guaranteed to find a feasiblepoint in F in a finite number of steps. Remark 2.
The parameterization of V ( k ) ( x ; P ) is inspired by the finite-step Lyapunov func-tion [32, 33] and the non-monotonic Lyapunov functions [34] which also construct Lyapunov func-tion candidates using the system states several steps ahead. V ( k ) ( x ; P ) allows us to parameterizecomplex function classes with a relatively small number of parameters. For example, when f ( x ) is a PWA function with N mode (possibly large) partitions, V (1) ( x ; P ) is a PWQ function with thesame partitions in the state space. As will be shown next, the proposed method always finds a feasible solution in the interior of F . Therefore,choosing α = 0 does not affect the positive definiteness of the solution. Learning Lyapunov functions from counterexamples
We recall from the previous section that finding a feasible point of the convex set F is computation-ally intractable since the condition ∆ V ( k ) ( x ; P ) < x ∈ X \ { } . To overcome thisintractability, we adopt a learning-based approach, in which we first select a set of finite samples S = { x , x , · · · , x N } ⊂ X and then enforce the linear constraint ∆ V ( k ) ( x ; P ) < x ∈ S . This results in an over-approximation of F given by˜ F = { P ∈ S ( k +1) n x | αI (cid:22) P (cid:22) βI, ∆ V ( k ) ( x, P ) ≤ , ∀ x ∈ S} . (22)We call ˜ F the localization set. Finding a feasible point in ˜ F now becomes a tractable convexfeasibility problem. However, there is no guarantee that a P ∈ ˜ F would correspond to a valid Lya-punov function, even if the number of samples approaches infinity. As one of our contributions, wepropose an efficient learning strategy based on the analytic center cutting-plane method (ACCPM)to iteratively grow the sample set and refine the set ˜ F until we find a feasible point in F . In thenext subsections, we describe the proposed approach and provide finite-step termination guaranteeswhen F is non-empty and satisfies the following assumption: Assumption 1.
The target set F defined in (20) is full-dimensional and there exists P center ∈ S ( k +1) n x such that { P ∈ S ( k +1) n x |k P − P center k F ≤ (cid:15) } ⊂ F where k·k F is the Frobenius norm. Cutting-plane methods [35, 36, 37] are iterative algorithms to find a point in a target convexset F or to determine whether F is empty. In these methods, we have no information on F except for a “cutting-plane oracle”, which can verify whether P ∈ F for a given P . Let ˜ F be alocalization set defined by a finite set of inequalities that over approximates the target set, i.e., F ⊆ ˜ F . If ˜ F is empty, then we have proof that the target set F is also empty. Otherwise,we query the oracle at a point P ∈ ˜ F . If P ∈ F , the oracle returns ‘yes’ and the algorithmterminates; if P / ∈ F , it returns ‘no’ together with a separating hyperplane that separates P and F . In the latter case, the cutting-plane method updates the localization set by ˜ F ← ˜ F ∩{ half space defined by the separating hyperplane } as shown in Fig. 1. This process continues untileither a point in the target set is found or the target set is certified to be empty.Based on how the query point is chosen, different cutting-plane methods have been proposedincluding the center of gravity method [38], the maximum volume ellipsoid (MVE) cutting-planemethod [20], the Chebyshev center cutting-plane method [36], the ellipsoid method [39, 40], and theanalytic center cutting-plane method [41, 42, 35]. In this paper, we use the analytic center cutting-plane method since it allows the localization set ˜ F to be described by linear matrix inequalities.In the ACCPM, the query point P is chosen as the analytic center of the localization set ˜ F . Definition 4 (Analytic center [43]) . The analytic center x ac of a set of convex inequalities andlinear equalities h i ( x ) ≤ , i = 1 , · · · , m, F x = g , is defined as the solution of the convex problem minimize x − m X i =1 log( − h i ( x )) subject to F x = g. (23)We design the learning strategy according to the ACCPM by constructing a learner, whichproposes Lyapunov function candidates based on a set of samples, and a verifier, which serves as acutting-plane oracle and updates the sample set with counterexamples.8 igure 1: In the i -th iteration, the query point P ( i ) is chosen as the “center” of the localization set ˜ F i toguarantee the removal of a portion of ˜ F i from the search space whenever a separating hyperplane (green solidline) is given. The localization set is then updated to ˜ F i +1 (shaded set) which is a finer over-approximationof the target set F . Let S = { x , x , · · · , x N } ⊂ X be a collection of samples from X . The localization set ˜ F in thespace of P is given by (22), which represents the learner’s knowledge about F by observing the k -step trajectories of the dynamical system (1) starting from S . According to the ACCPM, thelearner proposes a Lyapunov function candidate V ( k ) ( x ; P ac ) with P ac as the analytic center of ˜ F : P ac := argmin P − X x ∈S log( − ∆ V ( k ) ( x, P )) − log det( βI − P ) − log det( P − αI ) . (24)This is a convex program that can be solved efficiently through off-the-shell convex optimizationsolvers. We denote the objective function in (24) as φ ( ˜ F ) and call it the potential function on theset ˜ F . If (24) is infeasible, then we have a proof that no Lyapunov function exists in the functionclass of order k . Otherwise, the learner proposes V ( k ) ( x ; P ac ) = z ( k ) > P ac z ( k ) as a Lyapunov functioncandidate. Due to the log-barrier function in (24), P ac is in the interior of ˜ F . When the sampleset S is empty, the potential function simply becomes φ ( ˜ F ) = − log det( βI − P ) − log det( P − αI )with P ac = α + β I as the optimal solution. Suppose the learner proposes the Lyapunov function candidate V ( k ) ( x ; P ) by solving (24). Given V ( k ) ( x ; P ), the verifier either ensures that this function satisfies the constraints in (16a) and (17), orreturns a state where constraint (16a) or (17) is violated as a counterexample. Since the log-barrierfunction in (24) guarantees P (cid:31)
0, constraint (16a) is readily satisfied and the verifier must checkthe violation of constraint (17). This can be done by solving the optimization problemmaximize x ∈X \{ } ∆ V ( k ) ( x, P ) . (25)When f ( x ) has a mixed-integer formulation L gr( f ) = { ( x, x + , λ, µ ) ∈ R n x × R n x × R n λ × Z n µ | ‘ ( x, x + , λ, µ ) ≤ v } , V ( k ) ( x ; P ) of order k , we write problem (25) explicitly as amixed-integer quadratic program (MIQP):maximize { x i } , { λ i } , { µ i } x x ... x k +1 > P x x ... x k +1 − x x ... x k > P x x ... x k (26a)subject to x ∈ X (26b) k x k ∞ ≥ (cid:15) (26c) ‘ ( x i , x i +1 , λ i , µ i ) ≤ v i , i = 0 , , · · · , k (26d)where the objective function (26a) is equal to the Lyapunov difference ∆ V ( k ) ( x , P ), constraint (26b)restricts the search of counterexamples inside the ROI X , constraint (26c) excludes a small ‘ ∞ normball centered at the origin B (cid:15) = { x |k x k ∞ < (cid:15) } from the search space, and constraints (26d) enforcethe dynamical constraint ( x i , x i +1 ) ∈ gr( f ), or equivalently x i +1 = f ( x i ) for i = 0 , · · · , k .Denote p ∗ the optimal value and x ∗ the x -component of the optimal solution of (26). When p ∗ ≥ x ∗ is a counterexample for the Lyapunov function candidate V ( k ) ( x ; P ) since ∆ V ( k ) ( x ∗ , P ) ≥
0. Then, the separating hyperplane induced by the counterexample x ∗ is given as ∆ V ( k ) ( x ∗ , P ) = 0where ∆ V ( k ) ( x ∗ , P ) is interpreted as a linear function in the matrix variable P . When p ∗ < Corollary 1.
Let
Ω = { x | V ( k ) ( x ; P ) ≤ γ } be the largest sublevel set of V ( k ) ( x ; P ) inside X with γ =inf x ∈R\X V ( k ) ( x ; P ) . Define the successor set of B (cid:15) as suc ( B (cid:15) ) := { y | y = f ( x ) , x ∈ B (cid:15) } . Assume B (cid:15) ⊂ Ω , suc ( B (cid:15) ) ⊂ Ω , and Ω loc is the smallest sublevel set of V ( k ) ( x ; P ) such that suc ( B (cid:15) ) ⊆ Ω loc .If p ∗ < , then we have lim t →∞ x t ∈ Ω loc for all x ∈ Ω .Proof. First, note that all trajectories starting from x ∈ Ω and x / ∈ B (cid:15) reach B (cid:15) in a finite numberof steps since V ( k ) ( x t +1 ; P ) − V ( k ) ( x t ; P ) ≤ p ∗ < x t / ∈ B (cid:15) and V ( k ) ( x t ; P ) is lower-bounded by 0 for all x t . If x t never reaches B (cid:15) , we will have a contradiction that V ( k ) ( x N ; P ) < N . After the trajectories reach B (cid:15) , the subsequent states will remain inside Ω loc byconstruction. Therefore, we have lim t →∞ x t ∈ Ω loc for all x ∈ Ω.The alternation between the learner and the verifier is summarized in Algorithm 1. In the nextsubsection, we show that when the target set F satisfies Assumption 1, our proposed algorithm isguaranteed to find a feasible point in F in a finite number of steps. Remark 3.
In this paper, the verifier is constructed as a global optimization problem (25) . Alter-native constructions of the verifier do exist, e.g., by using SMT solvers as shown in [16, 17]. Thetermination results in the next subsection will always hold no matter how the verifier is formulated.
The convergence and complexity of the ACCPM have been studied in [35, 42, 44, 45, 46, 47] undervarious assumptions on the localization set, the form of the separating hyperplane, whether multiplecuts are applied, etc. Directly related to Algorithm 1 and the search for V ( k ) ( x ; P ) is [47] whichanalyzes the complexity of the ACCPM with a matrix variable and semidefiniteness constraints.Notably, it provides an upper bound on the number of iterations that Algorithm 1 can run beforetermination when the target set is non-empty. In [47], it is assumed that10 lgorithm 1 Learning-based Lyapunov function synthesis procedure LearningLyapunov S = ∅ . initialize the sample set i = 1 while True do Generate ˜ F i from sample set S i by (22). if ˜ F i = ∅ then Return: Status = Infeasible, P ∗ = ∅ . P ( i ) = arg min (24) with S i . find analytic center solve (26) with P ( i ) . query the verifier if max (26) < 0 then Return: Status = Feasible, P ∗ = P ( i ) . else x ( i ) ∗ = arg min x (26) . find a counterexample S i +1 = S i ∪ { x ( i ) ∗ } i = i + 1• A1: F is a convex subset of S n .• A2: See Assumption 1.• A3: F ⊂ { P ∈ S n | (cid:22) P (cid:22) I } .For the Lyapunov function candidate class of order k , we have n = ( k + 1) n x . Let the ACCPMstart with the localization set ˜ F = { P | (cid:22) P (cid:22) I } and initialize the first query point P (1) = I correspondingly. If at iteration i a query point P ( i ) is rejected by the oracle, a separating hyperplaneof the form h D i , P − P ( i ) i = 0 with k D i k F = 1 is given by the verifier. By induction, we have thatat iteration i >
1, the localization set ˜ F i is˜ F i = { P | (cid:22) P (cid:22) I, h D j , P i ≤ c j , j = 1 , · · · , i − } . with D j defining the separating hyperplane and c j = h D j , P ( j ) i . The query point at iteration i isgiven by P ( i ) = arg min φ ( ˜ F i ) with the potential function φ ( ˜ F i ) = − i − X j =1 log( c j − h D j , P i ) − log det( I − P ) − log det( P ) . Theorem 2.
Under Assumption 1 on the target set F , Algorithm 1 finds a feasible point in F inat most O ((( k + 1) n x ) /(cid:15) ) iterations.Proof. The proof follows from [47] which states for a general target set
F ⊂ S n under assump-tions A1 to A3, the analytic center cutting-plane method with separating hyperplanes of the form h D i , P i ≤ c i is shown to find a feasible point in at most O ( n /(cid:15) ) iterations. For the sequenceof localization sets ˜ F i , [47] computes an upper bound on the potential function φ ( ˜ F i ) which isapproximately i log( (cid:15) ), and a lower bound on φ ( ˜ F i ) which is proportional to i log( in ). Since theACCPM must terminate before the lower bound exceeds the upper bound of φ ( ˜ F i ), we obtain the O ( n /(cid:15) ) upper bound on the number of iterations.To show that the result in [47] applies to Algorithm 1, note that for each counterexample x ( i ) ∗ found by the verifier in iteration i , the separating hyperplane can be constructed as ∆ V ( k ) ( x ( i ) ∗ , P ) ≤ V ( k ) ( x ( i ) ∗ , P ( i ) ) since ∆ V ( k ) ( x ( i ) ∗ , P ( i ) ) ≥
0. We can rewrite this cutting plane in the form h D i , P i ≤h D i , P ( i ) i by setting ˆ D i = z + , ( i ) k, ∗ z + , ( i ) , > k, ∗ − z ( i ) k, ∗ z ( i ) , > k, ∗ , D i = ˆ D i / k ˆ D i k F , where z ( i ) k, ∗ denotes the basis (19)with x = x ( i ) ∗ , and z + , ( i ) k, ∗ denotes the basis (19) after setting x = f ( x ( i ) ∗ ). Then with α = 0 , β = 1in (20), Algorithm 1 satisfies assumptions A1 to A3 and it terminates in at most O ((( k + 1) n x ) /(cid:15) )iterations according to [47].Theorem 2 provides a finite-step termination guarantee for Algorithm 1 when the target set satis-fies Assumption 1. However, when F is empty, we do not have such a guarantee. Certification of thenon-existence of Lyapunov functions in V ( k ) ( x ; P ) relies on detecting that an over-approximation˜ F is empty. Hence, a quick expansion of the sample set S as shown in Section 4.5.5 is preferred. For a sample set S with N samples, the localization set ˜ F in (22) is described by N linear aswell as two semidefinite constraints on the variable P ∈ S ( k +1) n x . We first decide if ˜ F is emptyby solving an SDP feasibility problem. If the SDP is infeasible, then ˜ F = ∅ and so is F . If F isnon-empty, we move on to the analytic center problem (24) which can be solved, e.g., through aninfeasible start Newton’s method [43]. Since MIQP is well-known to be NP-hard, we use the number of binary variables, which we denoteby N µ in (26), as a rough measure of the complexity of (26). When f ( x ) is a PWA functionwith N mode modes, we have N mode binary variables in the MIL formulation of f ( x ). For the LCsystem (8), N µ equals the dimension of the orthogonal variables v and w . N µ is explicitly givenin the MLD system and is equal to the number of neurons in the mixed-integer formulation of theReLU networks. It follows that for the LC systems and ReLU networks, it is possible to use a smallnumber of integer variables to encode a PWA system with many more modes. Since we need toevaluate the hybrid dynamics f ( x ) for k times when V ( k ) ( x ; P ) is applied, the number of binaryvariables N µ in (26) is largely linear in the order k .The actual solving time of the MIQP has a complex dependence not only on the number ofvariables and constraints but also on how the constraints are formulated. The exploration of thenumerical performance of the proposed algorithm is left for future research. The optimization problem (26) is a nonconvex MIQP since the quadratic objective function is in-definite. Therefore, the relaxation of the problem after removing the integrality constraints wouldresult in a nonconvex quadratic program. Nonconvex MIQP can be solved to global optimalitythrough Gurobi v9 . .5.4 Exclusion of the origin In constraint (26c), we add a guard B (cid:15) at the origin to approximate the exact constraint x ∈X \{ } . Since constraint (26c) allows a mixed-integer linear formulation through the big-M method,problem (26) is an MIQP. Adding B (cid:15) bounds p ∗ off from 0 if V k ( x ; P ) is in fact a Lyapunov functionand we can decide the negativity of p ∗ by checking if p ∗ < − (cid:15) tol for some tolerance (cid:15) tol > x + = f ( x ) is linear inside B (cid:15) , i.e., x + = Ax for x ∈ B (cid:15) , we can showconvergence to the origin of the system trajectories starting inside B (cid:15) by checking the magnitudeof the eigenvalues of A . Combined with Corollary 1, asymptotic convergence inside the sublevel setΩ can be established. We can terminate the Branch & Bound [52] algorithm in solving (26) once a feasible solution isfound with non-negative objective function since in this case we already have a counterexample.
We demonstrate our method through two examples: closed-loop MPC systems and ReLU neuralnetwork controlled PWA systems. In particular, we compare the performances of our method withthe LMI-based Lyapunov function synthesis methods [7] on the MPC example. Algorithm 1 isimplemented in Python 3 . α = 0 . , β = 1 . (cid:15) tol = 10 − to decide negativity of the optimal value of the MIQP. The Gurobi solver is setwith feasibility tolerance 10 − , integer tolerance 10 − , and optimality tolerance 10 − . In addition,we terminate the MIQP once it finds a feasible solution that generates an objective ≥ − whichmeans a counterexample is already found. For an unstable linear system x + = Ax + Bu with A = " . .
20 1 . , B = " . , (27)we design an MPC controller with horizon T = 10, state constraint [ − − > ≤ x ≤ [5 5] > ,control input constraint − ≤ u ≤
1, stage cost p ( x, u ) = x > Qx + u > Ru with Q = 10 I, R = 1,and terminal cost q ( x ) = x > P ∞ x where P ∞ = DARE( A, B, Q, R ) is the solution to the discretealgebraic Riccati equation defined by (
A, B, Q, R ). The terminal set is chosen as the maximumpositive invariant set [2, Chapter 10] of the closed-loop system x + = ( A + BK ∞ ) x where K ∞ = − ( B > P ∞ B + R ) − B > P ∞ A .With the above setup, we obtain the explicit MPC controller, which is a PWA function of thestate x with 211 partitions, through the MPT3 toolbox in Matlab. The domain R of the explicitMPC is a polytope shown in Fig. 2 together with its partitions. After obtaining R , we verify thatit is positive invariant for the closed-loop dynamics. Through the LMI-based method, we are ableto synthesize a discontinuous PWQ Lyapunov function in MPT3 which verifies R is the ROA for13 igure 2: Domain R and the 211 modes of the closed-loop MPC system of (27). the closed-loop MPC system. The total running time is 38 .
890 seconds, with 0 .
25 spent in solvingthe constructed SDP and the rest in computing the transition map.
We import the PWA representation of the closed-loop MPC dynamics which we denote as x + = f cl ( x ) from MPT3 and run Algorithm 1 with the mixed-integer formulation of f cl ( x ) as shown inSection 2.1. Since the domain R is positively invariant, we set X = R . For Lyapunov functioncandidates of order k = 0 ,
1, Algorithm 1 certifies the non-existence of Lyapunov functions after1 .
709 seconds with 4 iterations and 28 .
333 seconds with 7 iterations, respectively. With V ( k ) ( x ; P )of order k = 2, Algorithm 1 terminates in 14 iterations with a valid Lyapunov function candidateand the total running time is 4379 .
716 seconds. We plot the counterexamples found in each iterationin Fig. 3a. The accumulated running time of Algorithm 1 in each iteration is shown in Fig. 3b. −4 −2 0 2 4x −4−2024 x = (a) Sequence of counterexamples found in Algorithm 1. A cc u m u l a t e d s o l v e r t i m e / s e c o n d s PWALCS (b) Accumulated running time of Algorithm 1.
Figure 3: (a) Sequence of counterexample states found by the verifier in Algorithm 1 with the PWA represen-tation of the closed-loop MPC system. (b) Accumulated running time of Algorithm 1 in each iteration withmixed-integer formulations induced by the PWA (blue) and LC (orange) representations of the closed-loopMPC system. .1.2 Algorithm 1 with LCS representation As shown in Section 2.2 and [26], we can obtain a mixed-integer formulation of the closed-loopMPC dynamics through its LC system representation instead of the PWA one. Based on thismixed-integer formulation, we run Algorithm 1 with Lyapunov function candidates V ( k ) ( x ; P ) oforder k = 2. The algorithm terminates in 368 .
618 seconds with 20 iterations and returns a validLyapunov function which certifies that the domain R is the ROA. The accumulated running timeat each iteration is summarized in Fig. 3b. Fig. 3b shows that Algorithm 1 depends on the mixed-integer representation of the system. Inthis example, we need 211 binary variables to describe the map x + = f cl ( x ) using the PWA repre-sentation, while with the LCS representation, we only need to use 64 binary variables to describethe map u = π MP C ( x ), and hence the closed-loop dynamics. However, when compared with theLMI-based method which synthesizes a PWQ discontinuous Lyapunov functions in 38 .
890 seconds,Algorithm 1 is rather inefficient. This is partly due to that the continuous Lyapunov functioncandidate class V ( k ) ( x ; P ) is more conservative than the discontinuous one [7]. To compensate forthe conservatism, we need to apply high order V ( k ) ( x ; P ) which increases the complexity of theMIQP. In the next subsection, we show that Algorithm 1 can synthesize a Lyapunov function whenthe LMI-based method fails. Consider model predictive control of a randomly generated 4-dimensional open-loop unstable linearsystem given by A = . − . − . . − . . − . . − . . − . . . . . . , B = − . − . . . With the state constraint k x k ∞ ≤ − ≤ u ≤
1, we design the MPCcontroller by choosing horizon T = 10, stage cost with Q = 10 I, R = 1, terminal cost with P ∞ andthe terminal set as in Section 5.1. The explicit MPC controller is constructed through MPT3 andhas 193 partitions in its domain R , which is validated to be positive invariant under the closed-loopdynamics. Although it took only 29 . V ( k ) ( x ; P ) with k = 1 and theLCS representation of the closed-loop MPC system. The ROI is chosen as X = R . Algorithm 1terminates in 9 iterations with a valid Lyapunov function candidate and therefore proves that theclosed-loop system is asymptotically stable in the domain R . The total running time is 680 . We use a hybrid system example from [55] and consider the inverted pendulum shown in Fig. 5 withparameters m = 1 , ‘ = 1 , g = 10 , k = 100 , d = 0 .
1. We denote the angle and the angular velocityof the pendulum by q and ˙ q , respectively, and define the system state as x = ( q, ˙ q ). By linearizingthe dynamics of the inverted pendulum around x = 0, we obtain a hybrid system which has two15 A cc u m u l a t e d s o l v e r t i m e / s e c o n d s Figure 4: Accumulated solver time of Algorithm 1 in each iteration for the 4-dimensional closed-loop MPCsystem. modes: not in contact with the elastic wall (mode 1) and in contact with the elastic wall (mode2). After discretizing the model using the explicit Euler scheme with a sampling time h = 0 .
01, aPWA model x + = ψ ( x, u ) of the form (4) is obtained with the following parameters A = " . . , B = " . , c = " , R = { x | [ − . − . > ≤ x ≤ [0 . . > } .A = " . − . , B = " . , c = " . , R = { x | [0 . − . > ≤ x ≤ [0 . . > } . (28) q l m g k d u x −0.20−0.15−0.10−0.050.000.050.100.150.20 x −1.5−1.0−0.50.0 0.5 1.0 1.5u −4−2024 Figure 5: A ReLU neural network (right) that approximates a hybrid MPC controller is applied on theinverted pendulum system with an elastic wall (left).
We then synthesize a hybrid MPC controller π MP C ( x ) for the PWA system [22] where thecontrol input constraints are given by − ≤ u ≤ T = 10. The16 −1.5−1.0−0.50.00.51.01.5 x V (0) (x; P) =0.86 estimate of ROA (a) Estimate of ROA of the neural network controlledsystem by a quadratic Lyapunov function in V (0) ( x ; P ). −0.2 −0.1 0.0 0.1 0.2x −1.5−1.0−0.50.00.51.01.5 x V (1) (x; P) =1.0 estimate of ROA (b) Estimate of ROA of the neural network controlledsystem by a PWQ Lyapunov function in V (1) ( x ; P ). Figure 6: Estimates of ROA found by Algorithm 1 with Lyapunov function candidates in V (0) ( x ; P ) and V (1) ( x ; P ). The partitions of the state space are marked by the black boxes. Simulated closed-loop trajec-tories with the neural network controller are plotted for a grid of initial conditions. stage and terminal costs are given by Q = I, R = 1, and P ∞ = DARE( A , B , Q, R ). We evaluate π MP C ( x ) on a uniform 40 ×
40 grid samples from the state space and let the reference ROI X bethe convex hull of all the feasible state samples. Then the ROI X = γ X with 0 < γ ≤ x, π MP C ( x )) aregenerated to train a ReLU neural network π ( x ) in Keras [56] to approximate the MPC controller.The neural network has 2 hidden layers with 20 neurons in each layer and its output layer biasterm is modified after training to guarantee π (0) = 0. We plot the neural network controller inFig. 5 and set (cid:15) = 0 . B (cid:15) .For Lyapunov function candidates of order k = 0 and k = 1, we run Algorithm 1 with ROI X = γ X of varying values of γ . For the quadratic function class V (0) ( x ; P ), the largest ROI isgiven by X = 0 . X through bisection with which Algorithm 1 terminates in 16 iterations witha total running time of 13 .
687 seconds. The corresponding estimate of ROA is shown in Fig. 6a.For the PWQ function class V (1) ( x ; P ), the largest ROI is given by X = 1 . X in which caseAlgorithm 1 terminates in 11 iterations with a total running time of 90 .
917 seconds. The estimateof ROA obtained by the found PWQ Lyapunov function candidate is shown in Fig. 6b. It showsthat the synthesized Lyapunov function in V (1) ( x ; P ) is less conservative compared with the onein V (0) ( x ; P ) and Algorithm 1 can obtain non-trivial estimates of ROA for the neural networkcontrolled hybrid systems. We have proposed a learning-based method to learn Lyapunov functions for autonomous hybridsystems that have a mixed-integer formulation, including piecewise affine, linear complementarity,mixed logical dynamical systems and ReLU neural networks. By designing the method accordingto the analytic center cutting-plane method, we show that the proposed algorithm is guaranteedto find a Lyapunov function in a finite number of steps when the set of Lyapunov functions is17ull-dimensional in the parameter space. Our method is an alternative to the LMI-based Lyapunovfunction synthesis approach which relies on the piecewise affine representation of hybrid systems.
References [1] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,”
Auto-matica , vol. 35, no. 3, pp. 407–427, 1999.[2] F. Borrelli, A. Bemporad, and M. Morari,
Predictive control for linear and hybrid systems . CambridgeUniversity Press, 2017.[3] W. Heemels, J. M. Schumacher, and S. Weiland, “Linear complementarity systems,”
SIAM journal onapplied mathematics , vol. 60, no. 4, pp. 1234–1269, 2000.[4] A. J. van der Schaft and J. M. Schumacher, “Complementarity modeling of hybrid systems,”
IEEETransactions on Automatic Control , vol. 43, no. 4, pp. 483–490, 1998.[5] E. Sontag, “Nonlinear regulation: The piecewise linear approach,”
IEEE Transactions on automaticcontrol , vol. 26, no. 2, pp. 346–358, 1981.[6] W. P. Heemels, B. De Schutter, and A. Bemporad, “Equivalence of hybrid dynamical models,”
Auto-matica , vol. 37, no. 7, pp. 1085–1091, 2001.[7] P. Biswas, P. Grieder, J. Löfberg, and M. Morari, “A survey on stability analysis of discrete-timepiecewise affine systems,”
IFAC Proceedings Volumes , vol. 38, no. 1, pp. 283–294, 2005.[8] H. Lin and P. J. Antsaklis, “Stability and stabilizability of switched linear systems: a survey of recentresults,”
IEEE Transactions on Automatic control , vol. 54, no. 2, pp. 308–322, 2009.[9] Z. Sun, “Stability of piecewise linear systems revisited,”
Annual Reviews in Control , vol. 34, no. 2,pp. 221–231, 2010.[10] R. Pascanu, G. Montufar, and Y. Bengio, “On the number of response regions of deep feed forwardnetworks with piece-wise linear activations,” arXiv preprint arXiv:1312.6098 , 2013.[11] M. K.-J. Johansson,
Piecewise linear control systems: a computational approach , vol. 284. Springer,2003.[12] M. Johansson and A. Rantzer, “Computation of piecewise quadratic lyapunov functions for hybridsystems,” in , pp. 2005–2010, IEEE, 1997.[13] S. Prajna and A. Papachristodoulou, “Analysis of switched and hybrid systems-beyond piecewisequadratic methods,” in
Proceedings of the 2003 American Control Conference, 2003. , vol. 4, pp. 2779–2784, IEEE, 2003.[14] A. Solar-Lezama, L. Tancau, R. Bodik, S. Seshia, and V. Saraswat, “Combinatorial sketching for finiteprograms,” in
Proceedings of the 12th international conference on Architectural support for programminglanguages and operating systems , pp. 404–415, 2006.[15] A. Solar-Lezama and R. Bodik,
Program synthesis by sketching . Citeseer, 2008.[16] D. Ahmed, A. Peruffo, and A. Abate, “Automated and sound synthesis of lyapunov functions with smtsolvers,” in
International Conference on Tools and Algorithms for the Construction and Analysis ofSystems , pp. 97–114, Springer, 2020.[17] A. Abate, D. Ahmed, M. Giacobbe, and A. Peruffo, “Formal synthesis of lyapunov neural networks,”
IEEE Control Systems Letters , 2020.[18] J. Kapinski, J. V. Deshmukh, S. Sankaranarayanan, and N. Arechiga, “Simulation-guided lyapunovanalysis for hybrid dynamical systems,” in
Proceedings of the 17th international conference on Hybridsystems: computation and control , pp. 133–142, 2014.
19] H. Ravanbakhsh and S. Sankaranarayanan, “Learning control lyapunov functions from counterexamplesand demonstrations,”
Autonomous Robots , vol. 43, no. 2, pp. 275–307, 2019.[20] I. E. S. Tarasov, L. G. Khachiyan, “The method of inscribed ellipsoids,”
Soviet Mathematics Doklady ,vol. 37, 1988.[21] L. Lindemann, H. Hu, A. Robey, H. Zhang, D. Dimarogonas, S. Tu, and N. Matni, “Learning hybridcontrol barrier functions from data,” in
Conference on Robot Learning , 2020.[22] T. Marcucci and R. Tedrake, “Mixed-integer formulations for optimal control of piecewise-affine sys-tems,” in
Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation andControl , pp. 230–239, 2019.[23] J. P. Vielma, “Mixed integer linear programming formulation techniques,”
Siam Review , vol. 57, no. 1,pp. 3–57, 2015.[24] T. Kleinert, M. Labbé, F. Plein, and M. Schmidt, “Technical note—there’s no free lunch: On thehardness of choosing a correct big-m in bilevel optimization,”
Operations Research .[25] L. Gurobi Optimization, “Gurobi optimizer reference manual,” 2020.[26] D. Simon and J. Löfberg, “Stability analysis of model predictive controllers using mixed integer linearprogramming,” in , pp. 7270–7275, IEEE,2016.[27] V. Tjeng, K. Xiao, and R. Tedrake, “Evaluating robustness of neural networks with mixed integerprogramming,” arXiv preprint arXiv:1711.07356 , 2017.[28] T.-W. Weng, H. Zhang, H. Chen, Z. Song, C.-J. Hsieh, D. Boning, I. S. Dhillon, and L. Daniel, “Towardsfast computation of certified robustness for relu networks,” arXiv preprint arXiv:1804.09699 , 2018.[29] M. Hein and M. Andriushchenko, “Formal guarantees on the robustness of a classifier against adversarialmanipulation,” in
Advances in Neural Information Processing Systems , pp. 2266–2276, 2017.[30] E. Wong and Z. Kolter, “Provable defenses against adversarial examples via the convex outer adversarialpolytope,” in
International Conference on Machine Learning , pp. 5286–5295, 2018.[31] W. M. Haddad and V. Chellaboina,
Nonlinear dynamical systems and control: a Lyapunov-based ap-proach . Princeton university press, 2011.[32] D. Aeyels and J. Peuteman, “A new asymptotic stability criterion for nonlinear time-variant differentialequations,”
IEEE Transactions on automatic control , vol. 43, no. 7, pp. 968–971, 1998.[33] R. Bobiti and M. Lazar, “A sampling approach to finding lyapunov functions for nonlinear discrete-timesystems,” in , pp. 561–566, IEEE, 2016.[34] A. A. Ahmadi and P. A. Parrilo, “Non-monotonic lyapunov functions for stability of discrete timenonlinear and switched systems,” in , pp. 614–621,IEEE, 2008.[35] D. S. Atkinson and P. M. Vaidya, “A cutting plane algorithm for convex programming that uses analyticcenters,”
Mathematical Programming , vol. 69, no. 1-3, pp. 1–43, 1995.[36] J. Elzinga and T. G. Moore, “A central cutting plane algorithm for the convex programming problem,”
Mathematical Programming , vol. 8, no. 1, pp. 134–145, 1975.[37] S. Boyd and L. Vandenberghe, “Localization and cutting-plane methods,”
From Stanford EE 364blecture notes , 2007.[38] A. Y. Levin, “An algorithm for minimizing convex functions,”
Soviet Mathematics Doklady , vol. 160,pp. 1244–1247, 1965.[39] L. G. Khachiyan, “Polynomial algorithms in linear programming,”
USSR Computational Mathematicsand Mathematical Physics , vol. 20, no. 1, pp. 53–72, 1980.
40] D. Y. A.S. Nemirovskii, “Problem complexity and method efficiency in optimization,” 1983.[41] J.-L. Goffin and J.-P. Vial, “On the computation of weighted analytic centers and dual ellipsoids withthe projective algorithm,”
Mathematical Programming , vol. 60, no. 1-3, pp. 81–92, 1993.[42] Y. Nesterov, “Cutting plane algorithms from analytic centers: efficiency estimates,”
Mathematical Pro-gramming , vol. 69, no. 1, pp. 149–176, 1995.[43] S. Boyd, S. P. Boyd, and L. Vandenberghe,
Convex optimization . Cambridge university press, 2004.[44] Y. Ye, “A potential reduction algorithm allowing column generation,”
SIAM Journal on Optimization ,vol. 2, no. 1, pp. 7–20, 1992.[45] Z.-Q. Luo and J. Sun, “A polynomial cutting surfaces algorithm for the convex feasibility problemdefined by self-concordant inequalities,”
Computational Optimization and Applications , vol. 15, no. 2,pp. 167–191, 2000.[46] J.-L. Goffin, Z.-Q. Luo, and Y. Ye, “Complexity analysis of an interior cutting plane method for convexfeasibility problems,”
SIAM Journal on Optimization , vol. 6, no. 3, pp. 638–652, 1996.[47] J. Sun, K.-C. Toh, and G. Zhao, “An analytic center cutting plane method for semidefinite feasibilityproblems,”
Mathematics of Operations Research , vol. 27, no. 2, pp. 332–346, 2002.[48] P. Belotti, C. Kirches, S. Leyffer, J. Linderoth, J. Luedtke, and A. Mahajan, “Mixed-integer nonlinearoptimization,”
Acta Numerica , vol. 22, p. 1, 2013.[49] S. Vigerske, “Decomposition in multistage stochastic programming and a constraint integer program-ming approach to mixed-integer nonlinear programming,” 2013.[50] M. Tawarmalani and N. V. Sahinidis,
Convexification and global optimization in continuous and mixed-integer nonlinear programming: theory, algorithms, software, and applications , vol. 65. Springer Science& Business Media, 2013.[51] P. Belotti, J. Lee, L. Liberti, F. Margot, and A. Wächter, “Branching and bounds tighteningtechniquesfor non-convex minlp,”
Optimization Methods & Software , vol. 24, no. 4-5, pp. 597–634, 2009.[52] L. A. Wolsey and G. L. Nemhauser,
Integer and combinatorial optimization , vol. 55. John Wiley &Sons, 1999.[53] M. Herceg, M. Kvasnica, C. Jones, and M. Morari, “Multi-Parametric Toolbox 3.0,” in
Proc. of theEuropean Control Conference , (Zürich, Switzerland), pp. 502–510, July 17–19 2013. http://control.ee.ethz.ch/~mpt .[54] M. ApS,
The MOSEK optimization toolbox for Python manual. Version 9.2 , 2020.[55] T. Marcucci, R. Deits, M. Gabiccini, A. Bicchi, and R. Tedrake, “Approximate hybrid model predictivecontrol for multi-contact push recovery in complex environments,” in , pp. 31–38, IEEE, 2017.[56] F. Chollet et al. , “Keras,” 2015. https://github.com/fchollet/keras ..