Linear Hybrid System Falsification With Descent
Linear Hybrid System Falsification Through Descent††thanks: This work was partially supported by a grant from the NSF Industry/University Cooperative Research Center (I/UCRC) on Embedded Systems at Arizona State University and NSF award CNS-1017074.
In this paper, we address the problem of local search for the falsification of hybrid automata with affine dynamics. Namely, if we are given a sequence of locations and a maximum simulation time, we return the trajectory that comes the closest to the unsafe set. In order to solve this problem, we formulate it as a differentiable optimization problem which we solve using Sequential Quadratic Programming. The purpose of developing such a local search method is to combine it with high level stochastic optimization algorithms in order to falsify hybrid systems with complex discrete dynamics and high dimensional continuous spaces. Experimental results indicate that indeed the local search procedure improves upon the results of pure stochastic optimization algorithms.
Keywords:Model Validation and Analysis; Robustness; Simulation; Hybrid systems
Despite the recent advances in the computation of reachable sets in medium to large-sized linear systems (about 500 continuous variables) [1, 2], the verification of hybrid systems through the computation of the reachable state space remains a challenging problem [3, 4]. To overcome this difficult problem, many researchers have looked into testing methodologies as an alternative. Testing methodologies can be coarsely divided into two categories: robust testing [5, 6, 7] and systematic/randomized testing [8, 9, 10, 11].
Along the lines of randomized testing, we investigated the application of Monte Carlo techniques  and metaheuristics to the temporal logic falsification problem of hybrid systems. In detail, utilizing the robustness of temporal logic specifications  as a cost function, we managed to convert a decision problem, i.e., does there exist a trajectory that falsifies the system, into an optimization problem, i.e., what is the trajectory with the minimum robustness value? The resulting optimization problem is highly nonlinear and, in general, without any obvious structure. When faced with such difficult optimization problems, one way to provide an answer is to utilize some stochastic optimization algorithm like Simulated Annealing.
In our previous work , we treated the model of the hybrid system as a black box since a global property, such as convexity of the cost function, cannot be obtained, in general. One question that is immediately raised is whether we can use “local” information from the model of the system in order to provide some guidance to the stochastic optimization algorithm.
In this paper, we set the theoretical framework to provide local descent information to the stochastic optimization algorithm. Here, by local we mean the convergence to a local optimal point. In detail, we consider the falsification problem of affine dynamical systems and hybrid automata with affine dynamics where the uncertainty is in the initial conditions. In this case, the falsification problem reduces to an optimization problem where we are trying to find the trajectory that comes the closest to the unsafe set (in general, such a trajectory is not unique). A stochastic optimization algorithm for the falsification problem picks a point in the set of initial conditions, simulates the system for a bounded duration, computes the distance to the unsafe set and, then, decides on the next point in the set of initial conditions to try. Our goal in this paper is to provide assistance at exactly this last step. Namely, how do we pick the next point in the set of initial conditions? Note that we are essentially looking for a descent direction for the cost function in the set of initial conditions.
Our main contribution, in this paper, is an algorithm that can propose such descent directions. Given a test trajectory starting from a point , the algorithm tries to find some vector such that gets closer to the unsafe set than . We prove that it converges to a local minimum of the robustness function in the set of initial conditions, and demonstrate its advantages within a stochastic falsification algorithm. The results in this paper will enable local descent search for the satisfaction of arbitrary linear temporal logic specifications, not only safety specifications.
2 Problem Formulation
The results in this paper will focus on the model of hybrid automata with affine dynamics. A hybrid automaton is a mathematical model that captures systems that exhibit both discrete and continuous dynamics. In brief, a hybrid automaton is a tuple
where is the state space of the system, is the set of control locations, is the set of control switches, assigns an invariant set to each location, defines the time derivative of the continuous part of the state, is the guard condition that enables a control switch and, finally, is a reset map. Finally, we let to denote the state space of the hybrid automaton .
Formally, the semantics of a hybrid automaton are given in terms of generalized or timed transition systems . For the purposes of this paper, we define a trajectory starting from a point to be a function . In other words, the trajectory points to a pair of control location - continuous state vector for each point in time: , where is the location at time , and is the continuous state at time . We will denote by the sequence of control locations that the trajectory visits (no repetitions). The sequence is finite when we consider a compact time interval and is not Zeno.
Assumptions: In the following, we make a number of assumptions. First, we assume that for each location the system dynamics are affine, i.e., , where and are matrices of appropriate dimensions. Second, we assume that the guards in a location are non-overlapping and that the transitions are taken as soon as possible. Thirdly, we assume that the hybrid automaton is deterministic, i.e., starting from some initial state, there exists a unique trajectory of the automaton. This will permit us to use directly results from . We also make the assumption that the simulation algorithms for hybrid systems are well behaved. That is, we assume that the numerical simulation returns a trajectory that remains close to the actual trajectory on a compact time interval. To avoid a digression into unnecessary technicalities, we will assume that both the set of initial conditions and the unsafe set are included in a single (potentially different) control location.
Let be an unsafe set and let be the distance function to , defined by
where is the projection to the set of locations, is the projection to the continuous state-space and
Definition 1 (Robustness)
Given a compact time interval , we define the robustness of a system trajectory starting at some to be . When is clear from the context, we’ll write .
Our goal in this paper is to find operating conditions for the system which produce trajectories of minimal robustness, as they indicate potentially unsafe operation. This can be seen as a 2-stage problem: first, decide on a sequence of locations to be followed by the trajectory. Second, out of all trajectories following this sequence of locations, find the trajectory of minimal robustness. This paper addresses the second stage. The central step is the solution the following problem:
Given a hybrid automaton , a compact time interval , a set of initial conditions and a point such that , find a vector such that , and .
An efficient solution to Problem 1 may substantially increase the performance of the stochastic falsification algorithms by proposing search directions where the robustness decreases. In summary, our contributions are:
We formulate Problem 1 as a nonlinear optimization problem, which we prove to be differentiable w.r.t. the initial conditions. Thus it is solvable with standard optimizers.
We developed an algorithm, Algorithm 1, to find local minima of the robustness function.
We demonstrate the use of Algorithm 1 in a higher-level stochastic falsification algorithm, and present experimental results to analyze its competitiveness against existing methods.
3 Finding a descent direction
Consider an affine dynamical system in ,
which we assume has a unique solution
where is the initial state of the trajectory
Let be the convex set of bad states, and its closure. Note that even for linear systems, is not necessarily differentiable or convex. Our goal is to find the trajectory of minimum robustness. That is done by a local search over the set of initial conditions.
Given an initial state and a trajectory that starts at , define the time of its closest proximity to , and the point which is closest to the trajectory:
3.1 Partial descent based at the nearest point
Given , choose an approach vector such that is closer to than . Such a vector always exists given that has a positive distance to . Moreover, it is not unique. Thus we have
Define . Then
and is a descent direction, provided that .
It is easy to see that for any and ,
so the new distance is achieved at the same time as the old one. This new distance is an upper bound on the new trajectory’s robustness. In general, the new trajectory’s robustness might be even smaller, and achieved at some other time .
As pointed out earlier, the approach is not unique. The requirement on is that be closer to than . So define the set of points that are closer to than (see Fig. 1):
Then must satisfy . Combined with the requirement that , we get
Any point in the above descent set is a feasible descent direction. As a special case, it is easy to verify that , for any , is a descent direction that leads to 0 robustness. Coupled with the requirement that must be in , it comes
If computing is too hard, we can approximate it with the following : imagine translating along the direction , so it is being drawn closer to , until it meets it. Then we claim that the union of all these translates forms a set of points closer to than :
Let be a convex set, a point outside it, and be the Minkowski sum of and . Then for any in ,
is convex by the properties of Minkowski sums. Let . Then for any , . So translates of boundary points are closer to than .
Now we show that all points in are translates of boundary points. Consider any point in : is in , but is not, so the line crosses for some value : . And, , so by what preceded, .
When , of course, .
We have thus defined 3 possible descent sets: .
The question we address here is: how do we obtain, computationally, points in the descent set , where or ? The following discussion is based on Chapters 8 and 11 of .
Since we’re assuming and to be convex, then the descent set is also convex. Describe with a set of inequalities where the are convex and differentiable, and for convex differentiable (the particular form of the will depend on the descent set at hand). We assume dom = dom .
Given an already simulated trajectory and its time of minimum robustness , we are looking for a feasible such that for some . Thus we want to solve the following feasibility problem
This is a convex program, which can be solved by a Phase I Interior Point method . A non-positive minimum means we found a feasible ; if , then our work is done: we have found an unsafe point. Else, we can’t just stop upon finding a non-positive minimum: we have merely found a new point whose robustness is less than ’s, but not (necessarily) 0. So we iterate: solve t-PDP to get , solve t-PDP to get , and so on, until , a maximum number of iterations is reached, or the problem is unsolvable. If the minimum is positive, this means that for this value of , it is not possible for any trajectory to enter at time .
The program suffers from an arbitrary choice of . One approach is to sample the trajectory at a fixed number of times, and solve (2) for each. This is used in the experiments of this section. A second approach, used in the next section, is to let the optimization itself choose the time, by adding it to the optimization variable. The resulting program is no longer necessarily convex.
3.3 Numerical Experiments
In this section, we present some numerical experiments demonstrating the practical significance of the previous theoretical results.
We consider the verification problem of a transmission line . The goal is to check that the transient behavior of a long transmission line has acceptable overshoot for a wide range of initial conditions. Figure 2 shows a model of the transmission line, which consists of a number of RLC components (R: resistor, L: inductor and C: capacitor) modeling segments of the line. The left side is the sending end and the right side is the receiving end of the transmission line.
The dynamics of the system are given by a linear dynamical system
where is the state vector containing the voltage of the capacitors and the current of the inductors and is the voltage at the sending end. The output of the system is the voltage at the receiving end. Here, , and are matrices of appropriate dimensions. Initially, we assume that the system might be in any operating condition such that . Then, at time the input is set to the value .
The descent algorithm is applied to the test trajectory that starts from and it successfully returns a trajectory that falsifies the system (see Fig. 3).
4 Hybrid systems with affine dynamics
We now turn to the case of hybrid systems with affine dynamics in each location. The objective is still to find a descent direction in , given a simulated trajectory originating at point . Note that since we have assumed that is a singleton set, the problem reduces to finding a descent direction in .
Assumptions. At this point, we make the following assumptions:
b. For every transition , the resets are differentiable functions of their first argument.
c. Conditions 4 and 5 of Theorem III.2 in  are satisfied, namely: for all , there exists a differentiable function such that ; and, for all such that , the Lie derivative . This allows us to have a differentiable transition time of the trajectory starting at the initial point .
d. The sequence of locations enters the location of the unsafe set. This is required for our problem to be well-defined (specifically, for the objective function to have finite values). The task of finding such an is delegated to the higher-level stochastic search algorithm, within which our method is integrated.
4.1 Descent in the Robustness Ellipsoid
Consider a trajectory with positive robustness, with . This is provided by the simulation. Let the initial set be in location and let denote the location of . In order to solve Problem 1, we assume that appears in (see Assumption d above) - otherwise, and the problem as posed here is ill-defined. We search for an initial point (actually ), whose trajectory gets closer to the unsafe set than the current trajectory .
In order to satisfy the constraints of Problem 1, we need to make sure that the new point that we propose generates a trajectory that follows the same sequence of locations as . This constraint can be satisfied using the notion of robust neighborhoods introduced in . In , it is shown that for stable systems and for a given safe initial point , there exists an ‘ellipsoid of robustness’ centered on , such that any trajectory starting in the ellipsoid, remains in a tube around . The tube has the property that all trajectories in it follow the same sequence of locations as . Therefore, we restrict the choice of initial point to , where is the ellipsoid of robustness centered on , with shape matrix . Formally, in , the following result was proven.
Consider a hybrid automaton , a compact time interval , a set of initial conditions and a point . Then, we can compute a number and a bisimulation function , where is a positive semidefinite matrix, such that for any , we have .
(i) In , in the computation of , we also make sure that any point in the robust neighborhood generates a trajectory that does not enter the unsafe set. In this work, we relax this restriction since our goal is to find a point that generates a trajectory that might enter the unsafe set. (ii) In view of Theorem 1, the shape matrix for the ellipsoid is defined as .
We now proceed to pose our search problem as a feasibility problem. Let be the time at which is closest to . We choose as our descent set: recall that it is the set of all points which are closer to than (Def. 1). Therefore, if we can find such that for some , it follows that . To simplify notation, let be the descent set. As before, it is assumed that for differentiable . The search problem becomes:
Given , find and , such that . This is cast as an optimization problem over :
where and .
Note that Problem (3) is specific to a choice of initial point ; this will be important in what follows. In our implementation, the first constraint is specified as bounds to the optimization and so is always satisfied.
Later in this section, we discuss how to solve this optimization problem. For now, we show how solving this problem produces a descent direction for the robustness function. For convenience, for , we define the constraint functions
A point is feasible if it satisfies the constraints in Problem (3). Finally, define the objective function .
The objective function measures the slack in satisfying the constraints: a negative means all constraints are strictly satisfied, and in particular, . Thus, we have a trajectory that enters and, hence, gets strictly closer to . This reasoning is formalized in the following proposition:
Let be a minimum of in program (3). Then .
It is assumed that the optimizer is iterative and that it returns a solution that decreases the objective function. In what follows, for a vector , is the largest entry in .
We first remark that for a given and that satisfy the constraints in (3),
is feasible, and for any feasible . Therefore, we may only consider points with .
Let be the initial point of the optimization. Because is the center of , . And, because , . Thus . Therefore, at the minimum returned by the optimizer, . In particular, , and the new trajectory enters . Therefore, its robustness is no larger than that of the initial trajectory .
We now address how Problem 3 might be solved. Functions , and are differentiable in . It is not clear that , or equivalently, , as a function of , is differentiable. We now show that under some asumptions on the , for trajectories of linear systems, is in fact differentiable in both and , over an appropriate range of . This implies differentiability in . Therefore, standard gradient-based optimizers can be used to solve Problem 3.
For the remainder of this section, we will re-write as to emphasize the dependence on the initial point . will denote the point, at time , on the trajectory starting at , and evolving according to the dynamics of location . When appearing inside location-specific trajectories such as , the time variable will be denoted by the greek letter to indicate relative time: that is, time measured from the moment the trajectory entered , not from datum . (without superscript) will denote the hybrid trajectory, traversing one or more locations. We will also drop the from , and write it simply as .
We first prove differentiability in . Therefore, unless explicitly stated otherwise, the term ‘differentiable’ will mean ‘differentiable in ’. Start by noting that is a composite function of . Since is differentiable, it is sufficient to prove that is differentiable. The hybrid trajectory is itself the result of composing the dynamics from the visited locations .
Recall that is the ellipsoid of robustness centered at . As shown by Julius et al., the following times are well-defined:
Definition 2 (Transition times)
Given , let .
is the time at which trajectory transitions from into through guard .
is the maximal time for which the image of under the hybrid dynamics is contained in :
In other words, is the time at which occurs the first -to- transition of a point in .
is the minimal time for which the image of under the hybrid dynamics is contained in :
In other words, is the time at which occurs the last -to- transition of a point in .
For a given point , () is the absolute (relative) transition time of trajectory from into through guard . Thus, for example, and , with . When the transition is clear from context, we will simply write ().
We will first show differentiability of a trajectory that visits only 2 locations and :
We first present a simple 1D example to illustrate the definitions and the idea of the proof. Consider the hybrid system with three locations
where for , and the flow is defined by
The guards are and . is the identity map, so there are no resets. The initial set is . The solutions in the individual locations are then
We can solve, in this simple case, for : . Similarly for : .
We first show differentiability of the trajectory over locations 0 and 1. We then do the same for a trajectory over locations 1 and 2. Then we stitch the two together and show differentiability over 3 locations. For locations 0 and 1: .
Moving on the trajectory over locations 1 and 2, the procedure is the same: from an initial point , for a fixed (relative time) : .
Finally we stitch up the 2 portions of the trajectory: , . . Since .
We now prove the general case.
Let , and fix . Consider the hybrid trajectory over 2 locations in Eq.(5). If Assumptions a-d are satisfied, then is differentiable at .
In what follows, .
Terms 1 and 2 are clearly differentiable in . For term3, write so term3 = . is differentiable by the Fundamental Theorem of Calculus and its derivative is . As a consequence of Assumption c, itself is differentiable in (Lemma III.3 in ), and the chain rule allows us to conclude that term3 is differentiable in . Thus is differentiable over . Since is differentiable by Assumption b, then is differentiable over . Note that is open and is continuous, so is open. Since is continuous, then is open. Next,
Using the same argument as above, terms 4, 5 and 6 are differentiable in . In conclusion, is differentiable at over , and this ends the proof.
The following proposition generalizes Prop. 3 to trajectories over more than 2 locations.
Fix , and consider the hybrid trajectory over locations. Then is differentiable at for all .
We argue by induction over the number of locations . The base case is true by hypothesis, and the case has been proven in Prop. 3. For and , let be the trajectory over the first locations, so that . By the induction hypothesis, is differentiable at . Then and satisfy the conditions of the case .
Differentiability with respect to time is easily proven:
Let and , that is, a time at which the trajectory is in the last location. Consider the hybrid trajectory over locations. Then is differentiable in over .
. The location-specific trajectories are solutions of differential equations involving at least the first time derivative. Therefore, they are smooth over . This implies differentibility of the hybrid trajectory over the same interval. At , the trajectory is only left-differentiable, since it’s undefined from the right.
The following result is now a trivial application of the chain rule to :
Let , . If is differentiable for all , then is differentiable in over .
We choose Sequential Quadratic Programming (SQP), as a good general-purpose optimizer to solve Problem 3. SQP is a Q-quadratically convergent iterative algorithm. At each iterate, is computed by simulating the system at . This is the main computational bottleneck of this method, and will be discussed in more detail in the Experiments section.
4.2 Convergence to a local minimum
Solving Problem (3), for a given , produces a descent direction for the robustness function. However, one can produce examples where a local minimum of is not a local minimum of the robustness function . This section derives conditions under which repeated solution of Problem (3) yields a local minimum of the robustness function.
For , let , and let be the time when is closest to . Let be the descent set for this trajectory. For each , one can setup the optimization Problem (3) with , and initial point ; this problem is denoted by Prob3. (Recall from the proof of Proposition 2 that at the initial point of the optimization problem). Finally, let be the minimum obtained by solving Prob3.
Algorithm 1 describes how to setup a sequence of optimization problems that leads to a local minimum of . It is called Robustness Ellipsoid Descent, or RED for short.
Algorithm 1 (RED) terminates