Synthesis of Hybrid Automata with Affine Dynamics from Time-Series Data
Miriam GarcΓa Soto, Thomas A. Henzinger, Christian Schilling
SSynthesis of Hybrid Automata with Affine Dynamics fromTime-Series Data
Miriam GarcΓa Soto
IST AustriaKlosterneuburg, [email protected]
Thomas A. Henzinger
IST AustriaKlosterneuburg, [email protected]
Christian Schilling
University of KonstanzKonstanz, [email protected]
ABSTRACT
Formal design of embedded and cyber-physical systems relies onmathematical modeling. In this paper, we consider the model classof hybrid automata whose dynamics are defined by affine differ-ential equations. Given a set of time-series data, we present analgorithmic approach to synthesize a hybrid automaton exhibitingbehavior that is close to the data, up to a specified precision, andchanges in synchrony with the data. A fundamental problem in oursynthesis algorithm is to check membership of a time series in ahybrid automaton. Our solution integrates reachability and opti-mization techniques for affine dynamical systems to obtain both asufficient and a necessary condition for membership, combined in arefinement framework. The algorithm processes one time series ata time and hence can be interrupted, provide an intermediate result,and be resumed. We report experimental results demonstrating theapplicability of our synthesis approach.
KEYWORDS synthesis, hybrid automaton, linear dynamics, membership
Formal design and verification of embedded control systems requirea mathematical model capturing the dynamics of each componentin the system. In general, embedded systems combine analog anddigital components. The analog components evolve continuouslyin real time, while the digital components evolve in discrete time.An appropriate mathematical formalism for modeling systems withmixed continuous and discrete behavior is a hybrid automaton [16].In this paper we propose an automated approach to synthesizinga hybrid automaton with affine continuous dynamics (abbreviatedadha) from time-series data in an online fashion. The design ofmodels from observed data has been extensively studied in controltheory for autoregressive systems [4, 6, 11, 28, 34], which can beseen as discrete dynamical systems, in contrast to the continuous dy-namics captured by a hybrid automaton. Most of these approachesprocess a single time-series or all data at once. In a setting where notall data is available at once, it is desirable to have an online approachthat processes time-series data sequentially and iteratively updatesa model; only a few approaches support this feature [15, 32, 33, 35].Our synthesis approach operates in two phases. In the first phasewe transform a (discrete) time-series into a piecewise continuoustrajectory π , for which we present an optimization procedure thatallows to specify the error between the data and the trajectory.The trajectories π we consider are piecewise-affine (pwa) functionswhere each piece is the solution of an affine dynamical system ofthe form (cid:164) x = π΄ x + b . pwa trajectories can model a large class ofphysical processes and approximate generic nonlinear systems. In the second phase, which is independent of how the continuous pwa trajectory π has been obtained, we synthesize an adha from π .More precisely, we construct an adha from an existing adha (ini-tialized with the βemptyβ adha) in two stages: 1) (membership) wedetermine whether the new trajectory is already captured by anexecution of the model, up to a predefined precision, and 2) (modelupdate) if the trajectory is not captured, we modify the model suchthat, after the modification, the new model captures the trajectory(and all trajectories that had been captured before).We propose a three-step algorithm for the membership problem(βis a pwa trajectory captured by an adha?β). The first step is areachability analysis inside a tube around the trajectory that we useto provide a negative answer. This problem has been studied in [33]for the class of hybrid automata with piecewise-constant dynamics.The second step is an optimization-based analysis that we use toprovide a positive answer. The third step is a refinement procedureto deal with cases when the first two steps were not conclusive.If we find that the pwa trajectory is not captured by the model inthe membership query, we apply a model update by adding behaviorto the automaton. We first try to relax the continuous constraintsof the automaton (called invariants and guards). If this relaxationis not sufficient to capture the trajectory, we also apply structuralchanges to the automaton (adding transitions and locations).In summary, we present algorithms to solve the following prob-lems for pwa trajectories and adhas with a given precision: β’ transforming time-series data to pwa trajectories (Section 4) β’ membership of a pwa trajectory in an adha (Section 5) β’ synthesizing an adha from pwa trajectories (Section 6)Together, our algorithms form an end-to-end approach to the syn-thesis of an adha from time-series data with a given precision. Related work.
The synthesis of hybrid systems has been ex-plored previously in different fields and is known as identifica-tion in the area of control theory (see the surveys [11, 28]) and as process mining and model learning to a broader research commu-nity. Most of the techniques focus on input-output models, such asswitched autoregressive exogenous (SARX) [15, 25] and (PWARX)models [6, 9, 10, 17, 23, 31]. SARX models constitute a subclassof linear hybrid automata (which, unlike the adha, only has dy-namics with constant derivatives) with deterministic switchingbehavior and PWARX models are piecewise ARX models where theregressor space forms a state-space polyhedral partition. The afore-mentioned methods mainly consider single-input single-output(SISO) systems, whereas a few of them consider multiple-inputmultiple-output (MIMO) systems [4, 18, 34]. Other techniques iden-tify piecewise affine systems in state-space form [2, 22, 34]. Theidentification techniques can also be classified into optimization-based methods [19, 26] clustering-based procedures [9, 23] and a r X i v : . [ ee ss . S Y ] F e b iriam GarcΓa Soto, Thomas A. Henzinger, and Christian Schilling algebraic approaches [4, 24]. Most of these methods are proposedfor offline identification, with some exceptions [15, 32, 35]. We pro-pose an online approach that synthesizes hybrid automata withaffine dynamics, which are systems in state-space form.In the field of computer science, we find techniques for learningmodels from traces, which refers to approaches based on learningfinite-state machines [3] or other machine-learning techniques.Most approaches learn a (simpler) linear hybrid automaton. Thework in [20] describes an abstract framework, based on heuristics, tolearn offline from input-output traces by first learning the discretestructure and later adding continuous dynamics. Bartocci et al. learn shape expressions , which have a similar expressiveness [5]. A recentonline approach provides soundness and precision guarantees [33].However, that approach is restricted to linear hybrid automata,i.e., constant dynamics. We consider affine dynamics and follow aprincipled search algorithm for the automaton modification.We are not aware of approaches that transform time-series datato continuous affine dynamical systems. Some approaches considerdiscrete-time models, such as the work by Willems for LTI sys-tems [36], and other approaches for SARX models based on convexoptimization [27] or generalized principal component analysis [4]. Sets.
Let R , R β₯ , and N denote the set of real numbers, non-negativereal numbers, and natural numbers, respectively. Given a set π , the power set P ( π ) is the set of all subsets of π . We write x for points ( π₯ , . . . , π₯ π ) in R π . Given a point x β R π and π β R β₯ , we define the ball of radius π around x as π΅ π ( x ) : = { y β R π : β₯ x β y β₯ β€ π } , where β₯ Β· β₯ is the infinity norm. Given two sets π, π β² β R π , we definethe distance between π and π β² as π ( π, π β² ) : = inf {β₯ x β y β₯ : x β π, y β π β² } . Let a β R π and π β R be constant and x be a variable in R π , and let β¨ a , x β© denote the dot product of a and x ; then β¨ a , x β© βΌ π is a linear constraint where βΌ β { = , β€} , the set { x : β¨ a , x β© = π } is a hyperplane , and the set { x : β¨ a , x β© β€ π } is a half-space . A (convex)polytope is a compact intersection of linear constraints. Equivalently,a polytope is the convex hull of a set of vertices v , . . . , v π β R π ,written chull ({ v , . . . , v π }) . For a polytope π we denote the setof its linear constraints by constr ( π ) and the set of its vertices by vert ( π ) . Let cpoly ( π ) be the set of convex polytopes over R π . Trees. A tree is a directed acyclic graph T = ( N , E ) with finite setof nodes N , including a root node, and edges E β N Γ N . Given a node π β N , the child nodes are children ( π ) = { π β² β N : ( π, π β² ) β E } . Functions, dynamical systems, and trajectories.
Given a function π , let dom ( π ) denote its domain. Let π β π· denote the restrictionof π to domain π· β dom ( π ) . Given two functions π and π with dom ( π ) = dom ( π ) , the distance between π and π is denoted by π ( π , π ) and defined as max π‘ β dom ( π ) β₯ π ( π‘ ) β π ( π‘ )β₯ . We typically have dom ( π ) = [ ,π ] , where the initial and final states of π correspondto π ( ) and π ( π ) and are denoted by π and π end , respectively. A time series is a sampling π : π· β R π over a finite time domain π· .A function π : [ ,π ] β R π is a piecewise-affine (pwa) trajectory with π pieces if it is continuous and there is a tuple (I , A , B) where I is a finite set of consecutive time intervals [ π‘ , π‘ ] , . . . , [ π‘ π β , π‘ π ] with [ ,π ] = βͺ β€ π β€ π [ π‘ π β , π‘ π ] , A and B are π -tuples of matrices π΄ π β R π Γ π and vectors b π β R π , respectively, π = , . . . , π , and π β [ π‘ π β ,π‘ π ] is a solution of the affine dynamical system (cid:164) x = π΄ π x + b π , where (cid:164) x denotes the derivative of x with respect to π‘ . We assumethat pwa trajectories are given as the above tuple. We call π β [ π‘ π β ,π‘ π ] the pieces of π , and t s ( π ) : = { π‘ , . . . , π‘ π } the switching times of π .Each piece of π is called an affine trajectory. A linear trajectory π isa special case of an affine trajectory where b = . We consider a particular class of hybrid automata [16] with invari-ants and guards given by linear constraints and with continuousdynamics given by affine differential equations.
Definition 2.1. An π -dimensional hybrid automaton with affinedynamics ( adha ) is a tuple H = ( Q , E , X , Flow , Inv , Grd ) , where1) Q is a finite set of locations, 2) E β Q Γ Q is a transition relation,3) X = R π is the continuous state space, 4) Flow : Q β R π Γ π Γ R π is the injective flow function that returns a matrix π΄ and a vector b ,and we write Flow π΄ ( π ) β R π Γ π and Flow b ( π ) β R π to refer to eachcomponent, 5) Inv : Q β cpoly ( R π ) is the invariant function, and6) Grd : E β cpoly ( R π ) is the guard function.A path π in H of length π is a sequence of locations π , . . . , π π in Q such that ( π π , π π + ) β E for each 1 β€ π < π . We write paths (H) for the set of paths in H . Given a path π = π , . . . , π π , we define len ( π ) = π as the length of π and last ( π ) = π π as the last location.Next we define an execution of an adha, describing the evolutionof the continuous state subject to time passing and discrete switches. Definition 2.2. An execution π of an adha H is a pwa trajectory π : πΌ β R π such that there is a path π = π , . . . , π π in H and asequence of time points π‘ = , π‘ , . . . , π‘ π satisfying 1) πΌ = [ π‘ , π‘ π ] ,2) π ( π‘ ) β Inv ( π π ) for every 1 β€ π β€ π and π‘ β [ π‘ π β , π‘ π ] , 3) π ( π‘ π ) β Grd ( π π , π π + ) for every 1 β€ π < π , and 4) (cid:164) π ( π‘ ) = Flow π΄ ( π π ) Β· π ( π‘ )+ Flow b ( π π ) for every 1 β€ π β€ π and π‘ β ( π‘ π β , π‘ π ) .Thus switches between dynamics are state-dependent. We call t s ( π ) = { π‘ , . . . , π‘ π } the switching times of π . We say that π follows π , written π (cid:123) π , and denote the set of executions by exec (H) . Our overall goal is to synthesize a hybrid automaton from data,given in the form of time series, such that the synthesized automa-ton captures the dynamical behavior of the data up to a givenprecision. We split up this problem into two phases. In the firstphase, given a time series π and a value πΏ β R β₯ , we find a pwatrajectory π that is πΏ -close to all points in π . Definition 3.1.
Given a time series π with domain π· β [ ,π ] , apwa trajectory π with dom ( π ) = [ ,π ] , and a value πΏ β R β₯ , wesay that π πΏ -captures π if β₯ π ( π‘ ) β π ( π‘ )β₯ β€ πΏ for each π‘ β π· .In the second phase, given another value π β R β₯ , we constructa hybrid automaton from this pwa trajectory. Definition 3.2.
Given a pwa trajectory π and a value π β R β₯ ,we say that an adha H π -captures π if there exists an execution π β exec (H) such that π ( π , π ) β€ π .The definition extends to a set πΉ of piecewise-affine trajectories,i.e., H π -captures πΉ if H π -captures each π in πΉ . A possible problemto consider is: Given a set of pwa trajectories πΉ and π β R β₯ , constructan adha H such that H π -captures πΉ . The construction of a universal ynthesis of Hybrid Automata with Affine Dynamics from Time-Series Data automaton, describing every possible behavior, trivially satisfiesthe constraint but is not a useful model. Our goal is to constructa model with a reasonable amount of behavior by introducing aminimality criterion that we formally discuss later.Problem 1 (Synthesis).
Given a set of pwa trajectories πΉ and π β R β₯ , construct an adha H such that H π -captures πΉ and satisfiesa minimality criterion. We propose an approach that processes one trajectory π in πΉ at a time and proceeds in two stages. Given a hybrid automaton H and a pwa trajectory π , in the first stage we check whether H π -captures π , which we call a membership query . In the secondstage, if π is not π -captured, we modify H such that it π -captures π . This modification may consist of several changes to the model:increasing the invariants and guards, adding new transitions, andadding new locations. We prioritize the modifications in the ordergiven above to minimize the number of locations.In the next three sections we present algorithmic approaches totransforming time series to pwa trajectories, solving the member-ship query, and performing the model update. In the first phase of our algorithmic framework we construct a pwatrajectory from a time series π . Recall that π is supposed to be thesolution of a piecewise-affine dynamical system, i.e., of a sequenceof contiguous solutions of systems of the form (cid:164) x ( π‘ ) = π΄ π x ( π‘ ) + b π with x ( ) = x . We simplify the problem of finding π by onlyconsidering switching times of π from the domain of π .We thus need to solve the following simpler problem. Given atime series π with domain π· and a value πΏ β R β₯ , find an affinedynamical system (cid:164) x ( π‘ ) = π΄ x ( π‘ ) + b and an initial state x ( ) = x such that the solution π satisfies β₯ π ( π‘ ) β π ( π‘ )β₯ β€ πΏ for every π‘ β π· ,or determine that no such system exists. We pose the problem offinding π as a parameter identification problem where the param-eters are the coefficients of π΄ , b , and x . This can be written as aquery to an optimization tool in combination with an ODE solver(we refer to Section 7 for implementation details). Given concreteparameter values, i.e., instances of π΄ , b , and x , the ODE solvercan compute the solution π corresponding to the affine dynamicalsystem. We can hence evaluate the norm β₯ π ( π‘ ) β π ( π‘ )β₯ at all timepoints π‘ β π· . The optimization tool thus has to find a solution π such that this norm at those time points is less than πΏ .We can use the above algorithm for solving the original problemof finding a pwa trajectory. The main idea is to maximize the du-ration in which we can use the same dynamics. Denote the timepoints of π by π‘ < Β· Β· Β· < π‘ π . We first find the maximum time point π‘ π such that the above-described algorithm finds a solution (e.g., usingbinary search). Then we iteratively solve the same problem for thetime-series suffix from π‘ π to π‘ π , until finally π‘ π = π‘ π . Note that weonly need to identify x for the first piece, as for subsequent piecesthe initial state is determined by x and the previous dynamics. In this section we formalize and solve the membership query. Givenan adha H , a pwa trajectory π , and a value π β R β₯ , the fundamen-tal problem we need to solve is to determine if H π -captures π . Wereduce this problem to checking whether for a given pwa trajectory π and a given path π in H there exists an execution π following π such that π ( π , π ) β€ π . We apply this check to every path π in H oflength equal to the number of pieces in π . We provide a solutionby restricting π and π to switch synchronously, which allows us toevaluate the pieces consecutively. Definition 5.1.
An execution π of an adha H is synchronized with a pwa trajectory π , denoted by π β₯ π , if dom ( π ) = dom ( π ) and t s ( π ) = t s ( π ) .Problem 2 (Membership). Given a path π in an adha H , a pwatrajectory π , and π β R β₯ , determine if there exists a synchronizedexecution π of H with π (cid:123) π and π ( π , π ) β€ π . Our membership algorithm uses reachability analysis to approx-imate the states that the synchronized executions of H can reach. Definition 5.2.
Given an π -dimensional pwa trajectory π and π β R β₯ , an π -tube of π is the function T ( π , π ) : R β₯ β P ( R π ) such that T ( π , π )( π‘ ) = π΅ π ( π ( π‘ )) . Definition 5.3.
Given an adha H , a path π β paths (H) , a pwatrajectory π , and π β R β₯ , the synchronized reachable set , startingfrom a set π β R π and following π , is defined as SReach ( π, π, π , π ) : = { x β R π : β π β exec (H) , π (cid:123) π,π β₯ π , π β π, π ( π‘ ) β T ( π , π )( π‘ ) β π‘ β dom ( π ) and π πππ = x } . For Problem 2, an execution in H satisfying the correspondingconstraints exists if SReach ( π, π, π , π ) is nonempty. Note that theconverse is not true due to unsynchronized executions.Proposition 5.4. Let H be an adha, π be a pwa trajectory, π β R π , and π β R β₯ . If SReach ( π, π, π , π ) is nonempty for some π β paths (H) , then H π -captures π . We inductively construct the synchronized reachable set for apwa trajectory π by computing the synchronized reachable set foreach affine piece of π . Concretely, given an initial set π , a path π = π , . . . , π π in H , and a pwa trajectory π with t s ( π ) = π‘ , . . . , π‘ π ,we define the synchronized reachable sets π : = π, π π : = SReach ( π π β , π π , π β [ π‘ π β ,π‘ π ] , π ) for 1 β€ π β€ π. (1)Observe that π π is equal to SReach ( π, π, π , π ) . We now present a method to approximate the synchronized reach-able set for a pwa trajectory π with just one piece, starting froma polytope π and following a path π of length one in H , that is, SReach ( π, π, π , π ) . This is a special case of Problem 2 where π isan affine trajectory and the path π in H is a single location π . Asobserved before, checking emptiness of the synchronized reach-able set is equivalent to checking whether there exists of an affinetrajectory π in the π -tube of π , starting from the given polytope π ,with the same time domain as π , and following the dynamics of π .Remark 1. Without loss of generality we restrict ourselves to lineardynamics, which are equivalent to affine dynamics under an appro-priate transformation: Add an extra variable π¦ to an affine system (cid:164) x = π΄ x + b as (cid:164) x = π΄ x + b π¦ where π¦ is constant (i.e., (cid:164) π¦ = ). Hence wealso consider hybrid automata with linear dynamics (ldha), whichmeans that the flow function has the signature Flow : Q β R π Γ π . iriam GarcΓa Soto, Thomas A. Henzinger, and Christian Schilling . πππ -tube π -tube (a) An affine trajectory π (black) and the π -tube around π (gray). The light red tube con-sists of all possible executions π followingsome other affine dynamics emerging from T ( π , π ) ( ) . The execution π (yellow) alwaysstays inside the gray tube. (b) Tube partition at apoint in time. Green:under-approximationof states whose execu-tions stay inside. Red:over-approximationof states whose execu-tions eventually leave.Yellow: undecided. Figure 1: Illustration of reachability computations.
Figure 1(a) illustrates that computing the exact synchronizedreachable set is not trivial. Hence we settle for an approximatesolution by successive polytope refinements into three regions,corresponding to the respective executions emerging from thoseregions, as illustrated in Figure 1(b): an under-approximation ofthe states in π whose executions definitely stay inside the tube, anover-approximation of the states in π whose executions definitelyleave the tube, and the remaining states that are undetermined. Insummary, we want to achieve the following goals:(G1) determine whether SReach ( π, π, π , π ) is empty,(G2) (approximately) compute SReach ( π, π, π , π ) , and(G3) refine the polytope π to improve the approximation.We next discuss in detail how to achieve these goals. We now work toward an algorithm for achieving goal (G1). Asargued before, solving the emptiness problem exactly is not trivial.A sufficient condition is to compute an over-approximation andshow emptiness for that set.
SReach ( π, π, π , π ) is empty if and only ifthere exists a time point π‘ β [ ,π ] such that SReach ( π, π, π β [ ,π‘ ] , π ) is empty. We can generalize this observation to sets of points π β² β π .Observe that SReach ( π, π, π , π ) = SReach ( π β² , π, π , π ) βͺ SReach ( π \ π β² , π, π , π ) , so if for each point x in π β² there exists a time point π‘ such that the execution emerging from x leaves the tube, we canremove the set π β² from π . We recall a classic result. Definition 5.5.
The reachable region from π β R π following thelinear dynamics described by π΄ β R π Γ π at time π‘ is defined as Reach ( π, π΄, π‘ ) : = { π π΄π‘ Β· x : x β π } .With π΄ = Flow ( π ) β R π Γ π , we know that Reach ( π, π΄, π‘ ) includesthe points of all executions π at time π‘ such that π (cid:123) π startingfrom x β π . Moreover, π ( π‘ ) belongs to the π -tube around π at time π‘ . Therefore, Reach ( π, π΄, π‘ ) β© T ( π , π )( π‘ ) is an over-approximation of SReach ( π, π, π β [ ,π‘ ] , π ) , providing a sufficient emptiness check. Algorithm 1
Over-approximation of
SReach
Require:
Matrices
π΄, π΅ β R π Γ π , a point x β R π , a polytope π ,values π,π β R β₯ , and a number π β N + . πΏ := π / π {uniform time step for sampling} for π β [ , . . . , π ] do x := ExpMatrix ( π΅ ππΏ ) Β· x π΅ π := Ball ( x , π ) {tube at time point ππΏ } π π := Reach ( π π β , π΄, πΏ ) β© π΅ π if isempty( π π ) then return π π return π π Proposition 5.6.
Emptiness of
Reach ( π, π΄, π‘ ) β© T ( π , π )( π‘ ) im-plies emptiness of SReach ( π, π, π β [ ,π‘ ] , π ) , which implies emptinessof SReach ( π, π, π , π ) . Proposition 5.6 suggests an algorithm for showing emptiness of
Reach ( π, π΄, π‘ ) β© T ( π , π )( π‘ ) at sampled time points π‘ β [ ,π ] . Ob-serve that a finer time sampling provides a more accurate approxi-mation, and a better chance to show emptiness if SReach ( π, π, π , π ) = β . For a uniform sequence of π time points of delay πΏ , Algorithm 1performs the above sufficient check numerically. Recall that π΄ isthe flow of location π , Flow ( π ) , and the linear trajectory π is givenas the tuple ({[ ,π ]} , { π΅ }) . Algorithm 1 takes as input two ma-trices π΄, π΅ β R π Γ π , a point x β R π , a polytope π β R π , twovalues π,π β R β₯ , and a natural number π >
0. For the π -thtime step πΏ = π / π , the algorithm computes π π΅ππΏ Β· x , where π π΅ππΏ is obtained with the function ExpMatrix ( π΅ ππΏ ) . Then Ball ( x , π ) constructs the ball π΅ π ( x ) , and Reach ( π π β , π΄, πΏ ) computes the set Reach ( π π β , π΄, πΏ ) , which is intersected with the ball for construct-ing π π .Proposition 5.7 (Soundness). Algorithm 1 returns an empty setonly if
SReach ( π, π, π , π ) is empty. Proof. Assume that the algorithm returns an empty set but
SReach ( π, π, π , π ) is nonempty. Then there is a point x β π with π π΄π‘ Β· x β T ( π , π )( π‘ ) for every π‘ = ππΏ , 1 β€ π β€ π . Hence π π΄ππΏ Β· x β π π ,which contradicts the condition in line 6. β‘ Proposition 5.8 (Robust completeness).
Let π be a polytopein R π , π΄ β R π Γ π , π a linear trajectory with dom ( π ) = [ ,π ] , and π > such that for every x β π there exists π‘ β [ ,π ] with π ( Reach ({ x } , π΄, π‘ ) , T ( π , π )( π‘ )) > π . Then there exists a finite num-ber π such that Algorithm 1 returns an empty set. Proof. Fix x β π and π‘ β [ ,π ] such that π ( Reach ({ x } , π΄, π‘ ) , T ( π , π )( π‘ )) > π . Then, by continuity of the distance function,there exists a time π‘ ππ₯ππ‘ β [ , π‘ ] such that π ( Reach ({ x } , π΄, π‘ ππ₯ππ‘ ) , T ( π , π )( π‘ ππ₯ππ‘ )) = π‘ β² β [ π‘ ππ₯ππ‘ , π‘ ] , π ( Reach ({ x } , π΄, π‘ β² ) , T ( π , π )( π‘ β² )) >
0. Let us denote π‘ β π‘ ππ₯ππ‘ = πΏ x . Compute the infimumof πΏ x for every x β π , denoted as πΏ β . Then, choose π > π / πΏ β . β‘ Remark 2.
The assumption on Proposition 5.8 about π is necessaryin general because π and the π -tube image are compact. Since π and π β² : = { x β π : Reach ({ x } , π΄, π‘ ) β T ( π , π )( π‘ ) β π‘ } are topologicallyclosed, π β² \ π is not topologically closed. ynthesis of Hybrid Automata with Affine Dynamics from Time-Series Data Algorithm 2
Synchronization check for linear trajectories
Require:
Two matrices
π΄, π΅ β R π Γ π , states x , y β R π , and values π,π β R β₯ . π ( π‘ ) := ExpMatrix ( π΄π‘ ) Β· y π ( π‘ ) := ExpMatrix ( π΅π‘ ) Β· x β ( π‘ ) := π ( π‘ ) β π ( π‘ ) if β₯ x β y β₯ β€ π and β₯ π ( π ) β π ( π )β₯ β€ π then π£ := 0 for β€ π β€ π do π£ πππ₯ := Max ( Abs ( Proj ( β ( π‘ ) , π ))) , [ ,π ]) π£ := max ( π£ πππ₯ , π£ ) if π£ β€ π then return True return
False
Algorithm 1 is a sufficient check: the result is empty only if
SReach ( π, π, π , π ) is empty. Next we consider membership of π inthe π -tube of π where π starts from a fixed point x in π . We can achieve goal (G2) (and hence goal (G1)) for a singleton set π = { x } . In other words, for a fixed starting point π ( ) = x we candecide if π ( π‘ ) β T ( π , π )( π‘ ) for every π‘ β [ ,π ] . We consider the casewhere x β T ( π , π )( ) . We can easily determine if π ( π ) β T ( π , π )( π ) (e.g., by executing Algorithm 1 with π = π ( π ) β T ( π , π )( π ) , the goal is to compute the maximum of π ( π ( π‘ ) , T ( π , π )( π‘ )) over time interval [ ,π ] . Our approach to thatproblem involves solving 2 π optimization problems.Proposition 5.9 (Theorem 4 in [14]). Let x β R π be a point and π΄ β R π Γ π a matrix with rational coefficients. Then Reach ({ x } , π΄, π‘ ) is computable for every time π‘ . We summarize the procedure in Algorithm 2. The inputs are twomatrices
π΄, π΅ β R π Γ π , states x , y β R π , and values π,π β R β₯ .Recall that π΄ is the flow of location π , Flow ( π ) , and the linear trajec-tory π is given as the tuple ({[ ,π ]} , { π΅ }) . Initially, the algorithmdefines the linear trajectories π ( π‘ ) and π ( π‘ ) and their difference β ( π‘ ) . If the norm of this difference is less than π for π‘ = π‘ = π ,the algorithm computes the maximum (Max function) over [ ,π ] ofthe absolute values (Abs function) for each coordinate π of β ( π‘ ) , thatis, Proj ( β ( π‘ ) , π ) . The algorithm returns True if the maximum dis-tance between π and π is less than π , and False otherwise. Thus thealgorithm determines emptiness of
SReach ({ x } , π, π , π ) for x β π .Proposition 5.10. Algorithm 2 returns
False if and only if
SReach ({ x } , π, π , π ) is empty. We assume a numerically sound optimization tool in practice.Algorithm 2 gives us a way to obtain an under-approximation of
SReach ( π, π, π , π ) : apply Algorithm 2 to every vertex of π and con-struct the convex hull of the vertices for which Algorithm 2 returns True . Next we prove that this set is contained in
SReach ( π, π, π , π ) .Proposition 5.11. Let π be a convex polytope, π΄ = Flow ( π ) , π be a linear trajectory with domain [ ,π ] , and π be a value in R β₯ .Then, SReach ( π, π, π , π ) = Reach ( π, π΄,π ) if SReach ({ v } , π, π , π ) isnot empty for every v β vert ( π ) . Proof. The inclusion
SReach ( π, π, π , π ) β
Reach ( π, π΄,π ) is ob-vious. Let SReach ( v , π, π , π ) β β for every v β vert ( π ) . We wantto show that for every point x β π , Reach ({ x } , π΄, π‘ ) belongs to T ( π , π )( π‘ ) for all π‘ β [ ,π ] . Assume there exist x β π and π‘ β [ ,π ] with Reach ({ x } , π΄, π‘ ) β T ( π , π )( π‘ ) , so SReach ({ x } , π, π , π ) = β .We know that Reach ({ x } , π΄, π‘ ) β Reach ( π, π΄, π‘ ) . So Reach ( π, π΄, π‘ ) β T ( π , π )( π‘ ) . Moreover, T ( π , π )( π‘ ) is convex for each π‘ β [ ,π ] .For any polytope π and convex set πΆ it holds that π β πΆ if andonly if vert ( π ) β πΆ . Therefore, Reach ( π, π΄, π‘ ) β T ( π , π )( π‘ ) if andonly if Reach ( vert ( π ) , π΄, π‘ ) β T ( π , π )( π‘ ) , i.e., Reach ({ v } , π΄, π‘ ) β T ( π , π )( π‘ ) for each v β vert ( π ) and π‘ β [ ,π ] . By assumption, SReach ({ v } , π, π , π ) β β for each v β vert ( π ) . Using Proposi-tion 5.10, Reach ({ v } , π΄, π‘ ) β T ( π , π )( π‘ ) for each v β vert ( π ) .Hence Reach ( π, π΄, π‘ ) β T ( π , π )( π‘ ) for each x β π : a contradic-tion. β‘ Corollary 5.12.
If Algorithm 2 returns
True for all vertices v β vert ( π ) , then SReach ( π, π, π , π ) = Reach ( π, π΄,π ) . Recall that π is a polytope, π is a linear trajectory π ( π‘ ) = π π΅π‘ Β· x with time domain [ ,π ] , π is a location in some ldha H with Flow ( π ) = π΄ , π is a value in R β₯ , and π > π of the synchronized reachable set. If π is nonempty,we can use Algorithm 2 from Section 5.3 for every vertex of π , andif the algorithm returns True for some vertex, we have a nonemptyunder-approximation and can conclude membership of π in H . IfAlgorithm 2 returns False for all vertices, we cannot conclude.Next we propose a new procedure, which together with Proposi-tion 5.11 suggests an algorithm for computing a more precise under-approximation of
SReach ( π, π, π , π ) . Finally, these procedures to-gether induce an algorithm to refine the over- and under-approximations.Intuitively, recalling Figure 1(b), this refinement narrows the dis-crepancy between the the over-approximation (yellow) and theunder-approximation (green).First we observe that the over-approximation π is a convex poly-tope. The idea is to contract this polytope to a new polytope. Givena value πΏ β R β₯ , we define the πΏ -contraction of π as follows. Definition 5.13.
Let π be a polytope, πΏ be a value in R β₯ , and constr ( π ) = { a x βΌ π , . . . , a π x βΌ π π } . The πΏ -contraction of π is the polytope π πΏ : = { x β R π : a x βΌ π , . . . , a π x βΌ π π } where π π = π π if βΌ is β = β and π π = π π β πΏ β₯ π π β₯ if βΌ is β β€ β, for every 1 β€ π β€ π .We can hence take the over-approximation π , compute the πΏ -contraction π β² , and then apply Algorithm 2 to all vertices of π and π β² . Ultimately we may have to repeat this contraction several times(at most β π / πΏ β times, where π is the diameter of π ). In the end, sincewe know that the true synchronized reachable set is convex, wecan take the convex hull of all those vertices for which Algorithm 2returned True (i.e., these vertices belong to the synchronized reach-able set). We summarize the refinement in Algorithm 3, where theprocedure Contract applies a πΏ -contraction.In principle, now that we have two polytopes π and π β² over-approximating and under-approximating the synchronized reach-able set, respectively, a natural additional refinement procedure iriam GarcΓa Soto, Thomas A. Henzinger, and Christian Schilling Algorithm 3
Polytope refinement
Require:
A polytope π , two matrices π΄, π΅ β R π Γ π , a state x β R π ,and three values π,π, πΏ β R β₯ . π + := β while True do π β := β for v β vert ( π ) do if Algorithm 2(
π΄, π΅, x , v , π,π ) then π + := π + βͺ { v } else π β := π β βͺ { v } if π β = β then return chull ( π + ) π := Contract ( π, πΏ ) Algorithm 4
Membership query for a single piece
Require:
A polytope π , two matrices π΄, π΅ β R π Γ π , a state x β R π , three values π,π, πΏ β R β₯ , and a value π β N + . π := Algorithm 1( π΄, π΅, x , π , π,π, π ) if isempty( π ) then return β , β x := ExpMatrix ( π΄π ) Β· x π := Algorithm 3( π , β π΄ , β π΅ , x , π , π , πΏ ) return π , π can be conceived where one iteratively tries to enlarge the under-approximation or shrink the over-approximation. We did not in-vestigate this direction because the above scheme is already veryprecise in practice. (In fact, we rather observed that the approxima-tions become too precise; see the further discussion in Section 7.1.) Algorithm 4 summarizes the overall procedure for computing bothan under-approximation and an over-approximation of the synchro-nized reachable set for a linear trajectory π . We first use Algorithm 1to compute the over-approximation π . If the over-approximation isempty, we can conclude that π is not π -captured. Otherwise, takingthe end state x of π and inverting the dynamics ( (cid:164) π inv ( x ) = β (cid:164) π ( x ) ),we use Algorithm 3 to compute the under-approximation π .We illustrate the algorithm and the generalization to multiplepieces with the following parametric linear trajectories: (cid:164) x = (cid:18) β (cid:19) x , x ( ) = (cid:18) (cid:19) (2) (cid:164) y = (cid:18) β πΌ β (cid:19) y , y ( ) = y (3)System (2) is fixed and takes the role of the linear trajectory π while system (3) models the location of an ldha. In the following,we fix the parameter value πΌ and ask whether there exists an initialstate y such that the corresponding execution of system (3) issynchronized with π . In Figure 2 we plot the executions for πΌ = . y = x ( ) . It can be seen that for a timehorizon of 4 π we need to choose π larger than βΌ . t x(t)y(t)difference (a) Two executions of system (2) (blue) and system (3) (red) andtheir difference (green). t dim 1dim 2 (b) Difference of the trajectoriesprojected to time and one di-mension. Figure 2: Trajectories and difference of systems (2) and (3) with parameters πΌ = . and y = x ( ) . In Figure 3 we plot the results for π = . πΌ = πΌ = .
01, respectively. In Figure 3(a) we see the under-approximationcomputed by the algorithm in light green. The dark green set is asimplified under-approximation that we use to handle complexity,further described in Section 7.1. Similarly, the dark yellow set isthe over-approximation computed by the algorithm, while the lightyellow set is a simplified over-approximation. It can be seen thatthe gap between the over-approximation (dark yellow) and theunder-approximation (light green) is very narrow, indicating thatthe refinement procedure (Algorithm 3) is precise. Also note thatthe true synchronized reachable set in this case is a Euclidean ballbecause, while the executions all follow the same dynamics as π ,those executions starting from a state outside this ball rotate around π and eventually leave the tube (since the tube does not rotate).In Figure 3(c) we plot the intermediate sets for the same execu-tions but modeled as pwa trajectories with 23 pieces, starting withthe set at time 0. In theory, the settings with a single piece and 23pieces of the same dynamics are equivalent; however, due to thesimplifications of the approximations for each piece, the approxi-mations lose precision in the latter case. Still, the approximationsin the last piece are sufficiently precise to prove that the under-approximation (green set) is nonempty and hence we can concludewith a positive answer to the membership query. In the last sub-plot we depict a random sampling from the over-approximationwhere we apply Algorithm 2 to check whether the state indeedcorresponds to a synchronized execution (green dot) or not (reddot). Figure 3(b) shows the setting for πΌ = .
01 with similar results.
In this section we describe a procedure to solve Problem 1 andpropose a minimality criterion. The procedure tackles the problemby evaluating the given pwa trajectories in an online fashion.
For a given adha H , a pwa trajectory π , and a value π β R β₯ , ourprocedure searches for a path π in H such that the membershipquery of π in H is positive. If there is no such path in H , theprocedure modifies H such that the modified adha includes sucha path. The path selection and the corresponding modifications ofthe adha are chosen in the following order: 1) increasing invariants ynthesis of Hybrid Automata with Affine Dynamics from Time-Series Data (a) Analysis for πΌ = and a singlepieces of duration π . (b) Analysis for πΌ = . and asingle pieces of duration π .
11 1.370.33 1.35-0.43 0.93-1.07 0.24-1.39-0.52-1.32 -1.13-0.85 -1.41-0.14 -1.280.61 -0.781.18-0.051.41 0.691.23 1.230.69 1.41-0.05 1.18-0.780.61-1.28 -0.14-1.41 -0.85-1.13 -1.32-0.52 -1.390.24-1.070.93 -0.431.35 0.331.37 11 1111 1.370.33 1.35-0.43 0.93-1.07 0.24-1.39-0.52-1.32 -1.13-0.85 -1.41-0.14 -1.280.61 -0.781.18-0.051.41 0.691.23 1.230.69 1.41-0.05 1.18-0.780.61-1.28 -0.14-1.41 -0.85-1.13 -1.32-0.52 -1.390.24-1.070.93 -0.431.35 0.331.37 11 11 (c) Analysis for πΌ = and pieces of uniform duration (only thefirst and last four intermediate results are shown). The last subplotshows , samplings from the final dark yellow region. Figure 3: Analysis results for systems (2) and (3) with π = . and πΌ = . or πΌ = , respectively. and guards, 2) adding new transitions, and 3) adding new locations.The rationale is to keep the number of locations as small as possible.Formally, we define a tuple mod = ( π π , π π‘ , π π ) that keeps trackof the above modifications, where π π tracks the number of newlocations, π π‘ tracks the number of new transitions, and π π tracksthe number of modified constraints (invariants and guards). Wewill use the tuple mod for path selection in a lexicographic orderwhere we aim to find the minimal tuple. For instance, the tuple ( , , ) , representing no modifications at all, is selected over anyother tuple; the tuple ( , , ) , representing two transition additions,is selected over the tuple ( , , ) , representing a location addition. For a given adha H , a pwa trajectory π , and a value π β R β₯ , weupdate H for each affine piece in π if required. We describe how H is modified for a concrete piece of π according to a location π thatmay either be part of H or be a new location to be added to H . Definition 6.1.
Consider an adha H , a path π in H with len ( π ) = π β
1, an existing or new location π , a pwa trajectory π , representedby the tuple ({[ , π‘ ] , . . . , [ π‘ π β , π‘ π ]} , { π΄ , . . . , π΄ π } , { b , . . . , b π }) , apolyhedron π , and a value π β R β₯ . A π -update of H with respectto π , π , and π is an adha H β² , denoted by Mod
π,ππ,π (H) , such that Q β² = Q βͺ { π } , E β² = E βͺ {( last ( π ) , π )} , and the remaining compo-nents depend on whether π exists in π or is a new location.If π β π : 1) Flow β² β‘ Flow , 2)
Inv β² β Q \{ π } β‘ Inv and
Inv β² ( π ) = chull ( Inv ( π ) βͺ R πΌ ) , and 3) Grd β² β E \{( last ( π ) ,π ) } β‘ Grd and
Grd β² ( last ( π ) , π ) = chull ( Grd ( last ( π ) , π ) βͺ R πΊ ) .If π β Q : 1) Flow β² β Q β‘ Flow and
Flow β² ( π ) = ( π΄ π , b π ) , 2) Inv β² β Q β‘ Inv and
Inv β² ( π ) = chull (R πΌ ) , and 3) Grd β² β E β‘ Grd and
Grd β² ( last ( π ) , π ) = R πΊ . Here R πΌ = (cid:208) π‘ β[ π‘ π β ,π‘ π ] SReach ( π π β , π Β· π,π β [ π‘ π β ,π‘ ] , π ) , with π π β as defined in (1), and R πΊ = T ( π , π )( π‘ π ) . Observe that exec (H) β exec ( Mod
π,ππ,π (H)) . Next we define atree capturing the adha updates for every affine piece in π . Definition 6.2.
Given an π -dimensional adha H and a pwatrajectory π with π pieces, an exploration tree for H and π is T = ( N , E ) with π layers (not counting the root node as a layer). Eachnode π β N is represented as a tuple ( π, H , mod , s ) where π is apath in an adha H , mod is a triple of integers ( π π , π π‘ , π π ) , and s isa four-valued variable called status (with meanings 0 : βunexploredβ,1 : βactivatedβ, 2 : βexploredβ, 3 : βdeactivatedβ).Observe that exploration trees for H and π can only differ inthe status. The set of all exploration trees for H and π is denotedby S (H , π ) , and we call all trees belonging to S (H , π ) similar . Wemay add a subscript to the elements in the node π = ( π, H , mod , s ) (i.e., write π π etc.) for clarity. We define, for an initial polyhedron π and a value π β R β₯ , an exploration tree T π,π β S (H , π ) suchthat the root node is ([ ] , H , ( , , ) , ) , where [ ] is the empty path.Each node ( π, H , ( π π , π π‘ , π π ) , ) in layer π β
1, for 0 < π β€ π , where Q H = { π , . . . , π π } , has π + child nodes . The first π nodes are: ( π Β· π , Mod
π,ππ ,π (H) , ( π π , π π‘ + π , π π + π ) , ) , . . . , ( π Β· π π , Mod
π,ππ π ,π (H) , ( π π , π π‘ + π π , π π + π π ) , ) , where π π = ( last ( π ) , π π ) β E H and π π = π π is the number of constraint modifications for invariants and guardswith respect to H , for every 1 β€ π β€ π . The last child node is: ( π Β· π, Mod
π,ππ,π (H) , ( π π + , π π‘ + , π π ) , ) , where π is a new location with Flow ( π ) = ( π΄ π , b π ) .The paths from root to leaves in an exploration tree representall possible paths in updated adhas, given the initial adha H ,for exploring membership of π . An upper bound on the numberof paths is ( π + π ) π , where π = | π | is the number of locationsin H and π is the number of pieces in π . The complexity for themembership query is in O( π ( π )) for some polynomial π in thedimension π . Hence the complexity for a membership check in eachpath of the exploration tree is upper-bounded by O(( π + π ) π ππ ( π )) .We introduce a strategy for partial exploration that minimizesautomaton modifications (according to mod ). Given H and π , a de-cision strategy is a function π· : S (H , π ) Γ N β N that determinesthe next node to be analyzed in an exploration tree. A decision strat-egy is combined with a tree update in order to activate and explorenodes or discard useless nodes. We say that a node is unexplored when its status is 0. We can explore a node when it is activated (status 1). After exploration, if the membership query is positive,we set the status to 2 ( explored ) and otherwise to 3 ( deactivated ).Child nodes of deactivated nodes need not be explored further. Definition 6.3.
Given an adha H , a pwa trajectory π , a poly-hedron π , and a value π β R β₯ , an π -tree update function, upd π : S (H , π ) Γ N β S (H , π ) , maps a tree T and a node π in the π -thlayer to a similar tree such that s π = SReach ( π, π π , π β [ ,π‘ π ] , π ) = β , and s π = s π β² = π β² β children ( π ) otherwise, and leaves the status of all other nodes unchanged.Given a set of nodes π , we denote by Act ( π ) the set of nodeswith activated status, i.e., { π β π : s π = } . Our decision strategy iriam GarcΓa Soto, Thomas A. Henzinger, and Christian Schilling π (cid:164) π₯ = π₯ β€ π₯ β€ π (cid:164) π₯ = β π₯ β€ π₯ β€ β€ π₯ β€ (a) Input adha H . . . . . . π‘ π₯ (b) pwa trajectory π . ([ ] , H , ( , , ) , )( π , H , ( , , ) , ) ( π , H , ( , , ) , )( π π , H , ( , , ) , )( π π , H , ( , , ) , )( π π , H , ( , , ) , )( π , H , ( , , ) , ) (c) Partial view of a search tree for π = . . Red nodes are deactivated,blue nodes have been explored, and black nodes are activated (andwe omit their child nodes). The algorithm selects the path to thegreen node for a membership query. Figure 4: An adha H , a pwa trajectory π with two pieces,and a (partial) exploration tree for H and π . minimizing mod is π· (T , π ) = arg min π β² β Act ( N ) mod π β² , assuming thatarg min returns one node if several nodes minimize the mod value. Example 6.4.
Figure 4 shows an example of an intermediate stateof an exploration tree for a given adha and a pwa trajectory withtwo pieces. The root node has been described before. For the remain-ing tree nodes we represent the automata H π only symbolically.The first piece of π follows the dynamics (cid:164) π₯ = β π₯ for 0 . π , π , ora new mode π ; hence the root node expands to three new nodes.The node with path π requires a modification of the invariant of π because dwelling in that mode for 0 . π -tube (negative membership query), this node status is setto 3 (deactivated) and none of the child nodes are explored further.The node with the new location π has a βlocation entryβ in themodification tuple. The node with path π does not require anymodifications (i.e., H = H ) and is hence chosen as the next nodefor exploration. Now we consider the second piece with dynamics (cid:164) π₯ = π₯ for 0 . π . The exploration workslike before, only that this time we need to add a transition from π to the next location in all three cases (since π does not haveany outgoing transitions in H = H yet). The 2-leaf with the path π π has the highest priority and we perform a membership queryfor it. In this case, the query returns a positive answer and thealgorithm outputs the automaton H , which looks like H but withan additional transition and an extended invariant in location π .Algorithm 5 shows the procedure for a model update given aninitial adha H , a pwa trajectory π , and a value π β R β₯ . Thefunction InitTree constructs the exploration tree T π,π for thepolyhedron π = T ( π , π )( ) . Then the algorithm starts exploringfrom the root node (line 2) and subsequently explores nodes driven Algorithm 5
Hybrid model update
Require:
An adha H , a pwa trajectory π and π β R β₯ . Ensure:
An adha H β² such that it π -captures π . T := InitTree (H , π , π ) π := RootNode (T ) T := upd π (T , π ) while π not in bottom layer with s π = do π := π· (T , π ) T := upd π (T , π ) H β² := H π return H β² Algorithm 6
Synthesis of adha from pwa trajectories
Require:
A finite set of pwa trajectories πΉ and π β R β₯ . H := β {empty automaton with no location} for π in πΉ do H := Algorithm 5( H , π , π ) return H by the decision strategy (line 5), which chooses activated nodeswith minimum mod component and iteratively activates every childnodes and deactivates the current node or sets it to explored (line 6).The algorithm returns the updated model H β² . Finally, Problem 1is solved by iteratively running Algorithm 5 over every trajectory π β πΉ and modifying the adha, as shown in Algorithm 6.Proposition 6.5. Given an adha H , a pwa trajectory π , and avalue π β R β₯ , Algorithm 5 provides an updated adha π -capturing π and minimizing the number of modifications. Proof. Given an adha H , a pwa trajectory π , and a value π β R β₯ , Algorithm 5 proceeds as follows. First, the algorithmconstructs an initial exploration tree T π,π (line 1) whose nodescontain all the possible modifications of H with unexplored status(0). Then, the algorithm sets π as the root node (line 2) and appliesthe π -tree update over the initial exploration tree and π (line 3). Thisupdate sets the status for all nodes in the first layer to 1 and forthe root node to 2 because SReach ( π, β , π β [ , ] , π ) = π . Next, thealgorithm iterates (line 4) as follows. The decision strategy π· (T , π ) selects the node of the search tree with minimum mod value andactivated status (1). Then, the π -tree update requires to check if theadha H π in the π -th layer of the exploration tree π -captures the π first pieces of π . If these pieces are not captured, the status ofthe node is deactivated (set to 3) and all child nodes will remainunexplored (with status 0) forever. If H π π -captures π β [ ,π‘ π ] , thestatus of the node is set to 2 (explored) and the status of all childnodes is activated (1). The loop runs until π is a node at the bottomlayer with status 2, which means that H π π -captures π and that mod is minimum due to the decision strategy. The algorithm terminatesbecause, in the worst case, it will choose the path in the search treewhere a new mode is added for each piece of π ; clearly the adha atthe leaf of that path π -captures π . β‘ Theorem 6.6.
Given an adha H , a set πΉ of pwa trajectories, anda value π β R β₯ , Algorithm 6 solves Problem 1. ynthesis of Hybrid Automata with Affine Dynamics from Time-Series Data In this section we describe our implementation and evaluate it: inthe first two examples we obtain the pwa trajectories from randomexecutions with perturbed dynamics from a given adha model; ina third example we construct the pwa trajectories from time series.
We implemented our approach in
HySynth [1] where we wrotethe high-level synthesis algorithm in Python and the low-levelalgorithms in Julia. For the ODE optimization (both in Section 4and Algorithm 2) we use the libraries
Optim.jl [21] (which usesBrentβs method [8] to find a root in a bracketing interval and guar-antees convergence for functions computable within the interval)and
DifferentialEquations.jl [30] as follows (assuming lineardynamics without loss of generality). Given two π -dimensionalexecutions (cid:164) x = π΄ x , x ( ) = x and (cid:164) y = π΅ y , y ( ) = y , we constructa 3 π -dimensional execution (cid:164) z = πΆ z , z ( ) = z where πΆ = (cid:169)(cid:173)(cid:171) π΄ π΅ π΄ β π΅ (cid:170)(cid:174)(cid:172) z = (cid:169)(cid:173)(cid:171) x y x β y (cid:170)(cid:174)(cid:172) . We are interested in the projection of z ( π‘ ) onto the last π dimensions.Calling this projection w ( π‘ ) , the norm of w ( π‘ ) describes the distancebetween x ( π‘ ) and y ( π‘ ) , i.e., β₯ w ( π‘ )β₯ = π ( x ( π‘ ) , y ( π‘ )) . We query thesolver for each dimension of w ( π‘ ) to find the maximum distance.We use JuliaReach [7] for the set computations and reachabilityanalysis. As mentioned in Section 5.4, the polytopes over- andunder-approximating the synchronized reachable sets constructedduring the membership query grow in complexity, especially forinput trajectories with many pieces. We simplify the sets aftereach piece, i.e., we under-approximate an under-approximation (forwhich
JuliaReach computes a polytope from support vectors intemplate directions) and over-approximate an over-approximation(for which we implemented an algorithm from [13] to computea zonotope in template directions) with octagonal directions (i.e.,axis-parallel or diagonal constraints in two dimensions).
We consider an adha that models a heater with two locations β ON βand β OFF β, as depicted in Figure 5(a) with parameter value π = . . . π and an initial (continuous) state x from Inv ( π ) . Thenwe repeat the following loop. Given a location π and a state x , wefirst compute a matrix π΄ by perturbing the dynamics matrix Flow ( π ) (technically, we only perturb non-zero entries). Then we computethe discrete-time successor of x with the fixed time step π‘ (via x β² : = Reach ({ x } , π΄, π‘ ) ), and we check which of the outgoing transitions of π are enabled for this new state. We continue computing successorstates and collecting enabled transitions until either the state leaves Inv ( π ) or we exceed the maximum dwell time. Then we choosea random transition together with a random time point of those ON (cid:164) π₯ = β π ( π₯ β ) π₯ β€ OFF (cid:164) π₯ = β ππ₯π₯ β₯ π₯ β₯ π₯ β€ (a) Original adha. ON (cid:164) π₯ = β π ( π₯ β ) β€ π₯ β€ OFF (cid:164) π₯ = β ππ₯ β€ π₯ β€ . β€ π₯ β€ . β€ π₯ β€ . (b) adha synthesized from ten simulated executions ( π = . ). ON (cid:164) π₯ = β π ( π₯ β ) . β€ π₯ β€ . OFF (cid:164) π₯ = β ππ₯ . β€ π₯ β€ . . β€ π₯ β€ . . β€ π₯ β€ . (c) adha synthesized from simulated executions ( π = . ). Figure 5: adha models of the heater system. Numbers arerounded to one decimal place. (a) Heater model: Variable x (or-dinate) over time (abscissa). (b) Gearbox model: Phase por-trait. Figure 6: Random simulations: Three blue simulations areobtained from the original model and three green simula-tions are obtained from the synthesized model. that were enabled. The above loop terminates if either there is notransition enabled or we exceed the desired path length.We applied the above procedure to obtain 100 random executionsfrom the heater model. Then we first learned a model from the firstten executions and then continued modifying the resulting adhawith the remaining 90 executions, where we used a precision value π = .
1. (Note that our algorithmic framework behaves exactly thesame way as if we had learned an adha from the 100 executionsat once. The split into two stages is only for illustrative purposes.)We show the intermediate and the final result obtained with ourimplementation in Figure 5(b) and Figure 5(c) respectively, andrandom simulations in Figure 6(a).The first observation is that the discrete structure of the result-ing adha matches exactly the structure of the original model. Thereason why the dynamics matrices of the locations is the same asin the original model is because we did not perturb the dynamics ofthe very first execution in order to obtain a legible flow represen-tation. Still, even though the algorithm is confronted with slightly iriam GarcΓa Soto, Thomas A. Henzinger, and Christian Schilling (cid:164) x = π΄ x π£ β₯ (cid:164) x = π΄ x β€ π£ β€ (cid:164) x = π΄ x β€ π£ β€ (cid:164) x = π΄ x π£ β€ π£ = π£ = π£ = (a) Original adha. (cid:164) x = π΄ x β€ π£ β€ β€ π€ β€ (cid:164) x = π΄ x β€ π£ β€ β€ π€ β€ (cid:164) x = π΄ x β€ π£ β€ β€ π€ β€ (cid:164) x = π΄ x β β€ π£ β€ β€ π€ β€ β€ π£ β€ β€ π€ β€ β€ π£ β€ β€ π€ β€ β€ π£ β€ β€ π€ β€ (b) adha synthesized from ten simulated executions ( π = . ). Figure 7: adha models of the gearbox system. Numbers arerounded to integers; constraints are approximated by boxes. (a) πΏ = . . (b) πΏ = . . (c) Simulations. Figure 8: pwa trajectories (blue) for three time series (red)and different values of πΏ . The last plot shows three simula-tions from the synthesized adha (green). different dynamics in all other executions, it does not add furtherlocations to the adha, thanks to the precision value π . As can beseen, the invariants and guards in the final adha over-approximatethe original guards by π , which is expected by construction.We also applied the algorithm to a two-dimensional gearbox model with variables π£ and π€ from [29] and we refer to that refer-ence for further details about the model. We present the results for10 simulations, a maximum perturbation of 0 . β€ π£ ( ) β€ , π€ ( ) = π / πΏ run time | π | π = . , π = .
07 51 s 3 606 505 , π = .
04 63 s 5 755 6 , , π = .
01 162 s 13 2 ,
798 731 , , π = . , π = .
07 12 s 5 46 13 , π = .
04 12 s 6 51 17 , π = .
01 17 s 10 109 80 , πΏ = .
05 157 s 8 185 7 , , πΏ = .
02 3 ,
115 s 7 101 ,
145 721 , , , Table 1: Benchmark results. The second column shows theallowed error π between pwa trajectories and adha resp. theallowed error πΏ between time series and pwa trajectories.The last two columns show the number of explored treenodes resp. the total number of possible tree nodes. In another experiment we investigate the conversion of timeseries to pwa trajectories. We consider three
ECG signals fromthe PhysioBank database [12]. For the distance value πΏ = .
05 weobtained three pwa trajectories of length 7 in 158 seconds. For πΏ = .
02 we obtained pwa trajectories of respective lengths 10, 11,and 13 in 220 seconds. Using π = . π and πΏ ). As expected, decreasing π results in bigger adha since existing modes can be shared fordifferent pwa trajectory pieces less often. In the ECG benchmarkwe observe that decreasing πΏ can result in smaller adha since theconstructed pwa trajectories are less diverse, even though theyhave more pieces (up to 13 pieces ( πΏ = .
02) compared to 7 pieces( πΏ = . We have presented an automatic synthesis algorithm for comput-ing a hybrid automaton with affine differential dynamics H froma set of time series π respectively from a set of piecewise-affinetrajectories πΉ . Given precision parameters πΏ and π , the main featureof our algorithm is that every time series π in π is πΏ -captured bysome trajectory π in πΉ and that H is guaranteed to π -capture everyfunction π in πΉ , that is, H contains an execution that has distanceat most π from π . Another feature of our algorithm is that it worksonline, meaning that the functions π are processed sequentiallyand we only modify the intermediate automaton models.For future work, hardness of the membership problem for theclass of automata that we considered is open. We currently do notknow if that problem is decidable, and if so, what complexity isrequired to solve it exactly. Another interesting but challenging ynthesis of Hybrid Automata with Affine Dynamics from Time-Series Data extension of our work is to allow for transition switches not at asingle time point but in a whole time interval. ACKNOWLEDGMENTS
This research was supported in part by the Austrian Science Fund(FWF) under grant Z211-N23 (Wittgenstein Award) and the Euro-pean Unionβs Horizon 2020 research and innovation programmeunder the Marie SkΕodowska-Curie grant agreement No. 754411.
REFERENCES [1] 2021. HySynth. https://github.com/HySynth/HySynth.[2] Rajeev Alur and Nimit Singhania. 2014. Precise Piecewise Affine Models fromInput-Output Data. In
EMSOFT (New Delhi, India). Association for ComputingMachinery, Article 3.[3] Dana Angluin. 1987. Learning Regular Sets from Queries and Counterexamples.
Inf. Comput.
75, 2 (1987).[4] Laurent Bako and RenΓ© Vidal. 2008. Algebraic Identification of MIMO SARXModels. In
HSCC , Vol. 4981. Springer.[5] Ezio Bartocci, Jyotirmoy Deshmukh, Felix Gigler, Cristinel Mateis, Dejan Nickovic,and Xin Qin. 2020. Mining Shape Expressions From Positive Examples.
IEEETrans. Comput. Aided Des. Integr. Circuits Syst.
39, 11 (2020), 3809β3820.[6] Alberto Bemporad, Andrea Garulli, Simone Paoletti, and Antonio Vicino. 2005.A bounded-error approach to piecewise affine system identification.
IEEE Trans.Automat. Contr.
50, 10 (2005).[7] Sergiy Bogomolov, Marcelo Forets, Goran Frehse, Kostiantyn Potomkin, andChristian Schilling. 2019. JuliaReach: a toolbox for set-based reachability. In
HSCC . ACM.[8] Richard P. Brent. 1971. An Algorithm with Guaranteed Convergence for Findinga Zero of a Function.
Comput. J.
14, 4 (1971).[9] Giancarlo Ferrari-Trecate and Marco Muselli. 2003. Single-Linkage Clustering forOptimal Classification in Piecewise Affine Regression.
IFAC Proceedings Volumes
36, 6 (2003).[10] Giancarlo Ferrari-Trecate, Marco Muselli, Diego Liberati, and Manfred Morari.2001. A Clustering Technique for the Identification of Piecewise Affine systems.In
HSCC . Springer.[11] Andrea Garulli, Simone Paoletti, and Antonio Vicino. 2012. A survey on switchedand piecewise affine system identification.
IFAC Proceedings Volumes
45, 16(2012).[12] Ary L. Goldberger, Luis A. N. Amaral, Leon Glass, Jeffrey M. Hausdorff, Pla-men Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody, Chung-KangPeng, and H. Eugene Stanley. 2000. PhysioBank, PhysioToolkit, and PhysioNet.
Circulation
SODA . ACM/SIAM.[14] Emmanuel Hainry. 2008. Reachability in Linear Dynamical Systems. In
Logic andTheory of Algorithms . Springer Berlin Heidelberg.[15] Yasmin Hashambhoy and RenΓ© Vidal. 2005. Recursive identification of switchedARX models with unknown number of models and unknown orders. In
CDC .[16] Thomas A. Henzinger. 2000.
The Theory of Hybrid Automata . Springer.[17] A. L. Juloski, S. Weiland, and W. P. M. H. Heemels. 2005. A Bayesian approach toidentification of hybrid systems.
IEEE Trans. Automat. Control
50, 10 (2005).[18] Kun Huang, A. Wagner, and Yi Ma. 2004. Identification of hybrid linear time-invariant systems via subspace embedding and segmentation (SES). In
CDC ,Vol. 3.[19] Fabien Lauer, GΓ©rard Bloch, and RenΓ© Vidal. 2011. A continuous optimizationframework for hybrid system identification.
Automatica
47, 3 (2011).[20] Ramy Medhat, S. Ramesh, Borzoo Bonakdarpour, and Sebastian Fischmeister.2015. A framework for mining hybrid automata from input/output traces. In
EMSOFT . IEEE.[21] Patrick Kofod Mogensen and AsbjΓΈrn Nilsen Riseth. 2018. Optim: A mathematicaloptimization package for Julia.
Journal of Open Source Software
3, 24 (2018).[22] Eberhard MΓΌnz and Volker Krebs. 2005. Continuous Optimization Approachesto the Identification of Piecewise Affine Systems.
IFAC Proceedings Volumes
38, 1(2005).[23] Hayato Nakada, Kiyotsugu Takaba, and Tohru Katayama. 2005. Identification ofpiecewise affine systems based on statistical clustering technique.
Automatica
41, 5 (2005).[24] Sohail Nazari, Bahador Rashidi, Qing Zhao, and Biao Huang. 2016. An IterativeAlgebraic Geometric Approach for Identification of Switched ARX Models withNoise.
Asian J. Control
18, 5 (2016).[25] Necmiye Ozay. 2016. An exact and efficient algorithm for segmentation of ARXmodels. In
ACC . IEEE.[26] N. Ozay, C. Lagoa, and M. Sznaier. 2009. Robust identification of switched affinesystems via moments-based convex optimization. In
CDC . [27] Necmiye Ozay, Constantino M. Lagoa, and Mario Sznaier. 2015. Set membershipidentification of switched linear systems with known number of subsystems.
Automatica
51 (2015).[28] Simone Paoletti, Aleksandar Lj. Juloski, Giancarlo Ferrari-Trecate, and RenΓ© Vidal.2007. Identification of Hybrid Systems: A Tutorial.
Eur. J. Control
13, 2-3 (2007).[29] Pavithra Prabhakar and Miriam GarcΓa Soto. 2016. An algorithmic approach toglobal asymptotic stability verification of hybrid systems. In
EMSOFT . ACM.[30] Christopher Rackauckas and Qing Nie. 2017. Differentialequations.jl - a per-formant and feature-rich ecosystem for solving differential equations in Julia.
Journal of Open Research Software
5, 1 (2017).[31] Jacob Roll, Alberto Bemporad, and Lennart Ljung. 2004. Identification of Piece-wise Affine Systems via Mixed-Integer Programming.
Automatica
40, 1 (2004).[32] Anders Skeppstedt and Mille Jung, Lennart L. Anders Millnert. 1992. Constructionof composite models from observed data.
Int. J. Control
55, 1 (1992).[33] Miriam GarcΓa Soto, Thomas A. Henzinger, Christian Schilling, and LukaZeleznik. 2019. Membership-Based Synthesis of Linear Hybrid Automata. In
CAV ,Vol. 11561. Springer.[34] V. Verdult and M. Verhaegen. 2004. Subspace identification of piecewise linearsystems. In
CDC , Vol. 4.[35] RenΓ© Vidal and Brian D. O. Anderson. 2004. Recursive identification of switchedARX hybrid models: exponential convergence and persistence of excitation. In
CDC , Vol. 1.[36] Jan C. Willems. 1986. From time series to linear system - Part II. Exact modelling.